ex24ch2.m

these codes are for solving OED with matlab

function ex24ch2
% y(1) = y, y(2) = v, y(3) = phi.
y0 = [0; 0.5; 9.7456];
opts = odeset('Events',@events);
[x,y] = ode45(@odes,[0 10],y0,opts);
plot(x,y(:,1))
title(['Trajectory of cannon ball with \phi(0) = ',num2str(y0(3)),'.'])
fprintf('With phi(0) = %g, the range is %g.\n',y0(3),x(end))
%======================================================================
function dydx = odes(x,y)
% y(1) = y, y(2) = v, y(3) = phi.
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 [value,isterminal,direction] = events(x,y)
value = y(1);
isterminal = 1;
direction = 0;