%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.6
% (c) 08/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
lambda = 0.023051;      % coefficient (decrease) [1/years]
N0 = 20.7;              % initial condition (radioactivity) [Bq/g]
t_max = 200.0;          % time under consideration
dt = 20;                % time step size [years]
n_a = 200;              % number of sample points 
                        % (analytical solution)

%-----------------------------------------------------------------
% Calculation / time integration by forward differences
t_f(1) = 0.0;                      % first value (time)
N_f(1) = N0;                       % initial value (radioactivity)
for i=1:1:t_max/dt
  t_f(i+1) = i*dt;
  N_f(i+1) = (1.0-lambda*dt)*N_f(i);                % Eq. (3.97)
end

% Calculation / time integration by backward differences (implicit)
t_b(1) = 0.0;                      % first value (time)
N_b(1) = N0;                       % initial value (radioactivity)
for i=1:1:t_max/dt
  t_b(i+1) = i*dt;
  N_b(i+1) = N_b(i)/(1.0+lambda*dt);                % Eq. (3.101)
end

%-----------------------------------------------------------------
% Analytical solution
t_a(1) = 0.0;                    % first value (time)
N_a(1) = N0;                     % initial value (radioactivity)
for i=1:1:n_a
  t_a(i+1) = i*t_max/n_a;
  N_a(i+1) = N0*exp(-lambda*t_a(i+1));              % Eq. (3.92)
end

%-----------------------------------------------------------------
% Output / plot
figure(1)
plot(t_f,N_f,'r',t_b,N_b,'b',t_a,N_a,'k')
xlabel('t / years')
ylabel('Radioactivity N [ Bq / g ]')
title('red: forward diff      blue: backward diff     black: analytical')