function [ ] = prob_2_23
% 
% 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 q d gamma R L;

g = 9.81;
q = 20;
d = 5;
gamma = 1.4;
R = 4;
L = 3.75;

y0 = [ 0 0 ];

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

[ts,ys,te,ye,ie] = ode15s( @f, [0, 100], y0, options );
figure;
plot( ts, ys(:,1), '-og' ); 
xlabel( 'time (sec)' ); ylabel( 'position of cork (cm)' );
grid on;

if( ~isempty(ie) )
  disp( 'Position of the cork at the neck is' ); 
  ye(1)
  disp( 'Velocity of the cork at the neck is' ); 
  ye(2)
end

return;



function [dydt] = f(t,y)
%
% y is 2x1
%
global g q d gamma R L;

dydt = [ y(2); 
         g*(1+q)*( ( 1 + y(1)/d )^(-gamma) + ( R*t/100 ) - 1 + q*y(1)/(L*(1+q)) )
       ];

return;



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

global g q d gamma R L;

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

return;