ex30ch2.m

来自「these codes are for solving OED with mat」· M 代码 · 共 16 行

function ex30ch2
global epsilon kappa
epsilon = 1/98;
kappa = 3;
options = odeset('Mass',diag([1,epsilon]),'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);