trying to solve differential equations simultaneously

I am trying to build a code for chemical reactor design which is able to solve for the pressure drop, conversion, and temperature of a reactor. All these parameters have differential equations, so i tried to define them inside a function to be able to integrate them using ODEINT. However it seems that the function i've built has an error which i can't figure out which held me back from integration it.

the error that i'm encountering:

--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-32-63c706e84be5> in <module> 1 X0=[0.05,1200,2] ----> 2 y=func(X0,0) <ipython-input-27-6cfd4fef5ee2> in func(x, W) 8 kp=np.exp(((42311)/(R*T))-11.24) 9 deltah=-42471-1.563*(T-1260)+0.00136*(T**2 -1260**2)- 2,459*10e-7*(T**3-1260**3) ---> 10 ra=k*np.sqrt((1-X)/X)*((0.2-0.11*X)/(1-0.055*X)*(P/P0)-(x/(kp*(1-x)))**2) 11 summ = 57.23+0.014*T-1.94*10e-6*T**2 12 dcp=-1.5625+2.72*10e-3*T-7.38*10e-7*T**2 TypeError: unsupported operand type(s) for -: 'int' and 'list' 

and here is the full code:

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint fi = 0.45 gas_density = 0.054 #density Pres0 = 2 #pressure visc = 0.09 U = 10 Ac = 0.0422 T0 = 1400 #and 1200 also gc = 4.17 * 10**8 bed_density = 33.8 Ta = 1264.6 #wall temp fa0 = 0.188 def func(x,W): X = x[0] T = x[1] P = x[2] P0 = 2 R = 0.7302413 k = np.exp((-176008 / T) - (110.1 * np.log(T) + 912.8)) kp = np.exp(((42311) / (R * T)) - 11.24) deltah = -42471 - 1.563 * (T - 1260) + 0.00136 * (T**2 - 1260**2) - 2,459 * 10e-7 * (T**3 - 1260**3) ra = k * np.sqrt((1 - X) / X) * ((0.2 - 0.11 * X) / (1 - 0.055 * X) * (P / P0) - (x / (kp * (1 - x)))**2) summ = 57.23 + 0.014 * T - 1.94 * 10e-6 * T**2 dcp = -1.5625 + 2.72 * 10e-3 * T - 7.38 * 10e-7 * T**2 dxdw = 5.31 * k * np.sqrt((1 - X) / X) * ((0.2 - 0.11 * X) / (1 - 0.055 * X) * (P / P0) - (x / (kp * (1 - x)))**2) dpdw = (((-1.12 * 10**-8) * (1 - 0.55 * X) * T) / P) * (5500 * visc + 2288) dtdw = (5.11 * (Ta - T) + (-ra) * deltah) / (fa0 * (summ + x * dcp)) return [dxdw, dpdw, dtdw] X0 = [0.05, 1200, 2] y = func(X0, 0) 

thanks in advance