# Can NDSolve be initialized with a close solution?

I have a set of 4 first order ODE with 4 initial conditions to be met. It models a simple damped pendulum like the one described here. What changes is the damping term and I switched to first order through Pontryagin's principle. In the model, there is a particular parameter that is making things difficult to solve. It is the parameter called "k" in the code below. Stiffness seems to increase with its value.

Ideally, I would like for k to reach up to the value of 10. However, in my trials, whenever k > 0.7, NDSolve struggles too much and is not able to further solve the set of ODE successfully.

Therefore, I would like to know if it is possible to initialize NDSolve with the solution at k -> 0.7, which is successful, as a starting point to search for the solution when k -> 0.8 (which sould be not to far from the previous solution). This way I could be able to form a loop in order to gradually increase the value of k.

Here is the code sample:

params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2.0;
solF = NDSolve[{q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(
Izz + a^2 m),
p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m),
q[0.] == 0.,
z[0.] == 0.,
q[T] == 0.8,
z[T] == 0.} /. params /. k -> 0.7,
{q, z, p, x},
{t, 0, T},
Method -> {"TimeIntegration" -> "ExplicitRungeKutta"}
] // AbsoluteTiming


Can this be accomplished?

EDIT 1

I have modified the above code according to the first suggestion in the comments. However, the proposed solution does not work for this particular example because it is only 4 additional initial conditions. I first solved the system with the previous code above to obtain a successful solution for k -> 0.6. I then assign that solution to oldSolF in order to start the new algorithm.

oldSolF = solF; (*obtained with the previous code with k->0.6*)
vars = {q, z, p, x};
varsF[t_] = Through[vars[t]];
params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2.0;
solF = NDSolve[{q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(
Izz + a^2 m),
p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m),
q[0.] == 0.,
z[0.] == 0.,
q[T] == 0.8,
z[T] == 0.} /. params /. k -> 0.7,
{q, z, p, x},
{t, 0, T},
Method -> {"Shooting",
"StartingInitialConditions" ->
Thread[varsF[0] == (varsF[0] /. oldSolF[[2, 1]])]}
] // AbsoluteTiming


However, at the end I get two error messages indicating that NDSolve failed to solve the system:

FindRoot: Failed to converge to the requested accuracy or precision within 100 iterations.

NDSolve: The scaled boundary value residual error of 4.2168606782853276*^7 \
indicates that the boundary values are not satisfied to specified \
tolerances. Returning the best solution found.


Any suggestions?

COMMENT

I believe that @bbgodfrey contribution answers the initial question most straighforwardly: it uses standard Mathematica functions (makes use of the InitialSeeding option from NDSolve). However, it should be highlighted that @AlexTrounev and @xzczd solutions are more robust and allow for a higher value of the k parameter.

• Something like vars = {q, z, p, x}; varsF[t_] = Through[vars[t]]; and Method -> {"Shooting", "StartingInitialConditions" -> Thread[varsF[0] == (varsF[0] /. oldsolF[[1]])]} after saving a successful solF as oldsolF = solF. Aug 31, 2022 at 15:59
• This problem can be solved with the local discontinuous Galerkin method (LDG) in the range $0.7\le k\le 3.6$. At k=3.7 we have same message from FindRoot as above. Substituting data back in NDSolve we can see that NDSolve reproduces LDG solution up to k=1.5, and over k=1.5 the solution deviates from LDG. Sep 1, 2022 at 15:25
• @AlexTrounev Is the LDG method an option in NDSolve? I cannot find it in the documentation... Would you share your code? I think that already going up to k=3.6 is a very good performance increase. Sep 1, 2022 at 17:40
• @Meclassic Actually LDG method is not an option in NDSolve. It is not even mentioned in Mathematica while it widely use as an effective solver in combination with FEM for ODEs and PDEs. Sep 1, 2022 at 19:04
• I am interested in what these equations describe, what mechanical phenomenon. Provide a link to this phenomenon? Sep 2, 2022 at 9:16

This interesting problem can be solved, at least up to k = 4.6, using standard Mathematica functions. To begin, turn off two extraneous warning messages for convenience.

Off[NDSolveValue::femibcnd]
Off[NDSolveFEMInitializePDECoefficients::femcscd]


Then use the NDSolve finite element option with InitialSeeding to obtain a solution for k = 1.

