I wish to solve the following set of ODE. $ $ i\frac{d}{dt}B_{n}\left(t\right) =f\sqrt{\left(P-n\right)\left(n+1\right)}B_{n+1}\left(t\right)+f\sqrt{n\left(P-n+1\right)}B_{n-1}\left(t\right) + Y\left[\left(P-n\right)\left(P-n-1\right)+n\left(n-1\right)\right]B_{n}\left(t\right)$ $ The question is similar to the Solving n simultaneous differential equation.It’s working perfectly for N=2 case. But the method is not giving any results for larger N with larger time say t=0 to 3000, and varying constants.

Values of constants are $ f=20$ , $ Y=334$ . Also $ n=0,1,2,…P$

When we consider a simple case, $ P=2$ we get 3 simulatneous ODE as follows. $ $ \frac{d}{dt}B_{0}\left(t\right) =-iYB_{0}\left(t\right)-i\sqrt{2}fB_{1}\left(t\right)$ $

$ $ \frac{d}{dt}B_{1}\left(t\right) =-i\sqrt{2}fB_{0}\left(t\right)-i\sqrt{2}fB_{2}\left(t\right)$ $

$ $ \frac{d}{dt}B_{2}\left(t\right) =-i\sqrt{2}fB_{1}\left(t\right)-iYB_{2}\left(t\right)$ $

The idea is to solve the above equation by finding eigen values and eigen vectors. We can convert the above equation in a matrix form as follows: $ $ \begin{pmatrix}\frac{d}{dt}B_{0}\left(t\right)\ \frac{d}{dt}B_{1}\left(t\right)\ \frac{d}{dt}B_{2}\left(t\right) \end{pmatrix}=\left(\begin{array}{ccc} -iY & -i\sqrt{2}f & 0\ -i\sqrt{2}f & 0 & -i\sqrt{2}f\ 0 & -i\sqrt{2}f & -iY \end{array}\right)\begin{pmatrix}B_{0}\left(t\right)\ B_{1}\left(t\right)\ B_{2}\left(t\right) \end{pmatrix}$ $

Let $ $ S=\left(\begin{array}{ccc} -iY & -i\sqrt{2}f & 0\ -i\sqrt{2}f & 0 & -i\sqrt{2}f\ 0 & -i\sqrt{2}f & -iY \end{array}\right)$ $

For solving the set of differential equations we need to find the eigen values and eigenvectors of $ S$ . When we calcualte them, it turns out to be 3 complex and distinct eigen values and 3 eigenvectors corresponding to each eigenvalues. Thus the solution will be:

$ $ B_{n}\left(t\right)=\sum_{n}G_{n}\overline{j}e^{\lambda t}$ $

provided: $ \overline{j}$ are eigen vectors and $ \lambda$ is the corresponding eigenvalue.

$ G_{n}$ can be solved using the initial condition, $ $ B_{n}\left(0\right)=\begin{pmatrix}1\ 0\ 0 \end{pmatrix}$ $

The above is just $ P=2$ case. Any way to extend the above for $ n$ case?

After that I need to plot the modulus square of the coeficcients with respect to time in a single plot.

What I have tried using NDSolve is:

 NN = 5;  f = 20    Y = 334     N1 = Sqrt[(NN - n)*(n + 1)];  N2 = Sqrt[n*(NN - n + 1)];  N3 = ((NN - n)*(NN - n - 1)) + (n*(n - 1));  FT = f*N1    ST = f*N2    TT = (Y/2)*N3    odes = Table[    I ToExpression["M" <> ToString[n]]'[t] ==      FT*ToExpression["M" <> ToString[n + 1]][t] +       ST*ToExpression["M" <> ToString[n - 1]][t] +       TT*ToExpression["M" <> ToString[n]][t], {n, 0, NN}];  deps = Table[ToExpression["M" <> ToString[n]][t], {n, 0, NN}]   {M0[t], M1[t], M2[t], M3[t], M4[t], M5[t]}   ic = {M0[0] == 1, M1[0] == 0, M2[0] == 0, M3[0] == 0, M4[0] == 0, M5[0] == 0}   MH = NDSolve[{odes, ic}, deps, {t, 0, 3000}]       Plot[{Evaluate[(M0[t]*Conjugate[M0[t]]) /. MH],    Evaluate[(M1[t]*Conjugate[M1[t]]) /. MH],    Evaluate[(M2[t]*Conjugate[M2[t]]) /. MH],    Evaluate[(M3[t]*Conjugate[M3[t]]) /. MH],    Evaluate[(M4[t]*Conjugate[M4[t]]) /. MH],    Evaluate[(M5[t]*Conjugate[M5[t]]) /. MH]}, {t, 0, 3000}, PlotRange -> All]