How to use the sympy.diff function in sympy

# incorrect of derivatives for higher orders, so there is no reasons to 
    # implement it. Plus, no uses have yet been found for the higher 
    # derivatives/integrals. For integrals, there is no way to get the 
    # indefinite integral.
    
    # Test 160K points in vector form for order 1, 2, and 3
    # Use lambdify for fast evaluation
    _Ts = np.linspace(5,500,20)
    _Tcs = np.linspace(501,900,20)
    _Pcs = np.linspace(2E5, 1E7,20)
    _omegas = np.linspace(-1, 10,20)
    _Ts, _Tcs, _Pcs, _omegas = np.meshgrid(_Ts, _Tcs, _Pcs, _omegas)
    _Ts, _Tcs, _Pcs, _omegas = _Ts.ravel(), _Tcs.ravel(), _Pcs.ravel(), _omegas.ravel()
    for order in range(1,4):
        B0c = diff(B0, T, order)
        B1c = diff(B1, T, order)
        Br = B0c + omega*B1c
        BVirial = (Br*R*Tc/Pc).subs(R, _R)
        f = lambdify((T, Tc, Pc, omega), BVirial, "numpy")
        
        Bcalcs = f(_Ts, _Tcs, _Pcs, _omegas)
        Bcalc2 = BVirial_Pitzer_Curl(_Ts, _Tcs, _Pcs, _omegas, order)
        assert_allclose(Bcalcs, Bcalc2)
        

    # Check integrals using SymPy:
    for order in range(-2, 0):
        B0c = B0
        B1c = B1
        for i in range(abs(order)):
            B0c = integrate(B0c, T)
            B1c = integrate(B1c, T)

pc["t"][tempn - ntimedata : tempn],
            prof,
            r2,
            1,
            ntimedata,
            range(pc["npc"]),
            range(spec),
            sharedata,
        )

        for p in range(spec):
            for i in range(pc["npc"]):
                pwt = 0
                for j in range(tempn - ntimedata, tempn):
                    profile = sharedata["pwalamo"]["profiles"][i][p][pwt]
                    targets[j, p] = diff(profile, st).subs(st, pc["t"][j])

                    if j == sharedata["pwalamo"]["tup"][i][p][pwt]:
                        pwt = pwt + 1

        tempn -= ntimedata
    return targets, sharedata

def Simple_manifold_with_scalar_function_derivative():
    coords = (x,y,z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3',metric='[1,1,1]',coords=coords)
    # Define surface
    mfvar = (u,v) = symbols('u v')
    X = u*e1+v*e2+(u**2+v**2)*e3
    print X
    MF = Manifold(X,mfvar)

    # Define field on the surface.
    g = (v+1)*log(u)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g,u) + (MF.rbasis[1]/MF.E_sq)*diff(g,v)
    print 'Vector derivative =', VectorDerivative.subs({u:1,v:0})

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print 'Vector derivative =', dg.subs({u:1,v:0})
    return

def compute_derivatives(self):
        symb = sp.symbols('s:%d' % len(self.ufl_operands))
        r = sp.diff(self.expr(*symb), *zip(symb, self.derivatives))
        return sp.lambdify(symb, r)

def ode_source_term(u):
    """Return the terms in the ODE that the source term
    must balance, here u'' + w**2*u.
    u is symbolic Python function of t."""
    return sm.diff(u(t), t, t) + w**2*u(t)

def test1():
    f = 2
    u = model1(f, L, D)
    print 'model1:', u, u.subs(x, 0), u.subs(x, L)
    print sp.latex(u, mode='plain')
    u = model2(f, L, C, D)
    print 'model2:', u, sp.diff(u, x).subs(x, 0), u.subs(x, L)
    print sp.latex(u, mode='plain')
    u = model3(0, 1+x**2, L, C, D)
    print 'model3:', u, u.subs(x, 0), u.subs(x, L)
    print sp.latex(u, mode='plain')

def ode_source_term(u, damping):
    """Return the terms in the ODE that the source term
    must balance, here m*u'' + f(u') + c*u.
    u is a symbolic Python function of t."""
    if damping == 'zero':
        return m*sym.diff(u(t), t, t) + c*u(t)
    elif damping == 'linear':
        return m*sym.diff(u(t), t, t) + \
               b*sym.diff(u(t), t) + c*u(t)
    else:  # damping is nonlinear
        return m*sym.diff(u(t), t, t) + \
               b*sym.diff(u(t), t)*abs(sym.diff(u(t), t)) + c*u(t)

xl = symbols('xl')
yl = symbols('yl')
z = symbols('z')
t = symbols('t')
pi = symbols('pi')

# Define the manufactured solution

n = 0.9 + 0.9*x + 0.2*sin(5.*x**2 - 2*z) + cos(10*z)*cos(y)*sin(y*7+1.234)



# Calculate gradients for boundaries

dndx = diff(n, x)
dndy = diff(n, y)

# Calculate RHS function

dndt = (
        Dx * D2DX2(n)
      + Dy * D2DY2(n)
      + Dz * D2DZ2(n)
      )

# Calculate sources

Sn = diff(n, t) - dndt

#x and y-domains are Lx and Ly long, respectively
# change to corresponding coordinates. Standard BOUT++ coordinates 
#have Lx = 1, Ly = 2pi. 

def newton_raphson(func: str, a: int, precision: int = 10 ** -10) -> float:
    """ Finds root from the point 'a' onwards by Newton-Raphson method
    >>> newton_raphson("sin(x)", 2)
    3.1415926536808043
    >>> newton_raphson("x**2 - 5*x +2", 0.4)
    0.4384471871911695
    >>> newton_raphson("x**2 - 5", 0.1)
    2.23606797749979
    >>> newton_raphson("log(x)- 1", 2)
    2.718281828458938
    """
    x = a
    while True:
        x = Decimal(x) - (Decimal(eval(func)) / Decimal(eval(str(diff(func)))))
        # This number dictates the accuracy of the answer
        if abs(eval(func)) < precision:
            return float(x)

all_derivatives = []
    for e in expressions:
        indexeds = retrieve_indexed(e.rhs, mode='all', deep=True)
        adjoint_output_fn, adjoint_mapper = diff_indexed(e.lhs, adjoint_mapper)
        derivatives = []
        for i in indexeds:
            i_d, adjoint_mapper = diff_indexed(i, adjoint_mapper)
            state = extract_le_state(e)

            old_ds = state['dspace']

            new_parts = {i_d.function: old_ds.parts[i.function]}
            new_ds = DataSpace(old_ds.intervals, new_parts)
            state['dspace'] = new_ds

            d_eqn = LoweredEq(i_d, i_d+diff(e.rhs, i)*adjoint_output_fn, **state)
            derivatives.append(d_eqn)

        for i, d in enumerate(derivatives):
            subs = {}
            for e_ind, d_ind in zip(e.lhs.indices, d.lhs.indices):
                ind = retrieve_dimension(d_ind).pop()
                if ind.is_Time:
                    continue
                else:
                    if (d_ind - e_ind) != 0:
                        subs[ind] = e_ind - d_ind
            derivatives[i] = shift_le_index(d, subs)

        all_derivatives += derivatives
    return all_derivatives