% This script illustrates stiffness for the (nonlinear) burning match
% problem y'(t) = y(t)^2 - y(t)^3, y(0) = delta, t in [0, 2/delta],
% where y(t) is the radius of the ball of flame.
%
% Physical interpretation: When one lights a match, the ball of flame
% grows rapidly until it reaches a critical size. Then it remains at
% that size because the amount of oxygen being consumed by combustion
% in the interior of the ball (in a volume \propto y^3) balances the
% amount of oxygen available through the surface (\propto y^2). The
% critical parameter is the initial radius delta: the solution grows
% rapidly until t=1/delta, then a quick transition occurs to the
% steady state, y=1.
%
% The exact solution is y(t) = 1 / (W(ae^{a-t}) + 1), where a =
% 1/delta - 1 and W(z) is the Lambert W function (solution of W(z)
% e^{W(z)} = z).
% 
% Things to note:
%
% - the problem is initially not stiff, but becomes stiff when the
%   solution approaches the steady state
% - Euler's Method breaks down for delta small enough
% - explicit solvers with adaptive stepsize control become slower
%   and slower as delta decreases
% - implicit solvers do not need small stepsizes after the sharp 
%   transition to steady state
clear all
close all
f = inline('y^2-y^3','t','y');
t0=0;
tol = 1e-4;
opts = odeset('RelTol',tol);
dvalues = [0.01 0.001 0.0001];
% using Euler with a fixed number of steps, for different deltas
for delta=dvalues
    tmax=2/delta;
    y0=delta;
    N=1000;
    h=tmax/N;
    [t,y]=EulerSystem(t0,y0,f,h,N);
    plot(t,y,'o-');
    title(sprintf('Euler with delta=%f (1000 steps)',delta)) 
    pause
end
% using non-stiff explicit solvers (here, embedded 4th-5th order
% explicit Runge Kutta, usually called "Runge-Kutta-Fehlberg", with
% adaptive stepsize control), for different deltas
for delta=dvalues
    tmax=2/delta;
    y0=delta;    % initial condition
    ode45(f,[t0 tmax],y0,opts);
    title(sprintf('ode45 (embedded explicit RK 4-5) with delta=%f',delta)) 
    pause
end
% using stiff implicit solver ode23tb (first stage trapezoidal rule,
% second stage 2nd order BDF, with adaptive timestep control), for
% different deltas
for delta=dvalues
    tmax=2/delta;
    y0=delta;    % initial condition
    ode23tb(f,[t0 tmax],y0,opts);
    title(sprintf('ode23tb (trapezoidal + BDF 2) with delta=%f',delta)) 
    pause
end
% using stiff implicit solver ode23s (Rosenbrock method, directly
% related to implicit embedded 2nd order 3rd order Gauss-Legrendre
% Runge Kutta), for different deltas
for delta=dvalues
    tmax=2/delta;
    y0=delta;    % initial condition
    ode23s(f,[t0 tmax],y0,opts);
    title(sprintf('ode23s (Rosenbrock ~ implicit RK 2-3) with delta=%f',delta)) 
    pause
end
% using stiff implicit solver ode15s ("Gear method" = BDF with
% adaptive order and adaptive stepsize; monitors the largest and
% smallest eigenvalue of the Jacobian to do adaptivity), for different
% deltas (PS: BDF order > 2 is not a-stable, but almost)
for delta=dvalues
    tmax=2/delta;
    y0=delta;    % initial condition
    ode15s(f,[t0 tmax],y0,opts);
    title(sprintf('ode15s (Gear = adaptive order BDF) with delta=%f',delta)) 
    pause
end