% demo to illustrate some points about stiffness
% we solve the IVP 
%   y' = -100y + 10,    y(0) = 1, 0 <= t <= 10.

% exact solution: y(t) = 1/10 + 9/10 e^{-100t}

clear all
close all

% we begin by plotting the exact solution
t0 = 0; tf= 10; nIntervals = 1000;
tSpan=linspace(t0,tf,nIntervals+1);
yExact = 1/10 + 9/10*exp(-100*tSpan);
plot(tSpan,yExact,'r-')
ylabel('y(t)')
xlabel('t')
hold on
pause
% from here we can see the existence of a 
% boundary layer or initial transient in the solution

% let's see what ode45 does 
options = odeset('Stats','on');
tic,ode45('mildlyStiffODE',[t0,tf],[1],options), toc
% let's zoom in on the solution
axis([0 1 0.09 0.11])
pause
% we see a fair number of step rejections
% and a 'bumpy' solution
% this is typical of the process:
% problem looks easy, take a larger step, error gets too large, reject;
% repeat until integration finished

% how about ode15s?
hold off
tic,ode15s('mildlyStiffODE',tSpan,[1],options); toc
axis([0 1 0.09 0.11])
% we see fewer steps, only 1 rejection, 
% a shorter (but comparable) integration time, 
% and solution is not nearly as bumpy 