%* Matab script to calculate exponential decay using 
% 	adaptive Runge-Kutta method
%
% variables set by user:
% lambda:		decay constant
% c0:			initial concentration
% ti:			initial time
% tf:			end time
% myfun:		name of function containing RHS of ODE
%
% results:
% crk:			vector containing concentration at all steps
% trk:			vector containing time at all time steps

clear

lambda=1;					%decay constant
c0=1;						%initial concentration
ti=0;						%initial time
tf=5;						%end time
myfun='f_decay_rk';		%rhight-hand side written for ode45


%* numerical solution

% set some ode45 options
% epsilon:		error tolerance
%epsilon=0.01			
%options=odeset('RelTol',epsilon,'AbsTol',epsilon); 
options=[];				%default values
 
[trk crk]=ode45(myfun, [ti tf], c0,options,lambda);

%* analytic solution
g=c0 * exp(-lambda*trk);


%* cpmpare results
figure(1)
clf
axis([ti tf 0 1]), hold on 
plot(trk, g, 'r-')
plot(trk, crk, 'bo')
legend('analytic solution', 'adaptive Runge-Kutta solution')
title('concentration of a radioactive element over time')
xlabel('time'), ylabel('concentration')