function ShowFlame(delta,tol)
% Solves dy/dt = y^2 - y^3 , y(0) = delta across [0, 2/delta]
% Compares ode23 and ode23s, with relative error tolerance tol

% Interesting choices:
% delta = .001 tol = .0001
% delta = .001 tol = .00000001
% delta = .0001 tol = .0001
% delta = .0001 tol = .00000001
% delta = .00001 tol = .0001

opts = odeset('RelTol',tol,'OutputFcn',@odeplot);
close all
% Use this to count f-evals...
global fcount
fcount = 0;

% The ode23 solution...
[tVals,yVals] = ode23(@MyF,[0 2/delta],delta,opts);
title(sprintf('ODE23: f-evals = %1d tol = %10.3e',fcount,tol),'fontsize',16)
xlabel(sprintf('delta = %8.6f',delta),'fontsize',14)
set(gca,'fontsize',14)

% Use ode23s...
figure
fcount = 0;
[tVals,yVals] = ode23s(@MyF,[0 2/delta],delta,opts);
title(sprintf('ODE23s: f-evals = %1d tol = %10.3e',fcount,tol),'fontsize',16)
xlabel(sprintf('delta = %8.6f',delta),'fontsize',14)
set(gca,'fontsize',14)
shg

function yp = MyF(t,y)
global fcount
fcount = fcount+1;
yp = y^2 - y^3;