ex16ch3.m

these codes are for solving OED with matlab

function ex16ch3
% y(1) = y, y(2) = v, y(3) = phi.
solinit = bvpinit(linspace(0,5,5),[0; 0.5; 0.1]);
sol = bvp4c(@odes,@bcs,solinit);
xint = linspace(0,5);
Sxint1 = deval(sol,xint);
solinit = bvpinit(linspace(0,5,5),[0; 0.5; 10]);
sol = bvp4c(@odes,@bcs,solinit);
Sxint2 = deval(sol,xint);
plot(xint,Sxint1(1,:),'b',xint,Sxint2(1,:),'r')
%===========================================================
function dydx = odes(x,y)
g = 0.032;
nu = 0.02;
dydx = zeros(3,1);
dydx(1) = tan(y(3));
dydx(2) = - (g*sin(y(3)) + nu*y(2)^2)/(y(2)*cos(y(3)));
dydx(3) = - g/y(2)^2;

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