ex28ch2.m

these codes are for solving OED with matlab

function ex28ch2
% The ball is at (x(t),y(t)).
% Here y(1) = x(t), y(2) = x'(t), y(3) = y(t), y(4) = y'(t).
clc
% Set initial height and coefficient of restitution:
y0 = [0; 0; 2; 0]; 
k = 0.70;
k = 0.35;
% Plot the initial configuration:
fill([0 0 1],[0 1 0],[0.8 0.8 0.8]); 
axis([-0.1 1.0 0 y0(3)]) 
hold on
plot(0,y0(3),'ro');  % Mark the initial point.
options = odeset('Events',@events);
% Accumulate the path of the ball in xplot,yplot.
xplot = []; 
yplot = []; 
tstar = 0; 
while 1
    tspan = linspace(tstar,tstar+1,20);
    [t,y,te,ye,ie] = ode23(@odes,tspan,y0,options);

    % Accumulate the path.
    xplot = [xplot; y(:,1)];
    yplot = [yplot; y(:,3)];
    if isempty(ie)          % Extend the interval.
        tstar = t(end);
        y0 = y(end,:);
    elseif ie(end) == 1     % Ball bounced.
        plot(ye(end,1),ye(end,3),'r*'); % Mark the bounce point.
        if (te(end) - tstar) < 0.01*tstar
            fprintf('Bounces accumulated at x = %g.\n',ye(end,1))
            break;
        end
        if abs(ye(end,1) - 1) < 0.01
            break;
        end
        tstar = te(end);
        y0 = [ye(end,1); -k*ye(end,4); ye(end,3); k*ye(end,2)];
    elseif ie(end) == 2     % Ball hit the wall.
        plot(ye(end,1),ye(end,3),'w');
        if ye(end,3) < 0.01  % Ball is in the corner.
            fprintf(' The ball is in the corner.')
            break;
        end
        % Reset the initial conditions and continue integrating
        % as it bounces back onto the ramp.
        tstar = te(end);
        y0 = [ye(end,1); -k*ye(end,2); ye(end,3); ye(end,4)];
    end
end
% Plot the solution.
plot(xplot,yplot);
title(['Ball on ramp with k = ', num2str(k)])

function dydt = odes(t,y)
% g = 9.81.
dydt = [y(2); 0; y(4); - 9.81];

function [value,isterminal,direction] = events(t,y) 
value = [ y(3) - (1 - y(1))
              y(1) - 1     ];
isterminal = [1; 1]; 
direction  = [-1; 0];