ex18ch3.m

these codes are for solving OED with matlab

function ex18ch3
global P_B 
P_B = 0.001;
L = 10;
solinit = bvpinit(linspace(0,L,10),[1;1;1;1]);
sol = bvp4c(@odes,@bcs,solinit);
x = sol.y(3,:);
y = sol.y(4,:);
plot(x,y)
axis([0 11 0 1.1])
%=================================================
function dPHIds = odes(s,PHI)
global P_B
dPHIds = [  PHI(2); -P_B*cos(PHI(1)); ...
            cos(PHI(1)); -sin(PHI(1)) ];

function res = bcs(PHIa,PHIb)
res = [ PHIa(1); PHIa(3); PHIa(4)-1; PHIb(2)];