params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2.0;
{q0[t], z0[t], p0[t], x0[t]} =
NDSolveValue[{q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(Izz +
a^2 m), p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m),
DirichletCondition[{q[t] == 0, z[t] == 0}, t == 0],
DirichletCondition[{q[t] == 0.8, z[t] == 0}, t == T]} /.
params /. k -> 1, {q[t], z[t], p[t], x[t]}, {t} \[Element]
Line[{{0}, {T}}],
InitialSeeding -> {q[t] == 0.8 t/T, z[t] == 0, p[t] == 0, x[t] == 0},
Method -> {"FiniteElement", "MeshOptions" -> MaxCellMeasure -> 0.001}];
Join[%[[1 ;; 2]] /. t -> 0, %[[1 ;; 2]] /. t -> T]

(* {0., 0., 0.8, 0.} *)

GraphicsGrid[
{{Plot[%%[[1]], {t, 0, T}, AxesLabel -> q,
LabelStyle -> {12, Bold, Black}, PlotRange -> All],
Plot[%%[[2]], {t, 0, T}, AxesLabel -> z,
LabelStyle -> {12, Bold, Black}, PlotRange -> All]},
{Plot[%%[[3]], {t, 0, T}, AxesLabel -> p,
LabelStyle -> {12, Bold, Black}, PlotRange -> All],
Plot[%%[[4]], {t, 0, T}, AxesLabel -> x,
LabelStyle -> {12, Bold, Black}, PlotRange -> All]}}]


For larger values of k, use InitialSeeding based on the k = 1 solution just obtained. Specifically,

InitialSeeding -> {q[t] == q0[t] 3/4 + .2 t/T, z[t] == z0[t] 3/4,
p[t] == p0[t] 3/4, x[t] == x0[t]/2}


Now, we might expect to use {q[t] == q0[t], z[t] == z0[t], p[t] == p0[t], x[t] == x0[t]}, but for some unknown reason it does not work at all well. With the choice above, k up to 3.8 can be solved with the result,

Then, using the solution for k = 3.8, we can obtain results up to k = 4.6.

InitialSeeding -> {q[t] == q00[t] 3/4 + .2 t/T, z[t] == z00[t] 3/4,
p[t] == p00[t] 3/4, x[t] == x00[t]/3}


I have been unsuccessful so far in obtaining solutions for larger k.

• I did not know about the InitialSeeding option and it is what I was looking for from the beginning. However, could you comment as to why you scale down the previous solutions by 3/4 for example and particularly, why use q[t] == q0[t] 3/4 + .2 t/T. Are these arbitrary choices in order to ensure NDSolve success? Or are these factors specifically chosen for some reason? Could you elaborate on this? Sep 2, 2022 at 14:50
• @Meclassic These arbitrary choices in order to ensure NDSolve success. I plan to look at this situation further over the weekend. InitialSeeding is a fairly new option. Sep 2, 2022 at 22:39
• Another idea: in the hyperplastic solid mechanics tutorial is a section on solving extremely non liner problems. A second alternative shown there is the method of using a pseudo time integration. Sep 3, 2022 at 4:34
• @bbgodfrey Very nice solution (+1). Probably we need Do loop to update InitialSeeding same as I do it with FindRoot. Sep 3, 2022 at 4:58

Let consider the local discontinuous Galerkin method in application to solve system of ODEs. The theory is discussed here. The implementation is very straightforward. In this example we use Euler polynomials and Gauss formula for integration

params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2.0; kmax = 30;

Get["NumericalDifferentialEquationAnalysis"];
UT[mm_, t_] := EulerE[mm, t];

