ex09ch3.m

these codes are for solving OED with matlab

function sol = ex09ch3
global d D
left = [1 0.1 0.01];
right = [20 40 60];
hold on
for i = 1:3
    d = left(i);
    D = right(i);
	if i == 1
        solinit = bvpinit(linspace(d,D,5),@guess,-1.6);
	else
		solinit = bvpinit(sol,[d,D]);
	end
	sol = bvp4c(@odes,@bcs,solinit);
	plot(sol.x,sol.y(1,:),sol.x(end),sol.y(1,end),'ro',...
         sol.x(1),sol.y(1,1),'ko');
	axis([0 60 0 1]);
	drawnow
end
hold off
figure
x = [0 sol.x];
y = [1 sol.y(1,:)];
plot(x,y)
axis([0 15 0 1])
fprintf('y''(0) is about %6.4f.\n',sol.parameters);
%===========================================================
function v = guess(x)
if x <= 1
    v = [1; 0];
else
    v = [ 144/x^3; -432/x^4 ];
end

function dydx = odes(x,y,p)
dydx = [ y(2); max(y(1),0)^(3/2) / sqrt(x) ];
 
function res = bcs(ya,yb,p)
global d D
res = zeros(3,1);
res(1:2) = ya - series(d,p);
res(3) =   yb(1) - 144/D^3 ;

function y = series(x,p)
yx = 1 + p*x + (4/3)*x^(3/2) + (2/5)*p*x^(5/2);
ypx = p + 2*x^(1/2) + p*x^(3/2);
y = [yx; ypx];