ex30ch3.m

these codes are for solving OED with matlab

function sol = ex30ch3
global beta_0
sol = bvpinit(linspace(0,1,10),ones(3,1));
for beta_0 = [10 100 500 1000 1575]
    sol = bvp4c(@odes,@bcs,sol);
    plot(sol.x,sol.y(1,:),sol.x,100*sol.y(2:3,:))
    title(['Solution with \beta_0 = ',num2str(beta_0),'.'])
    drawnow
%   pause
end
%=================================================
function dydt = odes(t,y)
global beta_0
mu = 0.02;
lambda = 0.0279;
eta = 0.01;
beta = beta_0*(1 + cos(2*pi*t));
dydt = [ mu - beta*y(1)*y(3)
         beta*y(1)*y(3) - y(2)/lambda
         y(2)/lambda - y(3)/eta ];

function res = bcs(ya,yb)
res = yb - ya;