CFD simulation (with multiple for loops and matrix operations) is very slow to run. Looking to replace with faster numpy functions (or alternative)

Question

As mentioned above, the function below works, however its very slow. I am very interested in using faster/optimised numpy (or other) vectorized alternatives. I have not posted the entire script here due to it being too large.

My specific question is - are there suitable numpy (or other) functions that I can use to 1) reduce run time and 2) reduce code volume of this function, specifically the for loop?

Edit: mass, temp, U and dpdh are functions that carry out simple algebraic calculations and return constants

def my_system(t, y, n, hIn, min, mAlumina, cpAlumina, sa, V):
    dydt = np.zeros(3 * n) #setting up zeros array for solution (solving for [H0,Ts0,m0,H1,Ts1,m1,H2,Ts2,m2,..Hn,Tsn,mn])
# y = [h_0, Ts_0, m_0, ... h_n, Ts_n, m_n]
# y[0] = hin
# y[1] = Ts0
# y[2] = minL

i=0

## Using thermo
T = temp(y[i],P) #initial T
m = mass(y[i],P) #initial m

#initial values
dydt[i] = (min * (hIn - y[i]) + (U(hIn,P,min) * sa * (y[i + 1] - T))) / m # dH/dt (eq. 2)
dydt[i + 1] = -(U(hIn,P,min) * sa * (y[i + 1] - T)) / (mAlumina * cpAlumina) # dTs/dt from eq.3
dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt (holdup variation) eq. 4b
dydt[i + 2] = min - dmdt # mass flow out (eq.4a)

for i in range(3, 3 * n, 3): #starting at index 3, and incrementing by 3 because we are solving for 'triplets' [h,Ts,m] in each loop

    ## Using thermo
    T = temp(y[i],P)
    m = mass(y[i],P)

    # [h, TS, mdot]
    dydt[i] = (dydt[i-1] * (y[i - 3] - y[i]) + (U(y[i-3], P, dydt[i-1]) * sa * (y[i + 1] - T))) /  m # dH/dt (eq.2), dydt[i-1] is the mass of the previous tank
    dydt[i + 1] = -(U(y[i-3], P, dydt[i-1]) * sa * (y[i + 1] - T)) / (mAlumina * cpAlumina) # dTs/dt eq. (3)
    dmdt = dydt[i] * dpdh(y[i], P) * V # Equation 4b
    dydt[i + 2] = dydt[i-1] - dmdt # Equation 4a

return dydt

The functions mass, temp, U, and dpdh used inside the my_system function all take numbers as input, perform some simple algebraic operation and return a number (no need to optimise these I am just providing them for further context)

def temp(H,P):
    # returns temperature given enthalpy (after processing function)
    T = flasher.flash(H=H, P=P, zs=zs, retry=True).T
    return T

def mass(H, P):
    # returns mass holdup in mol
    m = flasher.flash(H=H, P=P, zs=zs, retry=True).rho()*V
    return m

def dpdh(H, P):
    res = flasher.flash(H=H, P=P, zs=zs, retry=True)
    if res.phase_count == 1:
        if res.phase == 'L':
            drho_dTf = res.liquid0.drho_dT()
        else:
            drho_dTf = res.gas.drho_dT()
    else:
        drho_dTf = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='P')
    dpdh = drho_dTf/res.dH_dT_P()
    return dpdh

def U(H,P,m):
    # Given T, P, m
    air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
    mu = air.mu*1000/mWAir #mol/m.s
    cp = air.Cpm #J/mol.K
    kg = air.k #W/m.K
    g0 = m/areaBed #mol/m2.s
    a = sa*n/vTotal #m^2/m^3 #QUESTIONABLE
    psi = 1
    beta = 10
    pr = (mu*cp)/kg
    re = (6*g0)/(a*mu*psi)
 hfs = ((2.19*(re**1/3)) + (0.78*(re**0.619)))*(pr**1/3)*(kg)/diameterParticle

    h = 1/((1/hfs) + ((diameterParticle/beta)/kAlumina))

    return h

Reference Image: enter image description here