% program to numerically solve the two-element radioactive
% decay problem of equation 1.8


n0=1000;                % initial conditions!
m0=0;

tm=1000;                % e-folding times
tn=100;


totaltime=10000;        % total time over which we want to see the solution

dt=1;                   % delta time, small time step we are using

num=totaltime/dt;       % total number of time steps

m=zeros(num,1);         % set up arrays
n=zeros(num,1);
t=zeros(num,1);		% time


m(1)=m0;			% set up initial conditions
n(1)=n0;
t(1)=dt;

%now we're going to loop and solve the equation

for i=2:num
  n(i)=n(i-1)-n(i-1)/tn*dt;
  m(i)=m(i-1)-m(i-1)/tm*dt+n(i-1)/tn*dt;
  t(i)=i*dt;
end

plot(t,m,'b')		% plot numerical solution in blue
xlabel('time')

% now let's figure the exact solution
mexact=n0*tm/(tn-tm)*(exp(-t/tn)-exp(-t/tm));
hold on
plot(t,mexact,'r--')	% plot exact solution in red dashes

% one good way to compare the numeric and exact solutions
% is to see by what fraction does the numeric solution
% deviate from the analytic one

figure			% open a new plot window in MATLAB
fraction=(m-mexact) ./ mexact;    % calculation deviations
plot(t,fraction)
xlabel('time')
ylabel('fractional deviation')
axis([0 10000 -.01 .01])	% adjust the axes to show relevent parts