Write a computer program to solve a set of 3 decay equations (N_A N_B N_C, N_C is stable) numerically and analytically. Use a forward difference approximation (explicit scheme) for the numerical solution. The main program should call functions/subroutines to do the following:

Read input file with the input data:

decay constants: l_A, l_B

initial conditions: N_A0, N_B0, N_C0

numerical parameters: Dt, t_final

Perform the numerical solution of the equations for 0 < t < t_final.

Perform the analytical solution of the equations for 0 < t < t_final.

Write results to an output file that can later be used to plot the results. The output file should also contain the parameters read from the input file.

Submit a brief report with your results. See separate instructions on the report format.

Part 1, theory

Show differential equations you are solving (-10%).

Show analytic solution of the differential equations (-10%).

Show complete derivation of numerical solution (-10%).

Show complete derivation for the time of maximum N_B (-10%).

Part 2, radioactive decay chain

Half-life: t_1/2A = 1.1 h , t_1/2B = 9.2 h, t_1/2C = stable

Initial conditions: N_A0 = 100, N_B0 = 0, N_C0 = 0

Solution time: t_final = 50 h

For the numerical solution, start with Dt = 1 h and keep reducing it by the factor of two until you get a reliable solution (solution does not change significantly with deceasing Dt). Show the following results:

Plot numerical N_B(t) vs. time for 3 different values of Dt (coarse, medium, fine), all of them on the same graph. Add analytical solution on the same graph (-10%).

Plot numerical N_A(t), N_B(t), N_C(t) and N_A(t) + N_B(t) + N_C(t) as a function of time, all on the same graph, use Dt that gives reliable solution (-10%).

Using numerical solution, plot time of maximum N_B vs. 1/Dt for several different Dt. Use analytical solution to determine time of maximum N_B, add that value to the graph (-10%).

Note: numerical solution with large Dt might be unstable (oscillations between time-steps). If so, do not report such results. Continue reducing Dt until you get physically realistic (smooth) solution (-10%).