ex24ch3.m

these codes are for solving OED with matlab

function ex24ch3
% y(1) = G, y(2) = F, y(3) = F'
fprintf('\n   F''(0)   infinity\n')
clf reset
hold on
for infinity = 2:4
    if infinity == 2
        solinit = bvpinit(linspace(0,infinity,5),[0; 1; 0]);
    else
 		solinit = bvpinit(sol,[0,infinity]);
    end
	sol = bvp4c(@odes,@bcs,solinit);
	plot(sol.x,sol.y(1:2,:),sol.x(end),sol.y(1:2,end),'ro');
	axis([0 4 0 1.5]);
	drawnow
    fprintf('  %g      %i\n',sol.y(3,1),infinity)
end
hold off
%===========================================================
function dydx = odes(x,y)
dydx = [ x*y(3); y(3); 0];
dydx(3) = - 2*( (x*y(2) - y(1))*y(3) + 1 - y(2)^2 );

function res = bcs(ya,yb)
res = [ya(1); ya(2); yb(2)-1];