Numerically Solving 2nd Order System of Differential Equations $g, k_d, k_l$ are constants.
step size = $0.01$ sec = $t$
Intial Conditions:
x,y,z = 0 
$\frac{dx}{dt}=10.45$ 
$\frac{dy}{dt}=5.45$ 
$\frac{dz}{dt}=33.45$ 
$v_0=25$
 A: In python (and similar in other languages) you would encode the ODE function of the right side, for the first order system with $v_x=\dot x$ etc., as
def motioneqns(t,u):
    x,y,z,vx,vy,vz = u
    v = (vx**2+vy**2+vz**2)
    ax = -kd*v*vx + kl*v*vy
    ay = -kd*v*vy - kl*v*vx
    az = -kd*v*vz - g
    return vx,vy,vz,ax,ay,az

and pass it along to any general-purpose solver, which might be one implementing the Euler method, or RK4, or some higher order method with adaptive step size.
If you use implementations similar to those in Does fourth-order Runge-Kutta have an higher accuracy than the second-order one? you need to wrap the result list in a vector type, for instance numpy.array, either in the return line or in the function call. For further inspiration on getting the trajectories as list see https://stackoverflow.com/questions/61333687/how-do-i-solve-a-2nd-order-differential-equation-for-projectile-motion, https://stackoverflow.com/questions/53645649/cannot-get-rk4-to-solve-for-position-of-orbiting-body.