M0 = 5; nn = 12; ns = 4; np = 5; tmax = 2;
y[t_] = Table[Symbol["y" <> ToString[i]][t], {i, 1, ns}];
eqs = {q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(Izz + a^2 m),
p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m)}; F =
Table[eqs[[i, 2]], {i, ns}]; var = {q[t], z[t], p[t],
x[t]}; ini = {0, 0, pb, xb}; bc2 = {8/10, 0}; ff =
Table[F /. Table[var[[i]] -> y[t][[i]], {i, ns}] /.
params /. {k -> (7 + kap)/10}, {kap, 0, kmax}];
LDGODEs[M0_, nn_, ns_, np_, f_, ini_, tmax_, kap_] :=
Module[{dx = tmax/nn,   A = Array[aa, {M0 + 1, nn, ns}]},
xl = Table[ l*dx, {l, 0, nn}];
psi[mm_, k_, t_] :=
Piecewise[{{UT[
mm, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
gg = Table[
GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
dp = Table[
D[UT[mm, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
t], {k, nn}];
rul = Table[
y[t][[i]] ->
Sum[aa[mm, k, i] psi[mm - 1, k, t], {mm, 1, M0 + 1}, {k, 1,
nn}], {i, ns}];
eq = Flatten[
Table[-Sum[
aa[i + 1, n,
ks] (Table[(psi[i, n, t] If[j == 0, 0,
dp[[n]] /. mm -> j]), {t, gg[[n]][[All, 1]]}] .
gg[[n]][[All, 2]]), {i, 0,
M0}] - (Table[(f[[ks]] /. rul) psi[j, n, t], {t,
gg[[n]][[All, 1]]}] . gg[[n]][[All, 2]]), {j, 0, M0}, {n,
1, nn}, {ks, ns}] +
Table[Sum[
aa[i + 1, n,
ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0,
M0}] - Sum[
If[n < 2, ini[[ks]]/(M0 + 1),
aa[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n,
xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks, ns}]];
bc = Table[
Sum[aa[i + 1, k, s] psi[i, k, 2], {i, 0, M0}, {k, nn}] ==
bc2[[s]], {s, 1, 2}];
eqn = Table[eq[[k]] == 0, {k, Length[eq]}];
var = Join[Flatten[A], {pb, xb}];
soln = FindRoot[Join[eqn, bc],
Table[{var[[i]], If[kap == 0, 1/10, var[[i]] /. soln]}, {i,
Length[var]}]];
soln];


Solution

Do[
ldgsol[kap] =
LDGODEs[M0, nn, ns, np, ff[[kap + 1]], ini, tmax, kap];, {kap, 0,
kmax}] // AbsoluteTiming


In the end we have one message from FindRoot

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.


It means that step with k=3.7 is a limit of this code with current parameters. Visualization

leg = {"q", "z", "p", "x"};

Table[Plot[
Evaluate[
Table[Sum[aa[i + 1, k, s] psi[i, k, t], {i, 0, M0}, {k, nn}] /.
ldgsol[kap], {kap, 0, kmax, 1}]], {t, 0, tmax}, PlotRange -> All,
PlotLegends -> Table[(7. + n)/10, {n, 0, kmax, 1}], Frame -> True,
PlotLabel -> leg[[s]]], {s, ns}]



Using this solution we can generate initial data for NDSolve and generate new solution

data = Table[{.7 + kk/10, pb, xb} /. ldgsol[kk], {kk, 0, kmax, 1}];
sol1 = Table[
NDSolveValue[{q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(Izz +
a^2 m), p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m),
q[0] == 0, z[0] == 0, p[0] == data[[i, 2]],
x[0] == data[[i, 3]]} /. params /. k -> data[[i, 1]], {q, z, p,
x}, {t, 0, T}], {i, Length[data]}];


Finally we compare LDG numerical solution BVP and NDSolve solution with initial data only

Table[{.7 + i/10,
Table[Show[
Plot[sol1[[i + 1]][[s]][t], {t, 0, tmax}, PlotRange -> All],
Plot[Evaluate[
Sum[aa[i + 1, k, s] psi[i, k, t], {i, 0, M0}, {k, nn}] /.
ldgsol[i]], {t, 0, tmax}, PlotRange -> All,
PlotStyle -> {Red, Dashed}]], {s, ns}]}, {i, 0, kmax}]


Please note, that up to k=1.5 we have a good coincidence sol1 with ldgsol, but then solutions diverge, for example for k<1.5 and for k>1.5

Update 1. With option Method -> {"Newton", "StepControl" -> "TrustRegion"} for FindRoot we can extend solution to k=10.7, but there are oscillation around some points.

Code for this case

params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2.0; kmax = 50;

Get["NumericalDifferentialEquationAnalysis"];
UT[mm_, t_] := EulerE[mm, t];

M0 = 5; nn = 12; ns = 4; np = 5; tmax = 2;
y[t_] = Table[Symbol["y" <> ToString[i]][t], {i, 1, ns}];
eqs = {q'[t] == z[t],
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(Izz + a^2 m),
p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m)}; F =
Table[eqs[[i, 2]], {i, ns}]; var = {q[t], z[t], p[t],
x[t]}; ini = {0, 0, pb, xb}; bc2 = {8/10, 0}; ff =
Table[F /. Table[var[[i]] -> y[t][[i]], {i, ns}] /.
params /. {k -> (7 + 2 kap)/10}, {kap, 0, kmax}];
LDGODEs[M0_, nn_, ns_, np_, f_, ini_, tmax_, kap_] :=
Module[{dx = tmax/nn,   A = Array[aa, {M0 + 1, nn, ns}]},
xl = Table[ l*dx, {l, 0, nn}];
psi[mm_, k_, t_] :=
Piecewise[{{UT[
mm, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
xl[[k]] <= t <= xl[[k + 1]]}, {0, True}}];
gg = Table[
GaussianQuadratureWeights[np, xl[[i]], xl[[i + 1]]], {i, nn}];
dp = Table[
D[UT[mm, (2 t - xl[[k + 1]] - xl[[k]])/(xl[[k + 1]] - xl[[k]])],
t], {k, nn}];
rul = Table[
y[t][[i]] ->
Sum[aa[mm, k, i] psi[mm - 1, k, t], {mm, 1, M0 + 1}, {k, 1,
nn}], {i, ns}];
eq = Flatten[
Table[-Sum[
aa[i + 1, n,
ks] (Table[(psi[i, n, t] If[j == 0, 0,
dp[[n]] /. mm -> j]), {t, gg[[n]][[All, 1]]}] .
gg[[n]][[All, 2]]), {i, 0,
M0}] - (Table[(f[[ks]] /. rul) psi[j, n, t], {t,
gg[[n]][[All, 1]]}] . gg[[n]][[All, 2]]), {j, 0, M0}, {n,
1, nn}, {ks, ns}] +
Table[Sum[
aa[i + 1, n,
ks] (psi[i, n, xl[[n + 1]]] psi[j, n, xl[[n + 1]]]), {i, 0,
M0}] - Sum[
If[n < 2, ini[[ks]]/(M0 + 1),
aa[i + 1, n - 1, ks] psi[i, n - 1, xl[[n]]]] psi[j, n,
xl[[n]]], {i, 0, M0}], {j, 0, M0}, {n, 1, nn}, {ks, ns}]];
bc = Table[
Sum[aa[i + 1, k, s] psi[i, k, 2], {i, 0, M0}, {k, nn}] ==
bc2[[s]], {s, 1, 2}];
eqn = Table[eq[[k]] == 0, {k, Length[eq]}];
var = Join[Flatten[A], {pb, xb}];
soln = FindRoot[Join[eqn, bc],
Table[{var[[i]], If[kap == 0, 1/10, var[[i]] /. soln]}, {i,
Length[var]}],
Method -> {"Newton", "StepControl" -> "TrustRegion"}];
soln];

Do[
ldgsol[kap] =
LDGODEs[M0, nn, ns, np, ff[[kap + 1]], ini, tmax, kap];, {kap, 0,
kmax}] // AbsoluteTiming


Visualization

leg = {"q", "z", "p", "x"};

Table[Plot[
Evaluate[
Table[Sum[aa[i + 1, k, s] psi[i, k, t], {i, 0, M0}, {k, nn}] /.
ldgsol[kap], {kap, 0, kmax, 1}]], {t, 0, tmax}, PlotRange -> All,
PlotLegends -> Table[(7. + 2 n)/10, {n, 0, kmax, 1}],
Frame -> True, PlotLabel -> leg[[s]]], {s, ns}]

• Perhaps try the Affine Covariant Newton Method Sep 3, 2022 at 4:23
• @user21 Yes, you are right, with option Method -> {"Newton", "StepControl" -> "TrustRegion"} we can compute solution up to k=10, but x has oscillation in some points at k>4.5. Probably it is why solution with FEM proposed by bbgodfrey has a bound $k\le 4.6$. Sep 3, 2022 at 4:52
• @AlexTrounev I see that you updated your answer by adding the solutions given by the Affine Covariant Newton method. Could you post the code for the record? Sep 4, 2022 at 19:22
• @Meclassic Yes, code has been added. Sep 4, 2022 at 21:14

It's well-known the default "Shooting" method often doesn't work well for nonlinear BVP of ODE, and I'm tired of looking for proper initial guess of initial conditions for shooting method, so let me solve it with finite difference method (FDM).

First of all, it's easy to notice that z[t], x[t], p[t] can be easily eliminated, and further analysis suggests the original first order system seems to be troublesome. (To be more specific, it results in singular Jacobian in FindRoot). So let's simplify the ODE system to one ODE:

params = {a -> 0.4, m -> 2.4, Izz -> 1.0, g -> 9.81, fc -> 13.6829};
T = 2; (* Necessary change *)
Clear[z, x, p]
{z[t_], x[t_], p[t_]} =
Solve[##][[1, 1, -1]] & @@@
z'[t] == (-a g m Sin[q[t]] - fc Tanh[k z[t]] + x[t])/(Izz + a^2 m),
x'[t] == -p[t] + (fc k Sech[k z[t]]^2 x[t])/(Izz + a^2 m)},
{z[t], x[t], p[t]}}]

newsys = {p'[t] == (a g m Cos[q[t]] x[t])/(Izz + a^2 m),
q[0] == 0., z[0] == 0., q[T] == 0.8, z[T] == 0.} // Simplify


Next step is to generate the difference equation system, I'll use pdetoae for the task:

points = 200; difforder = 4; grid = Array[# &, points, {0, T}];
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[q[t], grid, difforder];

aesys = {#[[3 ;; -3]] &@ptoafunc[newsys[[1]]],
ptoafunc@newsys[[2 ;;]]} /. params // Flatten;


Then let's solve. For nonlinear FDM, we still need an initial guess, but proper initial guess for FDM is generally easier to find than that of shooting method (according to my limited experience). In your case, for k == 0.5, a naive guess is enough:

guess[t_] = 0;
sollst0 = FindRoot[aesys /. k -> 0.5 // Rationalize[#, 0] &,
Table[{q[t], guess[t]}, {t, grid}],
MaxIterations -> 500, WorkingPrecision -> 16][[All, -1]];


WorkingPrecision option can actually be taken away to speed up the code. You'll see FindRoot spit out warning FindRoot::lstol then, but this leads to no observable change in the solution.

Then we use this solution as the new guess:

Dynamic[kvalue]

sollstlst =
FoldList[FindRoot[aesys /. k -> (kvalue = #2) // Rationalize[#, 0] &,
Transpose[{q /@ grid, #}],
WorkingPrecision -> 16][[All, -1]] &,
sollst0,
Range[0.6, 10, 0.1]]; // AbsoluteTiming
(* {96.1869, Null} *)


Visualization:

ListLinePlot[sollstlst[[;; ;; 5]], DataRange -> {0, T},
PlotLegends -> Range[0.5, 10, 0.1][[;; ;; 5]]]


BTW, the last 50 solutions are virtually the same, so perhaps you don't need that large k:

ListLinePlot[sollstlst[[-50 ;; -1]], DataRange -> {0, T}]


The left work is to obtain z[t], x[t] and p[t]. This is harder than I thought: I have to use a very dense grid to surpress the non-physical oscillation of derivative (perhaps techniques like WENO will help in this case, but currently there's no built-in support for them):

qsol =
ListInterpolation[sollstlst[[-1]], grid, InterpolationOrder -> difforder];

pointsfinal = 10^4; gridfinal = Array[# &, pointsfinal, {0, T}];
ptoafuncfinal = pdetoae[q[t], gridfinal, difforder];

aesysfinal = {#[[3 ;; -3]] &@ptoafuncfinal[newsys[[1]]],
ptoafuncfinal@newsys[[2 ;;]]} /. params // Flatten;

sollstfinal =
FindRoot[aesysfinal /. k -> 10,
Table[{q[t], qsol[t]}, {t, gridfinal}]][[All, -1]]; // AbsoluteTiming
(* {20.7576, Null} *)

qsolfinal =
ListInterpolation[sollstfinal, gridfinal, InterpolationOrder -> difforder];

MapThread[Plot[#[t] /. q -> qsolfinal /. k -> 10 /. params //
Evaluate, {t, 0, T}, PlotRange -> All, PlotLabel -> #2] &, {#,
ToString /@ #}] &@{q, z, p, x} // GraphicsRow


In the final step I choose to live with the FindRoot::lstol warning, because otherwise the calculation is too slow.

• I find your solution outstanding. However, I tried to follow your code but I am having issues after the lines guess[t_] = 0; sollst0 = FindRoot[aesys /. k -> 0.5.... I have 2 errors: FindRoot: "The function value {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<190>>} is not a \ list of numbers with dimensions {200} at"... and FindRoot: Search specification 500 should be a list with 1 to 5 elements.. I imported all the necessary code related to your pdetoae function from the link. Any pointers? Sep 2, 2022 at 20:27
• @Meclassic Oops, I forgot I've modified T to T=2. Without this, the variable names won't match. (After Rationalize[……, 0]&, terms like q[0.], q[0.010050251256281407] become q[0], q[2/199]. ) Sep 2, 2022 at 22:41
• @xzczd It is nice solution (+1). I always trust in pdetoae, but try to test new method LDG. Sep 3, 2022 at 5:02
• @xzczd I am amazed by the robustness on your code. I can reach up to k->18.3 as it is, and when I use 400 points, it reaches up to k->39.1. Did you publish this method in a journal or something of the sort? Sep 15, 2022 at 16:06
• @Meclassic Nope, I haven't, you see, pdetoae is essentially an automation of the built-in NDSolveFiniteDifferenceDerivative :) . If you want to cite the code, you may directly cite this web page. (There's a related discussion in meta: mathematica.meta.stackexchange.com/q/968/1871 ) Sep 16, 2022 at 1:30