%* Euler single-step solution for decay series system of ODEs 
% lambda  decay constants, incl. stable daughter
% c    concentration (dependent variable)
% c0   initial conditions on c
% t    time (independent variable)
% ti   initial time
% tf   ending time for calculation
% h    step size
clear

lambda=[1 1 0];
c0=[1 0 0]; 	
ti=0;		
tf=5;		
h=0.05;         

t=[ti:h:tf]';	
N=length(t)-1;	   % number of steps, N+1 = number of nodes

C=zeros(N+1, length(c0));	
C(1,:)=c0;	   % insert initial conditions

for j=1:N
  C(j+1,:) = C(j,:) + h * f_decayseries(lambda, C(j,:));
end


%* information for plotting
ymin=min(min(C));
ymax=max(max(C));

figure(1)
clf
axis([min(t) max(t) ymin ymax]), hold on
plot(t, C(:,1), 'b-')
plot(t, C(:,2), 'r-')
plot(t, C(:,3), 'r-.')
legend('parent', 'daughter 1, radioactive', 'daughter 2, stable')
ylabel('concentration')
xlabel('time')
title('decay series')