Most of my knowledge about numerically solving differential equations is long forgotten. Unfortunately I stumbled upon a physics problem where I need to do exactly that.

This is the initial second order equation which I would like to solve for d(t):

$ d =  \frac{ 10 C_B V_{TLC} }{m d'' - C_R (d')² - C_B V_{diver} + mg} - 10   $

following a number of tutorials I rewrote this as a system of two first order equations as

$ d' = v $  
$ v' = d'' = \frac{  C_R v² + \frac{ 10 C_B V_{TLC} }{(d + 10)}  - mg + C_B V_{diver}   }{m}  $

then I tried to solve this in python:

    import numpy as np
    import math
    import scipy as sp
    from scipy.integrate import odeint
    from scipy.integrate import solve_ivp

    # Physical constants
    rho = 1023.6               # kg/m³ density of saline water
    g = 9.807                  # m/s²  gravitational acceleration on earth

    # Assumptions about diver
    V_diver = 0.062            # m³    volume of diver
    V_tlc = 0.006              # m³    total lung capacity
    m = 66                     # Kg    weight of diver
    A = 0.07                   # m²    crossectional area of diver in diving direction              
    C_D = 0.3                  # -     Drag coefficient

    # Derived
    C_B = rho * g              # Buoyency coefficient
    C_R = 0.5 * rho * C_D * A  # Resistive coefficient

    # equation
    def dSdd(d, S):
        d, v = S
        return [
            v,
            C_B * V_diver + 10 * C_B * V_tlc / ( d + 10 ) - m * g + C_R * v**2
        ]

    # initial conditions
    d_0 = 20
    v_0 = 1
    S_0 = [d_0, v_0]

    # time interval
    t = np.linspace(0, 60, 1000)
    
    # solution odeint
    odeint(dSdd, y0 = S_0, t=t, tfirst=True, full_output = 1) 

    # solution solve_ivp
    solve_ivp(dSdd, t_span=(0, max(t)), y0=S_0, t_eval=t)


from odeint I get the following output:

    /home/marc/.cache/pypoetry/virtualenvs/data-science-6CF2GDM8-py3.9/lib/python3.9/site-packages/scipy/integrate/odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
    warnings.warn(warning_msg, ODEintWarning)
    
and solve_ivp returns:

          message: 'Required step size is less than spacing between numbers.'
         nfev: 626
         njev: 0
          nlu: 0
          sol: None
       status: -1
      success: False
            t: array([0.       , 0.0060006, 0.0120012, 0.0180018, 0.0240024, 0.030003 ,
           0.0360036, 0.0420042, 0.0480048, 0.0540054, 0.060006 , 0.0660066,
           0.0720072, 0.0780078, 0.0840084, 0.090009 , 0.0960096, 0.1020102,
           0.1080108])
     t_events: None
            y: array([[20.        , 20.00611389, 20.01247168, 20.0191072 , 20.02606478,
            20.03339925, 20.04117593, 20.04947141, 20.05838901, 20.06805067,
            20.07863631, 20.09039598, 20.10366405, 20.11892066, 20.1369306 ,
            20.15915028, 20.18809947, 20.23005987, 20.30820872],
           [ 1.        ,  1.03820542,  1.08182954,  1.13179425,  1.18950952,
             1.25687343,  1.33627217,  1.43064426,  1.54448509,  1.6825812 ,
             1.85388098,  2.07263638,  2.36066427,  2.75061513,  3.2992731 ,
             4.1588621 ,  5.64191939,  8.821076  , 20.38529206]])
     y_events: None

Can anybody point me in the right direction to solve this?

Thank you