function ShowChem(tol)
% Applies several Matlab's stiff ODE solvers to an n=3 chemical
% reaction system.

% Interesting choices: tol = .01, .001 , .0001

% tol is the relative error tolerance passed to the solvers
opts = odeset('RelTol',tol);
close all
global fcount % for f-eval counting
y0 = [1 0 0]'; % initial condition
tSpan = [0 5]; % range of integration

% ode23s
fcount = 0;
[tOut23s,yOut23s] = ode23s(@Chem,tSpan,y0,opts);
% Display the three solutions y1(t), y2(t), and y3(t)...
subplot(3,1,1)
plot(tOut23s,yOut23s(:,1))
title(sprintf('y1, y2, y3 via ode23s f-evals = %1d ',fcount),'fontsize',14)
subplot(3,1,2)
plot(tOut23s,yOut23s(:,2))
subplot(3,1,3)
plot(tOut23s,yOut23s(:,3))

% ode23t
figure
fcount=0;
[tOut23t,yOut23t] = ode23t(@Chem,tSpan,y0,opts);
% Display the three solutions y1(t), y2(t), and y3(t)...
subplot(3,1,1)
plot(tOut23t,yOut23t(:,1))
title(sprintf('y1, y2, y3 via ode23t f-evals = %1d ',fcount),'fontsize',14)
subplot(3,1,2)
plot(tOut23t,yOut23t(:,2))
subplot(3,1,3)
plot(tOut23t,yOut23t(:,3))

% ode23tb
figure
fcount=0;
[tOut23tb,yOut23tb] = ode23tb(@Chem,tSpan,y0,opts);
% Display the three solutions y1(t), y2(t), and y3(t)...
subplot(3,1,1)
plot(tOut23tb,yOut23tb(:,1))
title(sprintf('y1, y2, y3 via ode23tb f-evals = %1d ',fcount),'fontsize',14)
subplot(3,1,2)
plot(tOut23tb,yOut23tb(:,2))
subplot(3,1,3)
plot(tOut23tb,yOut23tb(:,3))

% ode15s
figure
fcount=0;
[tOut15s,yOut15s] = ode15s(@Chem,tSpan,y0,opts);
% Display the three solutions y1(t), y2(t), and y3(t)...
subplot(3,1,1)
plot(tOut15s,yOut15s(:,1))
title(sprintf('y1, y2, y3 via ode15s f-evals = %1d ',fcount),'fontsize',14)
subplot(3,1,2)
plot(tOut15s,yOut15s(:,2))
subplot(3,1,3)
plot(tOut15s,yOut15s(:,3))

% ode23
figure
fcount=0;
[tOut23,yOut23] = ode45(@Chem,tSpan,y0,opts);
% Display the three solutions y1(t), y2(t), and y3(t)...
subplot(3,1,1)
plot(tOut23,yOut23(:,1))
title(sprintf('y1, y2, y3 via ode23 f-evals = %1d ',fcount),'fontsize',14)
subplot(3,1,2)
plot(tOut23,yOut23(:,2))
subplot(3,1,3)
plot(tOut23,yOut23(:,3))

% Compare solutions at tFinal:
fprintf('RelTol = %8.1e\n\n',tol)
fprintf('ode23s %13.10f %13.10f %13.10f\n',yOut23s(length(tOut23s),:))
fprintf('ode23t %13.10f %13.10f %13.10f\n',yOut23t(length(tOut23t),:))
fprintf('ode23tb %13.10f %13.10f %13.10f\n',yOut23tb(length(tOut23tb),:))
fprintf('ode15s %13.10f %13.10f %13.10f\n',yOut15s(length(tOut15s),:))
fprintf('ode23 %13.10f %13.10f %13.10f\n',yOut23(length(tOut23),:))

function yp = Chem(t,y)
global fcount
alpha = .04;
beta = 10000.00;
gamma = 30000000.00;
yp = [ -alpha*y(1) + beta*y(2)*y(3);...
alpha*y(1)-beta*y(2)*y(3) - gamma*y(2)*y(2);...
gamma*y(2)*y(2)];
fcount = fcount+1;