We know the exact equations of motion of the ball, except during the collision. However, by assuming the collisions take place instantaneously, we can solve the equations of motion up to the exact instant of time the ball touches the floor, and then, switch the velocity using the coefficient of restitution, and then, solve the exact analytical equations again until the next exact time the ball hits the floor. In this case, use non uniform 𝑑𝑡. Since you can calculate the time between bounces, always plot 40 points between the times the ball bounces. In this way, you will have an exact solution under the simplifying assumptions of the model. Make a code in MatLab by modifying the code below:

%% simulating a bouncing ball
% This script simulates a ball dropped or thrown in a vertical direction
% under constant gravitational acceleration. The ball will hit the floor
% and it will bounce a multiple number of times based on the coefficient
% of restitution given by epsilon. We have just before/after a collision.
% epsilon = |just after velocity| divided by |just before velocity|
%
% INPUT
% yo = initial height of ball
% vyo = initial velocity of ball in the y-direction
% dt = time resolution
% tMax = total time of the process
% epsilon = coefficient of restitution
%
% PROCESS
% solve differential equation using Euler’s method
% use crude approximation to determine when ball hits floor
%
% OUTPUT
% height versus time graph
% velocity versus time graph
clear all % never not do this
%% set initial condition
to = 0; % initial time on clock (seconds)
yo = 2.5; % meters (m)
vyo = 5; % velocity --> meters/second (m/s)
dt = 0.001; % time resolution --> seconds (s)
tMax = 20; % maximum time range --> (s)
epsilon = 0.7; % coefficient of restitution (unitless)
%% set constants
ay = -9.8; % acceleratin due to gravity m/s^2
%% solve differential equation using Euler’s method
%% Create interpolation data for velocity versus height
heights = [0, 2.5]; % heights at which velocity is known (m)
velocities = [0, vyo]; % velocities at corresponding heights (m/s)
v_interpolation = @(h) interp1(heights, velocities, h); % interpolation function

%% solve differential equation using Euler’s method
tArray = to; % time range as an array (s)
y = yo; % preallocate the y array (height) (m)
vy = vyo; % preallocate the vy array (height) (m)

% ------------------------------- calculate position as a function of time
while tArray(end) < (to+tMax)
% ----------------------------------------------------------- detect floor
t = tArray(end);
vy_next = vy(end) + ay*dt;
y_next = y(end) + vy(end)dt;

if (y_next <= 0) % => hit the floor with dt time resolution
% Linear interpolation for velocity just before collision at y = 0
vy_collision = v_interpolation(y_next);

% Update the velocity and height after collision
vy_next = - vy_collisionepsilon;
y_next = 0;

% Append the interpolated velocities and heights to the arrays
tArray = [tArray, t+dt/2, t+dt];
vy = [vy, vy_collision, vy_next];
y = [y, 0, 0];
else
% Append the next velocity and height to the arrays
tArray = [tArray, t+dt];
vy = [vy, vy_next];
y = [y, y_next];
end
end

%% plot results
figure(1);
clf;
plot(tArray,y,‘b’,‘lineWidth’,1.8);
xlabel(‘time (seconds)’);
ylabel(‘height (meters)’);

figure(2);
plot(tArray,vy,‘b’,‘lineWidth’,1.8);
xlabel(‘time (seconds)’);
ylabel(‘velocity (meters)’);