Python odeint solver for 2nd Newtons law

Question

I'm just starting to understand how things in Python work so I encountered a small problem. I have to simulate skydiving in Python. The data is in the code. Problem is that air resistance force is depending on velocity, but I don't know the velocity, I just know the position of parachute. So, when I try to solve the differential equation for initial velocity and position, for position it gives me some random numbers, but I want to get numbers from 0 to 4000 m because that is the ordinary position for jumping off plane. I also have drag coefficient as a function of time, and area as function of time because parachute opens after 60 sec and these variables are not the same anymore. Also, the velocity should, as for my understanding, start to increase when the parachute is closed, decrease when the parachute is opened and then accomplish terminal velocity as the parachute is approaching the ground. I don't know what am I doing wrong and how to solve this problem. Every answer would be helpful, thanks in advance.

Make a simulation of skydiving, F = ma, where m is weight of the parachute, a acceleration and F is a sum of gravitational force Fg and air resistance Fd ( drag force).

Drag force depends on vellocity, area, and drag coefficient. Drag coefficient is different for opened and closed parachute kz if 0 ≤ t ≤ t0 ko if ≥ t0

where kz is drag coefficient for closed parachute and ko for opened parachute. Skydiving is described using these differential equations:

m*(d2*z)/(dt2) = m*(dv/dt) = Fg + Fd

where z is the position of parachute Make animation for skydiving

m1 = 14  parachute weight, kg
m2 = 75  human weught, load, kg
m = m1 + m2 , kg
to = 60  time when the parachute opens, s
tmax = 600  time of gliding and free fall, s
z1 = 2500  position when the parachute opens, the path parachute falls like free fall, m
z2 = 4000  route that the parachute makes all together, m
g = 9.81  gravitation, m/s
rho = 1.2466  air density at 10°C
d1 = 11  parachute radius T-10, m
Ao = 95.033  parachute area when opened, m**2
Az = 1   parachute area when closed, m**2
ko = 1.75
kz = 1.

t = np.linspace(0, 600, 101)

# z0 = np.ones_like(t)
# z0[np.logical_and(0. <= t, t <= 60)] = np.linspace(0, 2500, 11)  
# z0[np.logical_and(60 < t, t <= 600)] = np.linspace(2501, 4000, 90)
# def z0(t):
#     for i in range(0, 60):
#         z0 = z1
#         return z1
#         for j in range(61, 600):
#             z0 = z2
#             return z2

k = np.ones_like(t)
k[np.logical_and(0. <= t, t <= 60)] = kz  
k[np.logical_and( 60 < t, t <= 600)] = ko

A = np.ones_like(t)
A[np.logical_and(0. <= t, t <= 60)] = Az  
A[np.logical_and( 60 < t, t <= 600)] = Ao

Fg = m*g # N

z0 = 0 # m
v0 = 50. # , m/s

Fd = (rho*k[0]*A[0]*(v0**2)) / 2

U0 = np.array ([z0(0), v0])

def D(U, t):
    z, v = U
    Dz = v 
    Dv = (1/m) * (Fg - Fd) 
    return np.array ([Dz, Dv])

R = odeint(D, U0, t)
Z = R[:20, 0]
V = R[:20, 1]

plt.subplot(2, 1, 1)
plt.plot(t, R[:, 0], label='Path', color='r')
plt.subplot(2, 1, 2)
plt.plot(t, R[:, 1], label='Velocity', color='b')
plt.xlim((0, 600))
plt.ylim=(0, 4000)
plt.show()
plt.legend()
plt.grid()