Question: I would like to solve the following formulas numerically in Python;

[![enter image description here](https://i.stack.imgur.com/fiOSU.png)](https://i.stack.imgur.com/fiOSU.png)

In Mathematica, you can input multiple differential equations and solve it at the same time. Is there a way to do the similar thing with Scipy?

edit: yes, I have looked at scipy.integrate.odeint already, but I'm still not sure how I can solve multiple equations that correlates each other at the same time. Does anyone have suggestion for that?

Answer: 
Eventually I figured out myself and I'm writing down for other people who might be clueless like me;
in `odeint`, you give three parameters `model`, `y`, and `t` where `model` is a function that takes in `y` and `t` then return dydt, `y` is the initial value of `y`, and `t` is the variable you're trying to take integral over. If you have multiple differential equations that are dependent on each other, you can just pass in all of them to `odeint`. In my case,

```py
t = np.linspace(0, 20) # range of t
y0 = [No_0, Na_0, Ni_0, Nn_0] # initial condition for each Ns
def model(y, t):
    No, Na, Ni, Nn = y
    dNodt = -k_oa * No
    dNadt = k_oa * No - k_ai * Na
    dNidt = k_ai * Na - k_in * Ni
    dNndt = k_in * Ni
    return [dNodt, dNadt, dNidt, dNndt]
y = odeint(model, y0, t)

```

You can define multiple differential equations you want to solve within the model you define, and pass in to `odeint` along with the initial values for all species.
Have a loot at SciPy's Integration package:

<https://docs.scipy.org/doc/scipy/reference/integrate.html>

specifically at `odeint` to solve systems of ODE's:

<https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint>