How to solve spatially discretised PDEs (method of lines) in solve_ivp or ODEint?

Question

I can discretise the spatial domain of a system of PDEs using the method of lines, converting the system of PDEs to a system of ODEs (with a time derivative only).

These equations (for context they are mass and energy balances) are here: https://imgur.com/gallery/MYdl5KI

It can be seen that each equation on the right hand side has a number of position dependent variables (where j denotes position, i denotes species which is irrelevant here). I have seen many examples of people solving time derivative ODEs in solve_ivp/ODEint. However, the difference here is that I need to solve these equations at a series of position values also.

From what I have learned to solve for multiple positions at a single time interval I will need a for loop to set up a vector for the ODEs at each position. This vector will then be passed to solve_ivp/odeint.

My question is - how can this for loop be set up in such a way that the ODEs can be solved at a number of positions?

Any help appreciated on this

Edit: My thoughts are the vector that will be accepted by solve_ivp/odeint will look something like this:

def my_system(t, x, a, b, c, d, etc): 
dxdt = np.zeros(4 * n) #setting up zeros array 
i=0

## Here the initial and boundary conditions are set (flow rate, temperature at inlet, initial temperature of the system, etc would be set) 
     Tb0 = ..
     M0 = ...
     Tg0 = ...
     mj0 = ...
     y0 = ...    

 ##The initial conditions are then used here to set up each ODE at the boundary position

##initial boundary values for the 4 ODEs
    dxdt[i] = ... 
    dxdt[i + 1] = ... 
    dxdt[i + 2] = ... 
    dxdt[i + 3] = ...

Then using the last mj equation to get the flowrate out of the boundary element
mj = ...

##Solution array y will look like [y0,M0,Tg,Tb,...yn,Mn,Tgn,Tbn] at all values of t

for i in range(2, 4 * n, 4): 
##starting at index 2, and incrementing by 4 because we are solving     for 4 values at each point
    dxdt[i] = ...
    dxdt[i + 1] = ... 
    dxdt[i + 2] = ... 
    dxdt[i + 3] = ...

return dxdt