import numpy as np
import matplotlib.pyplot as plt

y0 = 10.0
M = .015
R = 0.4
g = 9.81
eta = 1.002e-3
alpha = 6*np.pi*R*eta

t_final = 10.0
n = 100
t_vec = np.linspace(0.0, t_final, n+1)
h = t_vec[1]-t_vec[0]

y = np.zeros(n+1)
v = np.zeros(n+1)
y[0] = y0
v[0] = 0.0

for j in range(n):
    F = -M*g - alpha * v[j]
    y_new = y[j] + h * v[j]
    v_new = v[j] + h * F / M
    y[j+1] = y_new
    v[j+1] = v_new

plt.figure()
plt.plot(t_vec, y, 'o-', label='y')
plt.plot(t_vec, v, 'o-', label='v')
plt.legend()
plt.show()