I try to solve a system($n$) of coupled nonlinear 1-d PDE's which is basically coupled sine-Gordon equations (through 2nd spatial derivative) in time and space. I am interested in the dynamics when I increase the force gradually. So I use the final answer of smaller force $\beta$ and use it as initial condition for a larger force $\beta$. I wrote a code for smaller force and it produces the result in a decent amount of time (depending on the $n$). I use the final result as initial condition in the second part of my code. Since both parts of my code contain basically the same type of equations and only different in initial conditions and force magnitude, I expect it to run the same way as the first part. But, unfortunately no matter what I tried, the second part takes very long time and runs forever. Here is my code:
n = 3; \[Eta] = 0.1; \[Beta] = 0.4; L = 10; T = 100; \[Tau] = T/10;
af = Array[f, n];
txaf[x_, t_] := Through[af[x, t]];
eqs = Join[
Array[D[f[#][x, t], t, t] + \[Eta] D[f[#][x, t], t] -
D[Sum[i f[i][x, t], {i, 1, #}] + #*Sum[f[j][x, t], {j, # + 1, n}], x, x] +
Sin[f[#][x, t]] - \[Beta] Tanh[t/\[Tau]] == 0 &,n]];
init = {Union[Thread[Drop[txaf[x, 0], {(n + 1)/2}] ==
0], {txaf[x, 0][[(n + 1)/2]] == 4 ArcTan[Exp[x]]}],
Thread[Derivative[0, 1][txaf][x, 0] == 0]};
bc = {Thread[Derivative[1, 0][txaf][-L, t] == 0],
Thread[Derivative[1, 0][txaf][L, t] == 0]};
sol = NDSolve[{eqs, init, bc}, af, {x, -L, L}, {t, 0, T},
MaxStepSize -> 0.05]
using as intial condition:
f0[x_] := Thread[First[Evaluate[txaf[x, T] /. sol]]];
df0[x_] := Thread[First[Evaluate[Derivative[0, 1][txaf][x, T] /. sol]]];
second part of my code:
Clear[af, eqs, init, bc]
\[Beta]2 = 0.5;
af = Array[f, n];
txaf[x_, t_] =.;
txaf[x_, t_] := Through[af[x, t]];
eqs = Join[
Array[D[f[#][x, t], t, t] + \[Eta] D[f[#][x, t], t] -
D[Sum[i f[i][x, t], {i, 1, #}] + #*Sum[f[j][x, t], {j, # + 1, n}], x, x] +
Sin[f[#][x, t]] - \[Beta] - (\[Beta]2 - \[Beta])*
Tanh[t/\[Tau]] == 0 &, n]];
init = {Thread[txaf[x, 0] == f0[x]],
Thread[Derivative[0, 1][txaf][x, 0] == df0[x]]};
bc = {Thread[Derivative[1, 0][txaf][-L, t] == 0],
Thread[Derivative[1, 0][txaf][L, t] == 0]};
sol2 = NDSolve[{eqs, init, bc}, af, {x, -L, L}, {t, 0, T},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 1000, "MinPoints" -> 500, "DifferenceOrder" -> "Pseudospectral"}}]
you can copy and paste the code and just hit the run. In the mean time if you have any comment to make the code more efficient for higher number of $n$, it would be greatly appreciated.