ex19ch2.m

these codes are for solving OED with matlab

function ex19ch2
global epsilon kappa
epsilon = 1/98;
kappa = 3;
options = odeset('Stats','on');
[t,y] = ode15s(@f,[0,240],[1;0],options);
semilogx(t,y(:,2))
axis([0.01 100 0 1])
title('Reduced oxygen concentration y(t) when \epsilon = 1/98.');
%===========================================================
function dydt = f(t,y)
global epsilon kappa
dydt = [0; 0];
dydt(1) = - y(1) - y(1)*y(2) + epsilon*kappa*y(2);
dydt(2) =  (y(1) - y(1)*y(2) - epsilon*kappa*y(2))/epsilon;