# Solving ODE with RK method and plotting

I am trying to solve the following ode

$$i\dot{\phi_n}(t)=-\phi_{n+1}(t)-\phi_{n-1}(t)+\frac{g~ \phi_n(t)}{1+|\phi_n(t)|^2}$$ with I.C. $$\phi(0)=-2\cos(2.5)$$ and $$g=1$$.

However, I am unable to produce the "close to desired" plot of $$|\phi_n(t=fixed)|$$ y-axis vs $$n$$ on x-axis. What I wrote in the code is the following

k = 2.5;
n = 1000;
syms2 = Table[Subscript[\[Alpha], i], {i, n}];
For[j = 1, j < n + 1, j++, syms2[[j]] = 0];
syms2[[(n/2)]] = 1; syms2[[n/2 + 1]] = 1;
Subscript[\[Phi], 0][t_] := Subscript[\[Phi], 1][t];
Subscript[\[Phi], n + 1][t_] := Subscript[\[Phi], n][t];
eqns = Table[{Sqrt[-1]*
Subscript[\[Phi], i]'[t] == -Subscript[\[Phi], i + 1][t] -
Subscript[\[Phi], i - 1][t] + ((
syms2[[i]]*Subscript[\[Phi], i][t])/(
1 + Abs[Subscript[\[Phi], i][t]]^2)),
Subscript[\[Phi], i][0] == -2*Cos[k]}, {i, n}];
vars = Table[Subscript[\[Phi], i][t], {i, n}];
sol = NDSolve[eqns, vars, {t, 0, 250},
Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False},
MaxSteps -> \[Infinity], AccuracyGoal -> 8, PrecisionGoal -> 8];
Plot1 = ListPlot[Table[Abs[vars /. First[%]], {t, 0, 250}],
PlotRange -> All, Mesh -> False]


I don't know how to get a snapshot at a particular time, as I mentioned above.

• The equations in your code is inconsistent with the ones shown in $\LaTeX$, please double check them. Nov 19 '18 at 7:35
• For example, -Subscript[\[Phi], i + 1][t] - Subscript[\[Phi], i - 1][t] vs $-\phi_{n+1}(t)-\phi_{n+1}(t)$. Nov 19 '18 at 8:17
• i.stack.imgur.com/eA3aR.png Nov 19 '18 at 10:50
• I can't, because nobody except you knows which version is correct, and it's not impossible that both versions are incorrect. The only thing I can do is to point out the inconsistency. Nov 19 '18 at 11:12
• i.stack.imgur.com/4CNhm.png Nov 19 '18 at 11:36

tm = 1000; k = 2.5;
n = 1000;
Table[syms2[i] = 0, {i, n}];
syms2[n/2] = 1; syms2[n/2 + 1] = 1;
\[Phi][n + 1][t_] := \[Phi][n][t]
\[Phi][0][t_] := \[Phi][1][t]
eqns = Table[{I*
D[\[Phi][i][t], t] == -\[Phi][i + 1][t] - \[Phi][i - 1][t] +
syms2[i]*\[Phi][i][t]/(1 + Abs[\[Phi][i][t]]^2), \[Phi][i][
0] == -2*Cos[k]}, {i, n}];
vars = Table[Abs[\[Phi][i][t]], {i, n}];
sol = NDSolveValue[eqns, vars, {t, 0, tm}];

L = ListAnimate[
Table[ListPlot[sol, PlotRange -> {0, 2.5}, ColorFunction -> Hue,
Joined -> True, Frame -> True], {t, 0, tm, 5}]]


• How did you create the GIF movie of this in order to post it?
– Moo
Nov 19 '18 at 18:28
• @Moo I use Export["C:\\Users\\name\\Desktop\\Fig1.gif", L, AnimationRepetitions -> Infinity] Nov 19 '18 at 18:40
• @AtoZ Yes, x-axis is {i,1,n} and y-axis is Table[Abs[\[Phi][i][t]], {i, n}] at n=1000 Nov 20 '18 at 12:04
• @AtoZ Do you want a solution with broken symmetry, for example, n=500; 'Table[syms2[i] = 0, {i,-n, n}]; syms2[n/2] = 1; syms2[n/2 + 1] = 1;`? Nov 20 '18 at 12:16
• @AtoZ To reproduce your pictures you need to know the parameters and the initial conditions. Where did you get them from? Nov 21 '18 at 12:38