human
I am working on the file "02_numpy_md/recall_last_session.ipynb".
The first task is:
As you can see there is a difference between the numerical and analytical solution.
If you decrease the time step by a factor of 10 to:
python
time_step = 0.004
and rerun the numerical solution and then analytical solution you will see that for the given height the curves look the same
but we still have a error between them which increases with increasing time / height
Numerical solution second order
A more accurate solution can be provided by using a numerical integrator of higher order.
Taylor series
A Taylor series is a representation of a function as an infinite sum of terms that are calculated from the values of the function's derivatives at a single point. 
<div style="background:lightblue">
Note: see math script of Luik for more
</div>

One can use Taylor series to numericaly approximate functions by using a finite number of terms.
Here, we give as an example a 1-d Taylor series which was terminated after the third order term:
$
\begin{align}
x(t + \Delta t) &= x(t) + \dot{x}(t) \Delta t + \frac{1}{2} \ddot{x}(t) (\Delta t)^2 + \frac{1}{3!} \dddot{x}(t) (\Delta t)^3 + \mathcal{O}(\Delta t ^ 4) \
\end{align}
$
If we now take the second order form and use it for our example we get:
$
\begin{align}
x(t + \Delta t) &= x(t) + \dot{x}(t) \Delta t + \frac{1}{2} \ddot{x}(t) (\Delta t)^2 + \mathcal{O}(\Delta t ^ 3) \
x(t + \Delta t) &= x(t) + v(t) \Delta t + \frac{F}{2 m} \Delta t^2 \
\end{align}
$
Since we have a constant force in both direction we get for the x-direction:
$
\begin{align}
x(t + \Delta t) &= x(t) + v_x(t) \Delta t + \frac{F_x}{2 m} \Delta t^2 \
x(t + \Delta t) &= x(t) + v_x(t) \Delta t &\textrm{with: }F_x=0 \
\end{align}
$
and in y-direction: 
$
\begin{align}
y(t + \Delta t) &= y(t) + v_y(t) \Delta t + \frac{F_y}{2 m} \Delta t^2 \
y(t + \Delta t) &= y(t) + v_y(t) \Delta t - \frac{g}{2} \Delta t^2 
&\textrm{with: }F_y=-g m \
\end{align}
$
This equations are implemented in the code below.
Which will now update the positions with a second order integrator
python
self.position[0] +=...
Can you write Python program for it?
 ... Give me a the Python program you requested?
assistant
#Here is the Python program you requested.

class Ball:
    def __init__(self, start_position, start_velocity, radius):
        self.position = start_position
        self.velocity = start_velocity
        self.radius = radius
    
    def move(self, time_step):
        self.position[0] += self.velocity[0] * time_step
        self.position[1] += self.velocity[1] * time_step - 9.81 / 2.0 * time_step**2

        self.velocity[0] += 0
        self.velocity[1] -= 9.81 * time_step
        if self.position[1] <= 0.0:
            self.velocity[1] *= -1
        
tennis_ball = Ball([0, 1], [0.2, 0], 0.033)
time_step = 0.04
number_of_steps = 1500

trajectory_numerical_2nd = []  # this is an empty list
trajectory_time = [] # also a empty list to store the time
for step in range(number_of_steps):
    trajectory_numerical_2nd.append(tennis_ball.position[:])  # append adds something at the end of a list
    trajectory_time.append(step*time_step)
    if tennis_ball.position[1] <= 0.0: break # stop the calculation here
    tennis_ball.move(time_step)
        
# don't worry about the code below too much it just crates an animation
# of your calculated positions.
from animation import animate
%matplotlib notebook

#ani = animate(trajectory_numerical, time_step*100/0.04)

trajectory_analytical = [] # an empty list to store the analytic results
start_position=[0, 1]
start_velocity=[0.2, 0]
for point_in_time in trajectory_time:
    position_x = start_position[0] + start_velocity[0] * point_in_time
    position_y = start_position[1] - 0.5 * 9.81 * point_in_time**2
    trajectory_analytical.append([position_x,position_y])
    

error = []
error_x_sq=error_y_sq=0
for i in range(len(trajectory_analytical)):
    diff_x=trajectory_analytical[i][0]-trajectory_numerical_2nd[i][0]
    diff_y=trajectory_analytical[i][1]-trajectory_numerical_2nd[i][1]
    error_x_sq+=diff_x**2
    error_y_sq+=diff_y**2
    error.append([error_x_sq,error_y_sq])
    
from animation import plot
%matplotlib inline

plot(trajectory_time, trajectory_numerical_2nd,trajectory_analytical,error)