# Plotting basins of attraction of magnetic pendulum in spherical coordinates

I am new in Mathematica, and this is the first time I coding. Following is the code solving differential equations of motions in spherical coordinates of magnetic pendulum and plotting basins of attraction. I got idea from this code:

X[1] = 1;
X[2] = -(1/2);
X[3] = -(1/2);
Y[1] = 0;
Y[2] = Sqrt[3]/2;
Y[3] = -Sqrt[3]/2;
nf = Nearest[{{1, 0}, {-0.5, Sqrt[3]/2}, {-0.5, -Sqrt[3]/2}} ->
Automatic] /* First;
getStateData[k_, c_, h_, x0_, y0_] :=
First@NDSolveProcessEquations[{x''[t] + k x'[t] + c x[t] -
Sum[(X[i] -
x[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i,
3}] == 0,
y''[t] + k y'[t] + c y[t] -
Sum[(Y[i] -
y[t])/(h^2 + (X[i] - x[t])^2 + (Y[i] - y[t])^2)^(3/2), {i,
3}] == 0, x'[0] == 0, y'[0] == 0, x[0] == x0, y[0] == y0}, {x,

sd = getStateData[.15, .2, .2, 1, 1];

getBasin2[x0_, y0_] :=
Module[{state = sd, sol},
state = First@NDSolveReinitialize[state, {x[0] == x0, y[0] == y0}];
NDSolveIterate[state, 100.5];
sol = {x[100], y[100]} /. NDSolveProcessSolutions[state];
nf[sol]]
ArrayPlot[
ParallelTable[
getBasin2[xpos, ypos], {xpos, -2, 2, 0.1}, {ypos, -2, 2, 0.1}],
ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming


And following is my code.(It's quite long)

mags = {{1.2, 0, 0, 1, 0, 0, -1}, {1.2, Pi/6, 0, 1, -Sin[Pi/6],
0, -Cos[Pi/6]}, {1.2, Pi/6, 2*Pi/3,
1, -Sin[Pi/6]*Cos[2*Pi/3], -Sin[Pi/6]*
Sin[2*Pi/3], -Cos[Pi/6]}, {1.2, Pi/6, 4*Pi/3,
1, -Sin[Pi/6]*Cos[4*Pi/3], -Sin[Pi/6]*Sin[4*Pi/3], -Cos[Pi/6]}};

nf = Nearest[{{-0.082788, 0}, {0.041394,
0.071696}, {0.041394, -0.071696}} -> Automatic] /* First;

getStateData[k_, c_, M_, m_, L_, g_, tht0_, phi0_] :=
First@NDSolveProcessEquations[{tht''[t] ==
Sin[tht[t]]*Cos[tht[t]]*
phi'[t]^2 - (Sum[
c*mags[[i,
4]]*((mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
mags[[i, 7]]*
Sin[tht[t]]) + (3*
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Cos[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
Sin[phi[t]] -
Cos[mags[[i, 2]]]*
Sin[tht[
t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]])) -
3*(mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
mags[[i, 7]]*Sin[tht[t]])*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))/(((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^1) -
15*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Cos[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
Sin[phi[t]] -
Cos[mags[[i, 2]]]*
Sin[tht[t]])*((mags[[i, 5]]*Sin[tht[t]]*
Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^2)/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^1.5, {i,
4}] + (m + M/2)*g*L*Sin[tht[t]] +
k*(L^2)*tht'[t]/3)/((M/3 + m)*L^2),
phi''[t] == -Sum[
c*mags[[i,
4]]*((-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*
Cos[phi[t]]) + (3*
mags[[i,
1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Sin[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Cos[phi[
t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]])) -
3*(-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Cos[phi[t]])*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^1) -
15*mags[[i,
1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Sin[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*

Cos[phi[t]])*((mags[[i, 5]]*Sin[tht[t]]*
Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] + Cos[mags[[i, 2]]]*Cos[tht[t]]))/(1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^2)/(1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^1.5, {i,
4}]/(((M/3 + m)*L^2)*Sin[tht[t]]^2) -
2*tht'[t]*phi'[t]/Tan[tht[t]] -
k*(L^2)*phi'[t]/(3*((M/3 + m)*L^2)), tht'[0] == 0,
phi'[0] == 0, tht[0] == tht0, phi[0] == phi0}, {tht, phi}, t,

sd = getStateData[0.1, 0.0004, 0.05, 0.01, 0.3, 9.8, 1, 3];

getBaSin[tht0_, phi0_] :=
Module[{state = sd, sol},
state = First@
NDSolveReinitialize[state, {tht[0] == tht0, phi[0] == phi0}];
NDSolveIterate[state, 110.];
sol = {tht[100], phi[100]} /. NDSolveProcessSolutions[state];
nf[{0.3*Sin[sol[[1]]]*Cos[sol[[2]]],
0.3*Sin[sol[[1]]]*Sin[sol[[2]]]}]]

ArrayPlot[
ParallelTable[
getBaSin[tht, phi], {tht, 0.1, 1.5, 0.1}, {phi, 0.1, 6.2, 0.1}],
ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming


First I want to know what the red lines mean, which says such like:

NDSolveReinitialize::dvnoarg : -- Message text not found -- (0.1),

since I expected the result as picture colored evenly with three colors.

And second I want to know how the result will differ if I change the NDSolve method(ex)Adams -> ExplicitRungeKutta), and which method will be most appropriate.

• May 8, 2019 at 12:08

I fixed the code for my taste. There's only a typo when calling getBaSin[tht, phi]. You can not use the parameters tht, phi. We can reduce the time to t = 20. Fig.1 is different from what is obtained on the plane. In my opinion it is even more beautiful.

mags = {{1.2, 0, 0, 1, 0, 0, -1}, {1.2, Pi/6, 0, 1, -Sin[Pi/6],
0, -Cos[Pi/6]}, {1.2, Pi/6, 2*Pi/3,
1, -Sin[Pi/6]*Cos[2*Pi/3], -Sin[Pi/6]*
Sin[2*Pi/3], -Cos[Pi/6]}, {1.2, Pi/6, 4*Pi/3,
1, -Sin[Pi/6]*Cos[4*Pi/3], -Sin[Pi/6]*Sin[4*Pi/3], -Cos[Pi/6]}};

nf = Nearest[{{-0.082788, 0}, {0.041394,
0.071696}, {0.041394, -0.071696}} -> Automatic] /* First;
eq = {tht''[
t] - (Sin[tht[t]]*Cos[tht[t]]*
phi'[t]^2 - (Sum[
c*mags[[i,
4]]*((mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
mags[[i, 7]]*
Sin[tht[t]]) + (3*
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Cos[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
Sin[phi[t]] -
Cos[mags[[i, 2]]]*
Sin[tht[
t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]])) -
3*(mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
mags[[i, 7]]*Sin[tht[t]])*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))/(((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^1) -
15*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Cos[tht[t]]*
Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
Sin[phi[t]] -
Cos[mags[[i, 2]]]*
Sin[tht[t]])*((mags[[i, 5]]*Sin[tht[t]]*
Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^2)/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))^1.5, {i,
4}] + (m + M/2)*g*L*Sin[tht[t]] +
k*(L^2)*tht'[t]/3)/((M/3 + m)*L^2)) == 0,
phi''[t] - (-Sum[
c*mags[[i,
4]]*((-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*
Cos[phi[t]]) + (3*
mags[[i,
1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Sin[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*

Cos[phi[
t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]])) -
3*(-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Cos[phi[t]])*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]])))/((1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^1) -
15*mags[[i,
1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Cos[phi[t]])*((mags[[i, 5]]*Sin[tht[t]]*
Cos[phi[t]] +
mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
mags[[i, 7]]*Cos[tht[t]]) -
mags[[i,
1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
Cos[mags[[i, 3]]] +
mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] + Cos[mags[[i, 2]]]*Cos[tht[t]]))/(1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
Sin[tht[t]]*Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^2)/(1 +
mags[[i, 1]]^2 -
2*mags[[i,
1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
Cos[phi[t]] +
Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
Sin[phi[t]] +
Cos[mags[[i, 2]]]*Cos[tht[t]]))^1.5, {i,
4}]/(((M/3 + m)*L^2)*Sin[tht[t]]^2) -
2*tht'[t]*phi'[t]/Tan[tht[t]] -
k*(L^2)*phi'[t]/(3*((M/3 + m)*L^2))) == 0, tht'[0] == 0,
phi'[0] == 0, tht[0] == tht0, phi[0] == phi0};

getStateData[op_] :=
First@NDSolveProcessEquations[eq /. op, {tht, phi}, t,
op = {k -> .1, c -> .0004, M -> .05, m -> .01, L -> .3, g -> 9.8,
tht0 -> 1, phi0 -> 3};
sd = getStateData[op]

getBaSin[x_, y_] :=
Module[{state = sd, sol},
state = First[NDSolveReinitialize[sd, {tht[0] == x, phi[0] == y}]];
NDSolveIterate[state, 21];
sol = {tht[20], phi[20]} /. NDSolveProcessSolutions[state];
nf[{0.3*Sin[sol[[1]]]*Cos[sol[[2]]],
0.3*Sin[sol[[1]]]*Sin[sol[[2]]]}]]

ArrayPlot[
ParallelTable[
getBaSin[x, y], {x, 0.1, 1.5, 0.01}, {y, 0.1, 6.2, 0.025}],
ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue}] // AbsoluteTiming


• Thank you so much! May I ask you a few more questions? First, what is the principle of this code solving differential equations numerically? I mean I thought there should be somewhere defining the time step, which is not shown in my code. Second, does NDSolveIterate[state, 21] means that function works about 20 times? Is it related to the time step? Lastly, is Adams' method usually better than RK4 method?
– 최준호
May 10, 2019 at 16:44
• Sorry, I accidentally edited my question... I hope you are able to see my edited comment
– 최준호
May 10, 2019 at 17:19
• 1) Read the tutorial section Advanced Numerical Differential Equation Solving in the Wolfram Language -> Components and Data Structures. 2) 21 is a time of integration (t,0,21}, it is related to the time step rough as 21/number of steps. 3) Adams method is 20% faster than Automatic whereas with the Method -> "ExplicitRungeKutta"` you can not find a solution. May 10, 2019 at 17:26