Solving a system of many coupled differential equations numerically

Question

I have the following system of coupled differential equations:

$$ r \frac{d}{dr}\left(r \frac{d \varphi_{m}}{dr} \right) - m^2 \varphi_{m} + r^2(F_0 \varphi_{m}+F_1 \varphi_{m-1}+F_2 \varphi_{m-2})=0, $$ where $m$ is an integer that ranges from $(m_0-k)$ to $(m_0+k)$. The $F$'s are known and I have some boundary conditions that depend on $m$.

The $\varphi$'s vanish for $m<-k$ or $m>k$. So the equation for $m=m_0-k$ is simply: $$ r \frac{d}{dr}\left(r \frac{d \varphi_{m_0-k}}{dr} \right) - (m_0-k)^2 \varphi_{m_0-k} + r^2F_0 \varphi_{m_0-k}=0. $$ So I need to know $\varphi_{m-2}$ and $\varphi_{m-1}$ to be able to solve the equation for $\varphi_m$ and iterate up to $(m_0+k)$.

I was thinking of having a list of InterpolatingFunction's. However, I don't know to handle the indices in a neat way.

Can someone please suggest an approach or point me some directions?

Answer 1

Maybe you can use the following as a skeleton: I handle the indices via creating the variables by string manipulation and ToExpression. Wether this is neat or not is another question.

m0 = -10;
k = 3;
F = {r \[Function] 1, r \[Function] 2, r \[Function] 3};
Y = Table[ToExpression["y" <> ToString[i]], {i, 1, 2 k + 1}];
Y0 = RandomReal[{-1, 1}, {2 k + 1}];
Y1 = RandomReal[{-1, 1}, {2 k + 1}];
t0 = 0.1;
tmax = 2.;
eq = Simplify@Table[
    Plus[
      1/r D[1/r D[Y[[m - (m0 - k) + 1]][r], r], r],
      -m^2 Y[[m - (m0 - k) + 1]][r],
      r^2 Plus[
        F[[1]][r] Y[[m - (m0 - k) + 1]][r],
        If[m - 1 >= m0 - k, F[[2]][r] Y[[m - 1 - (m0 - k) + 1]][r], 0],
        If[m - 2 >= m0 - k, F[[3]][r] Y[[m - 2 - (m0 - k) + 1]][r], 
         0]
        ]
      ] == 0,
    {m, m0 - k, m0 + k}];
initeqs = Flatten[Table[{Y[[i]][t0] == Y0[[i]], Y[[i]]'[t0] == Y1[[i]]}, {i, 1, 2 k + 1}]];

NDSolve[Join[eq, initeqs], Evaluate[Table[yy[r], {yy, Y}]], {r, t0, tmax}]

This returns replacement rules involving InterpolationFunction objects as you expected.