import numpy as np

m, g, k = 2.0, 9.81, 0.0567           # massa (kg), zwaarteveld (N/kg), k = ½ · ρ · cd · A (kg/m)


def f(y, t):                          # y = [gevallen afstand x (m), snelheid v (m/s)], naar beneden positief
    x, v = y
    return np.array([v, g - k / m * v * abs(v)])


y, t, h = np.array([0.0, 0.0]), 0.0, 0.1
while y[0] < 40.0:                    # tot het pakket 40 m gevallen is
    y = y + h * f(y, t)               # dezelfde Euler-stap, nu met een vector
    t = t + h
print(f"na {t:.1f} s: {y[0]:.1f} m gevallen, v = {y[1]:.1f} m/s")