function [ ] = prob_2_24
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

global g nu;

g = 0.032;
nu = 0.02;

clc; close all; 

options = odeset('Events',@events);

y0 = [ 0 0.5 0.3782 ];
[ts,ys,te,ye,ie] = ode45( @f, [0, 10], y0, options );
figure;
plot( ts, ys(:,1), '-og' ); grid on; 

if( ~isempty(ie) )
  disp( 'Crossing the x-axis at' ); 
  ie, te, ye
end

y0 = [ 0 0.5 9.7456 ];
[ts,ys,te,ye,ie] = ode45( @f, [0, 10], y0, options );
figure;
plot( ts, ys(:,1), '-og' ); grid on; 

if( ~isempty(ie) )
  disp( 'Crossing the x-axis at' ); 
  ie, te, ye
end

return;



function [dydt] = f(t,y)
%
% y is 3x1 = [ y, v, phi ]
%
global g nu;

dydt = [ 
         +tan(y(3));
         -( g*sin(y(3)) + nu*(y(2)^2) )/( y(2)*cos(y(3)) );
         -g/(y(2)^2)
       ];

return;



function [value,isterminal,direction] = events(t,y)

global g nu;

value = y(1);
isterminal = 0;
direction = -1;

return;