Edit: Corrected the integration

This is a simulation I made for a ball falling (the principles can be used everywhere):

import numpy as np
import matplotlib.pyplot as plt


def calc_acceleration(pos, v):
    return np.array([0, 0, -9.8])


def integrate(pos, v, a_old, delta_t):
    pos_new = pos + v * delta_t + 0.5 * a_old * delta_t * delta_t

    a_new = calc_acceleration(pos_new, v)

    v_new = v + 0.5 * (a_old + a_new) * delta_t

    if pos_new[2] < 0.0:
        v_new[2] = -v_new[2]

    return pos_new, v_new, a_new


# Constants
delta_t = 0.01
t_max = 10

# Time evolution
times: np.ndarray = np.arange(0, t_max, delta_t)

# Initial state
v = np.array([15.0, -3.0, 15.0])
pos = np.array([0.0, 0.0, 7.0])
positions = np.zeros([times.size, 3])
a_old = calc_acceleration(pos, v)

# Simulating
for i in range(times.size):
    pos, v, a_old = integrate(pos, v, a_old, delta_t)
    positions[i, :] = pos

# Plotting
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.scatter(positions[:, 0], positions[:, 1], positions[:, 2])
plt.show()

Very simple, you can see a ball falling and getting up again from that starting point, so you just see at least some starting point on how to do this stuff, the basics are not that difficult. But every simulation is like this, just more complicated hahaha.