Is 2D ode integration possible in Stan?

I am attempting to fit several SIR models at once to model a Bayesian Hierarchical Model.

The current signature for the ODE integrators is
real[] ODEFunc(real t, real[] y, real[] params, real[] x_r, int[] x_i)

Is there any way to solve an ODE which looks something like this?
real[,] ODEFunc(real t, real[,] y, real[,] params, real[] x_r, int[] x_i)

Essentially, if we have a differential equation that looks like this:
dydt[1] = params[1] * y[2]
in the 1d case, dydt[1], params[1], and y[2] are all reals. Can the ODE solver work in the case where dydt[1], params[1], and y[2] are real[n], where n is the number of different fits I want?

Sorry if this was asked somewhere else or if I’m ignorant of something. I’m not really well versed in how the ODE solvers work , and I couldn’t find a related topic in the forum that answered this question.

So could you fit n totally independent odes? Independent in the sense that their states do not depend on one another.

If your odes are not independent, then you would have to flatten the structure yourself. In case the odes are independent, then solve them one by one.