% script to run & compare euler's method to runge-kutta 2nd order
% (mid-point method) using exponential decay as a model

% establish initial conditions & parameters
N_o = 1; % initial # of particles
tau = 1/log(2); % e-folding time of decay (assumed to be positive)
dt = 0.1; % step size in time
t_o = 0; % start time
t_f = 2; % end time

% preallocate results
t = (t_o:dt:t_f)';
nPts = length(t);
N = zeros(nPts,1);
N_euler = N;
N_midPt = N;
N_rk4 = N;

% first step is initial condition
N(1) = N_o;
N_euler(1) = N_o;
N_midPt(1) = N_o;
N_rk4(1) = N_o;

% function handle
dydx = @(t,N) expDec(t,N,tau);

% loop through integration for ode methods
for k = 2:length(t)
    N(k) = N_o*exp(t(k)/(-tau));  % theoretical
    N_euler(k) = odeStep_euler(t(k-1),N_euler(k-1),dt,dydx);
    N_midPt(k) = odeStep_midPt(t(k-1),N_midPt(k-1),dt,dydx);
    N_rk4(k) = odeStep_rk4(t(k-1),N_rk4(k-1),dt,dydx);
end

% get exact results on dense mesh for plot line
te = linspace(t_o,t_f,1e3)';
Ne = N_o*exp(te/(-tau));  % theoretical

% plot results
subplot(311)
plot(te,Ne,'b-',...
    t,N_euler,'ro',...
    t,N_midPt,'g+',...
    t,N_rk4,'kv')
xlabel('t')
ylabel('N(t)')
legend('exact','euler','rk2','rk4')
grid

% plot residuals
subplot(312)
plot(t,(N_euler-N)./N,'ro',...
    t,(N_midPt-N)./N,'g+',...
    t,(N_rk4-N)./N,'kv')
xlabel('t')
ylabel('(N(t)-N_{th}(t))/N_{th}(t)')

% plot log residuals
subplot(313)
semilogy(t,abs(N_euler-N)./N,'ro',...
    t,abs(N_midPt-N)./N,'g+',...
    t,abs(N_rk4-N)./N,'kv')
xlabel('t')
ylabel('(N(t)-N_{th}(t))/N_{th}(t)')


return