I have 15 ode equations need to be solved simultaneously and I want to solve them using solve\_ivp.

There are each 5 states for T, co2, and q. The initial conditions are T=20, co2 = 0, q=0

I tried to separate them into 3 lists, one for T, one for co2, and one for q.

I am not sure how to resolve this bug and have worked on it for couple hours.
Really appreciate your help!

```
import math
from scipy.integrate import solve_ivp
from bokeh.core.property.instance import Instance
from bokeh.io import save
from bokeh.layouts import column
from bokeh.model import Model
from bokeh.models import CustomJS, Slider, Callback
from bokeh.plotting import ColumnDataSource, figure, show
import numpy as np

# three plots , co2 as y and z as x

###############    User generated - Slider initial value   ############### 
V= 100.0 # volume
r = 5.0
T = 20.0
c_co2_0 = 5.0 # concentration
episl_r = 0.3 # void
v0 = 2.0 # initial vilocity

###############  ---          Static Parameters        ---  ############### 
b0 = 93.0 * (10**(-5))
deltH_0 = 95.3 # calculate b
Tw = -5.0 # room temperature
T0 = 353.15 # temeperature 
t0 = .37 # heterogeneity constant, in paper is denoted as t_h0
alpha = 0.33
chi = 0.0 
q_s0 = 3.40
R = 8.314
kT = 3.5*(10**3) #calculate rA
ρs = 880.0
deltH_co2 = 75.0 # calculate temeprature change 

# ------------------ For Equation 4 : Enegergy Ballance  --------------
ρg = 1.87 # ?
h = 13.8
Cp_g = 37.55 # J/molK
Cp_s = 1580.0 # J/molK

###############   -----  Parameters depend on input  -----  ############### 
L = V / (math.pi * r**2)
deltZ = L / 5.0 # 5 boxes in total
p_co2 = R * T * c_co2_0
a_s = deltZ / r
theta = (1-episl_r) * ρs * Cp_s +  episl_r * ρg * Cp_g

# Equations are calclulated in order 
def b(T):
    b = ( b0 ** ( (deltH_0/ (R * T0) ) * (T0/T - 1) ) )
    return b

def t_h(T):
    return ( t0 + alpha * (1 - T0 / T) )

def q_s(T):
    return ( q_s0 ** ( chi * (1 - T / T0)) )

# Calculate rco2_n (not ode)
# change it to q
def R_co2(T, c_co2, q): 
    b_var = b(T)
    t_var = t_h(T)
    qs_var = q_s(T)
    # print(qs_var)
    r_co2 =  kT * ( R * T * c_co2 * ( (1- ( (q / qs_var)**t_var) )**(1/t_var) ) - q / (b_var*qs_var) ) 
    # print(r_co2)
    return r_co2
# ODE Part 
# Repetitive shortcut

# Equation 2
ener_balan_part1 = v0  * ρg* Cp_g 
def ener_balan(theta, deltZ): # replace v0  * ρg* Cp_g / (theta * deltZ) 
    return(ener_balan_part1/ (theta*deltZ) )
def ener_balan2(episl_r):
    return( (1-episl_r) * ρs * deltH_co2)
def ener_balan3(a_s, Tw, T0):
    return (a_s * h *(Tw-T0))

# Equation 1 Mass Balance : find co2_n

def mass_balan(episl_r, deltZ):
    return ( v0/ (episl_r * deltZ) )
def masss_balan2(episl_r, ρs):
    return( (1-episl_r ) * ρs )
def deriv(t, y):
    T_n, co2_n, q_n = y
    # rco2_ first, rate of generation 
    T1 = -ener_balan(theta, deltZ) * T_n + ener_balan(theta, deltZ) * T0 + ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n))+ ener_balan3(a_s, Tw, T0)
    co2_1 = -mass_balan(episl_r, deltZ) * co2_n + mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
    q_1 = R_co2(T_n, co2_n, q_n)

    T2 = -ener_balan(theta, deltZ) * T_n + ener_balan(theta, deltZ) * T0 + ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n))+ ener_balan3(a_s, Tw, T0)
    co2_2 = -mass_balan(episl_r, deltZ) * co2_n + mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
    q_2 = R_co2(T_n, co2_n, q_n)

    T3 = -ener_balan(theta, deltZ) * T_n + ener_balan(theta, deltZ) * T0 + ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n))+ ener_balan3(a_s, Tw, T0)
    co2_3 = -mass_balan(episl_r, deltZ) * co2_n + mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
    q_3 = R_co2(T_n, co2_n, q_n)

    T4 = -ener_balan(theta, deltZ) * T_n + ener_balan(theta, deltZ) * T0 + ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n))+ ener_balan3(a_s, Tw, T0)
    co2_4 = -mass_balan(episl_r, deltZ) * co2_n + mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
    q_4 = R_co2(T_n, co2_n, q_n)

    T5 = -ener_balan(theta, deltZ) * T_n + ener_balan(theta, deltZ) * T0 + ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n))+ ener_balan3(a_s, Tw, T0)
    co2_5 = -mass_balan(episl_r, deltZ) * co2_n + mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
    q_5 = R_co2(T_n, co2_n, q_n)

    T_ls = np.array([T1, T2, T3, T4, T5])
    co2_ls = np.array([co2_1, co2_2, co2_3, co2_4, co2_5])
    q_ls = np.array([q_1, q_2, q_3, q_4, q_5])

    return T_ls, co2_ls, q_ls

t0, tf = 0, 10
############# initial condition 
T_initial = 20
c_co2_0 = 0
q0 = 0
init_cond = np.array([20, 0, 0])
N=5
soln = solve_ivp(deriv, (t0, tf), init_cond)

```

Here is the error message

```
helper.py:293: RuntimeWarning: divide by zero encountered in double_scalars
  r_co2 =  kT * ( R * T * c_co2 * ( (1- ( (q / qs_var)**t_var) )**(1/t_var) ) - q / (b_var*qs_var) )
Traceback (most recent call last):
  File "helper.py", line 350, in <module>
    soln = solve_ivp(deriv, (t0, tf), init_cond)
  File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/ivp.py", line 546, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
  File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/rk.py", line 96, in __init__
    self.h_abs = select_initial_step(
  File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py", line 104, in select_initial_step
    d1 = norm(f0 / scale)
ValueError: operands could not be broadcast together with shapes (3,5) (3,) 

```