% Demonstration of a stiff system
% AMATH301, 2015 Winter Term, U Washington
% Instructor: Dr. King-Fai Li

clear all; close all; clc;
% Flame model described in Mathworks.com
F = @(t,y) y.^2 - y.^3 ;
y0 = 0.0001;

% Analytic solution
t = linspace(0,2/y0,1000);
a = 1/y0 - 1;
plot(t,1./(1+lambertw(a*exp(a-t))),'-r','LineWidth',1.2);
axis([t(1) t(end) 0 1.2])
hold on

% Plot as ODE45 calculates
ode45(F,[0 2/y0],y0);

figure(2)
plot(t,1./(1+lambertw(a*exp(a-t))),'-r','LineWidth',1.2);
axis([t(1) t(end) 0 1.2])
hold on

% Plot as ODE23tb calculates
ode23tb(F,[0 2/y0],y0);
