% Test of several solvers against a stiff ODE (stiffeq.m).
% CCK 2009-Mar-01
% CCK 2014-Mar-01 Slightly revised.

% Initial conditions:
ode = @stiffeq;
q0 = [1, 0];
tfinal = 10;
Nsteps = 60;  % How many steps to save.
t = tfinal * 10.^( 6 .* ((0:Nsteps)/Nsteps - 1) ); 
   % Series of times to be saved. These are evenly spaced on a log(t) axis!
t0 = t(1);

% Analytic solution:
a = 1000;
ep = 0.001;
u0 = q0(1) + q0(2);
v0 = q0(1) - q0(2);
u = u0 * exp(-(t-t0));
v = v0 * exp(-2001*(t-t0));
x = (u + v)/2;
y = (u - v)/2;

display(" ") % Print a blank line to the terminal.

% Integrate with fixed-step rk23
stiffeq(0,0,1); % Reset the count.
stepper = @rk2;
order = 2;
solver = @(ode, x0, t0, h) twostep(ode, x0, t0, h, stepper, order);
time_begin_run = cputime;
qrk23 = ode_fixed2(solver, @stiffeq, q0, t);
time_seconds_rk23 = cputime - time_begin_run
number_calls_rk23 = stiffeq(0,0,1)
max_error_rk23 = max( abs( [qrk23(:,1)-x'; qrk23(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.

% Integrate with adaptive.m (RK23)
stiffeq(0,0,1); % Reset the count.
stepper = @rk2;
order = 2;
tolerance = [1e-5, 1e-5];
max_rec_depth = 20;
time_begin_run = cputime;
qrk23a = adaptive_stepper2(ode, q0, t, stepper, order,...
    tolerance, max_rec_depth);
time_seconds_rk23_adaptive = cputime - time_begin_run
number_calls_rk23_adaptive = stiffeq(0,0,1)
max_error_rk23_adaptive = max( abs( [qrk23a(:,1)-x'; qrk23a(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.

% Integrate with fixed step, backward_euler
stiffeq(0,0,1); % Reset the count.
backward_euler(0,0,0,0, @stiffeq_jacobian); 
   % set the Jacobian. This only needs to be done once.
solver = @(ode, x0, t0, h) backward_euler(ode, x0, t0, h);
time_begin_run = cputime;
qbe = ode_fixed2(solver, @stiffeq, q0, t);
time_seconds_be = cputime - time_begin_run
number_calls_be = stiffeq(0,0,1)
max_error_be = max( abs( [qbe(:,1)-x'; qbe(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.

% Integrate with adaptive step, doubled backward_euler
% (implicit RK12).
stiffeq(0,0,1); % Reset the count.
stepper = @backward_euler;
order = 1;
tolerance = [1e-4, 1e-4];
max_rec_depth = 20;
time_begin_run = cputime;
qirk12a = adaptive_stepper2(ode, q0, t, stepper, order,...
    tolerance, max_rec_depth);
time_seconds_irk12_adaptive = cputime - time_begin_run
number_calls_irk12_adaptive = stiffeq(0,0,1)
max_error_irk12_adaptive = max( abs( [qirk12a(:,1)-x'; qirk12a(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.


% Integrate with LSODE, explicit method (Adams)
stiffeq(0,0,1); % Reset the count.
lsode_options("integration method", "adams")
time_begin_run = cputime;
qadams = lsode(ode, q0, t);
time_seconds_adams = cputime - time_begin_run
number_calls_adams = stiffeq(0,0,1)
max_error_adams = max( abs( [qadams(:,1)-x'; qadams(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.


% Integrate with LSODE, implicit method (BDF)
stiffeq(0,0,1); % Reset the count.
lsode_options("integration method", "bdf")
time_begin_run = cputime;
qbdf = lsode(ode, q0, t);
time_seconds_bdf = cputime - time_begin_run
number_calls_bdf = stiffeq(0,0,1)
max_error_bdf = max( abs( [qbdf(:,1)-x'; qbdf(:,2)-y'] ) )
display(" ") % Print a blank line to the terminal.
% QUESTION: Does BDF need the Jacobian? If so, how does BDF get the Jacobian?

% PLOTS
figure(20)
semilogx(t, x, 'r-',...
         t, y, 'g-',...
         t,  qrk23(:,1), 'b-+',...
         t,  qrk23(:,2), 'm-+',...
         t, qrk23a(:,1), 'rx',...
         t, qrk23a(:,2), 'gx')
legend('x (analytic)', 'y (analytic)',...
       'x (rk23, fixed)', 'y (rk23, fixed)',...
       'x (rk23, adaptive)', 'y (rk23, adaptive)')
xlabel('t','fontsize',18)
axis([0.00001,10,0,1])
title('stiff\_test1.m: CCK explicit solvers','fontsize',18)

figure(21)
semilogx(t, x, 'r-',...
         t, y, 'g-',...
         t,     qbe(:,1), 'r+',...
         t,     qbe(:,2), 'g+',...
         t, qirk12a(:,1), 'rx',...
         t, qirk12a(:,2), 'gx')
legend('x (analytic)', 'y (analytic)',...
       'x (backward Euler, fixed)', 'y (backward Euler, fixed)',...
       'x (implicit RK12, adaptive)', 'y (implicit RK12, adaptive)')
xlabel('t','fontsize',18)
axis([0.00001,10,0,1])
title('stiff\_test1.m: CCK implicit solvers','fontsize',18)

figure(22)
semilogx(t, x, 'r-',...
         t, y, 'g-',...
         t, qadams(:,1), 'r+',...
         t, qadams(:,2), 'g+',...
         t,   qbdf(:,1), 'rx',...
         t,   qbdf(:,2), 'gx')
legend('x (analytic)', 'y (analytic)',...
       'x (LSODE, Adams)', 'y (LSODE, Adams)',...
       'x (LSODE, BDF)', 'y (LSODE, BDF)')
xlabel('t','fontsize',18)
axis([0.00001,10,0,1])
title('stiff\_test1.m: LSODE','fontsize',18)
