ex07ch3.m

these codes are for solving OED with matlab

function ex07ch3
% h = y(1), h' = y(2), f = y(3), f' = y(4)
global beta lambda E h0
beta = 10;
lambda = 10;
E = 1;
h0 = 1;
hold on
for infinity = 2:4
    if infinity == 2
        solinit = bvpinit(linspace(0,infinity,5),[1;0;1;0]);
    else
 		solinit = bvpinit(sol,[0,infinity]);
    end
	sol = bvp4c(@odes,@bcs,solinit);
	plot(sol.x,sol.y(3,:),sol.x(end),sol.y(3,end),'ro');
	axis([0 4 0 1]);
	drawnow
end
plot(0.1,0.7295,'k*')
hold off
%===========================================================
function dydz = odes(z,y)
global beta lambda E h0
g = E - y(1)/lambda;
dydz =  [ y(2); lambda*beta*g*y(1); y(4); beta*g*y(3) ];

function res = bcs(ya,yb)
global beta lambda E h0
res = [ ya(1)-h0; ya(3)-1; yb(1); yb(3) ];