ex20ch3.m

these codes are for solving OED with matlab

function sol = ex20ch3
Da = 0.05;
solinit = bvpinit(linspace(0,1,5),[1; 1]);
sol = bvp4c(@odes,@bcs,solinit,[],Da);
plot(sol.x,sol.y);
fprintf('The initial values are about (%3.1f,%3.1f)\n',...
        sol.y(1,1),sol.y(2,1));
%===========================================================
function dydx = odes(x,y,Da)
dydx = [0; 0];
dydx(1) = Da*(1 - y(1))*exp(y(2)/(1 + 0.05*y(2)));
dydx(2) = 6*dydx(1);
	 
function res = bcs(ya,yb,Da)
res = ya - 0.5*yb;