```function ball(initial_velocity, timestep)
% Program to compute the position and velocity of a ball
%    which is thrown vertically upward from ground level.
%    We make simplifying assumptions:
%          - no air resistance
%          - constant force of gravity
%
% User supplies the initial velocity of the ball and
%    the value of the timestep for numerical integration.
%    We use simple Euler's method to compute the motion.
%
% Print out lines which include mechanical energy during each timestep.
%
% Usage: ball (initial_velocity, timestep)
%
%    where
%             initial_velocity       is initial speed of ball
%                                    vertically upward (m/s)
%
%             timestep               is size of steps in time (s)
%
% Returns: nothing
%
% MWR 11/28/2012

% initialize variables
%    current height above the ground (m)
y = 0;
%    new height above the ground (m)
y_new = 0;
%    current velocity (positive upward) (m/s)
vy = initial_velocity;
%    new velocity (m/s)
vy_new = 0;
%    time since start of simulation (s)
t = 0;
%    acceleration due to gravity (positive up) (m/s^2)
g = -9.8;
%    mass of ball (kg)
m = 1.0;

% compute initial mechanical energy
ke = 0.5*m*vy*vy;
gpe = m*abs(g)*y;
initial_energy = ke + gpe;
% fprintf(1, ' ke %9.3f  gpe %9.3f  tot %9.3f \n', ke, gpe, energy);

% loop over time; continue as long as ball is above ground
while (y >= 0)

%  compute mechanical energy
ke = 0.5*m*vy*vy;
gpe = m*abs(g)*y;
energy = ke + gpe;
frac_error = (energy - initial_energy) / initial_energy;

%  print current time, position velocity, energy, frac error in energy
fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  E %12.6e  f %9.4e \n', t, y, vy, energy, frac_error);

%  use const accel to compute new velocity
vy_new = vy + g*timestep;

%  use current velocity to compute new position
y_new = y + vy*timestep;

%  update position and velocity for next time through loop
vy = vy_new;
y = y_new;

%  update time
t = t + timestep;

% end of loop over time
end

%  print current time, position velocity, energy, frac error in energy
fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  E %12.6e  f %9.4e \n', t, y, vy, energy, frac_error);

% end of program
end
```