ex23ch2.m

these codes are for solving OED with matlab

function ex23ch2
% x = y(1), x' = v = y(2)
options = odeset('Events',@event);
[t,y,te,ye,ie] = ode45(@odes,[0,100],[0;0],options);
if isempty(ie)
    error('Cork did not leave bottle by t = 100.')
end
fprintf('Cork left bottle at t = %g with speed v(t) = %g. \n',...
         te(end),ye(end,2))
%========================================================
function dydt = odes(t,y)
g = 9.81;
q = 20;
d = 5;
gamma = 1.4;
R = 4;
L = 3.75;
dydt = [y(2); 0];
dydt(2) = g*(1+q)*((1 + y(1)/d)^(-gamma) + R*t/100 - 1 + ...
                    q*y(1)/(L*(1+q)));

function [value,isterminal,direction] = event(t,y)
L = 3.75;
value = y(1) - L;
isterminal = 1;
direction = 0;