Solve_ivp initial condition set up for 2D problem/issue

Question

So I am trying to solve a project from MIT Numerical Methods course, basically need to solve a now uncoupled equation for supersonic flow over an airfoil which is just a parabolic shape.but for some reason I can't make it work, feel very dumb. Here is the project: https://learning-modules.mit.edu/service/materials/groups/149130/files/1ccf396b-e6ee-4867-b648-2c0c5eb313eb/link?errorRedirect=%2Fmaterials%2Findex.html&download=true

So basically I have a 100x200 mesh and need to solve a PDE (now two ODE) at each node. In the mit lectures they actually combine a finite difference matrix for the spatial part (so they actually solve using finite difference) but the time is solved using matlab ODE45. I am using python since it is what it is more used now.

So I have a parabolic shape and have a different value of y at each point. So there should be 100 ODEs each representing one y. And have 2 ODEs so a shape of (2,100) but when I set up the initial condition if I set up like an array of shape (2,100) it needs to be one dimensional but if I flatten it, it says the shape doesnt match. What am I supposed to do so it matches?

import numpy as np
import itertools
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


## General Variables
N=100
r = np.linspace(0, 1, N)
s = np.linspace(-1, 1, N*2)
R, S = np.meshgrid(r, s)
y_max = 1
y_min = -1
y_upper = 0.1*r-0.1*r**2
y_lower = -0.1*r+0.1*r**2

Mach = 2 #given in problem

## Equation to solve def
def SSonicFlowFun(t, P):
    phix = P[0]
    phiy = P[1]
    dpxdt =(y_max/(y_max-y_upper))/(1-Mach**2) *phiy
    dpydt = (y_max/(y_max-y_upper))*phix
    return np.array([dpxdt, dpydt])

t_span = (0, 1)
Phi0 = np.zeros((2, 100))
print(Phi0.shape)


## Solving ODE
sol= solve_ivp(SSonicFlowFun, t_span, Phi0)

I eliminated the finite difference because it is throwing me off, still unsure how to incorporate it, so will try that next but since solve_ivp is solving the space too it should theoretically give same answer right? Sorry just trying to learn and everything here may or may not make sense so please correct me.

Your help is appreciated

Tried flattening the initial condition but it says I need to have a (2,100) shape but if I use the (2,100) shape it tells me initial condition needs to be 1-dimensional.

If I solve for a single ODE (shape 1,100) it works easily...