# Numerical solution for coupled ODE with variable coefficients

I'm trying to solve the system

Subscript[m, p] = 1.22*10^(19) * 1.52*10^(24);
Ag = Derivative[1][b][t]-1.5608310063810777811971842694*10^-43 b[t] Sqrt[(1 + 2.9079917009173806522983649735*10^-87 \[Psi][t]^2) (1.18253*10^149 f^4 (1 - (2.90799*10^-87 \[Psi][t]^2)/h^2)^2 + 1/2 Derivative[1][\[Psi]][t]^2)]

Bg = -((1.37552*10^63 f^4 \[Psi][t] (1 - (2.90799*10^-87 \[Psi][t]^2)/h^2))/h^2) - (3 \[Psi][t] Derivative[1][b][t]^2)/(4 b[t]^2 (1 + 2.9079917009173806522983649735*10^-87 \[Psi][t]^2)^2) + (3 Derivative[1][b][t] Derivative[1][\[Psi]][t])/ b[t] + Derivative[2][\[Psi]][t]


Using

sol = ParametricNDSolve[{Bg == 0, Ag == 0, \[Psi][0] == 2.1*(Subscript[m, p]),b[0] == 1, \[Psi]'[0] == 0}, {\[Psi], b}, {t, 0.0, 2.2*10^(-31)}, {f,h}, MaxSteps -> 100000, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,"Coefficients" -> FehlbergCoefficients,"StiffnessTest" -> False}];
Manipulate[ParametricPlot[{\[Psi][f, h][t], b[f, h][t]} /. sol, {t, 0,6*10^(-32)}, AspectRatio -> 1], {f, 0.1, 10}, {h, 0.1, 10}]


Where the FehlbergCoefficients is:

Fehlbergamat = {{1/4}, {3/32, 9/32}, {1932/2197, -7200/2197,
7296/2197}, {439/216, -8, 3680/513, -845/4104}, {-8/27,
2, -3544/2565, 1859/4104, -11/40}};
Fehlbergbvec = {25/216, 0, 1408/2565, 2197/4104, -1/5, 0};
Fehlbergcvec = {1/4, 3/8, 12/13, 1, 1/2};
Fehlbergevec = {-1/360, 0, 128/4275, 2197/75240, -1/50, -2/55};
FehlbergCoefficients[4, p_] :=
N[{Fehlbergamat, Fehlbergbvec, Fehlbergcvec, Fehlbergevec}, p];


But it fails with $f=5$ and $h=1$, (I got "\$abort"). Using$f = 5$and$h = 1$before integration with jut "NDsolve" the result is what I'm looking for. What I need is a plot$\psi \times t$and$b \times t$with variable coefficients h and f. • The quantity Subscript[m, p] is undefined. In general, it is wise not to use superscript or subscript variables in computations. – bbgodfrey Oct 11 '17 at 16:32 • Sorry: Subscript[m, p] = 1.22*10^(19) * 1.52*10^(24) – Kamog Oct 11 '17 at 16:53 • With v 11.2, I cannot reproduce the error you describe. InterpolatingFunction does complain that some numerical values of t lie outside its range, but that is not serious and can be fixed. – bbgodfrey Oct 11 '17 at 17:17 • @bbgodfrey, it is wierd, setting$f = 5$and$h=1$before integration the numerical values of$t$lies inside its range. My version is 10. – Kamog Oct 11 '17 at 17:28 ## 1 Answer The computation as described in the question fails, because FehlbergCoefficients is undefined. It also seems unnecessary. To produce the desired results, replace the second block of code by sol = ParametricNDSolve[{Bg == 0, Ag == 0, ψ[0] == 2.1*(Subscript[m, p]), b[0] == 1, ψ'[0] == 0}, {ψ, b}, {t, 0.0, 2.2*10^(-31)}, {f, h}]; Manipulate[ParametricPlot[{ψ[f, h][t], b[f, h][t]} /. sol, {t, 0, 6*10^(-32)}, AspectRatio -> 1], {{f, 5}, 0.1, 10}, {{h, 1}, 0.1, 10}]  If it is necessary for some reason to employ FehlbergCoefficients, use the code provided under "A Method of Fehlberg" here. Addendum Defining FehlbergCoefficients, as just added to the question, allows sol = ParametricNDSolve[{Bg == 0, Ag == 0, ψ[0] == 2.1*(Subscript[m, p]), b[0] == 1, ψ'[0] == 0}, {ψ, b}, {t, 0.0, 2.2*10^(-31)}, {f, h}, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, "Coefficients" -> FehlbergCoefficients, "StiffnessTest" -> False}];  to be used instead. It gives the same result. • My FehlbergCoefficients is defined... this plot is for$\psi$or$b\$ ? – Kamog Oct 12 '17 at 1:23
• @Kamog Where in your question is FehlbergCoefficients defined? The plot is of b as a function of ψ, as you specified by using ParametricPlot. If, instead, you want b and ψ as functions of t, use Plot instead. – bbgodfrey Oct 12 '17 at 1:29
• Sorry, I forgot to put it in my question, but I did it now.... thank you – Kamog Oct 12 '17 at 1:31
• Using: Manipulate[ Plot[{[Psi][f, h][t]} /. sol, {t, 0, 6*10^(-32)}, AspectRatio -> 1], {f, 0.1, 10}, {h, 0.1, 10}] I got the same error... =\ – Kamog Oct 12 '17 at 1:41
• I cannot reproduced the error. I suggest that you copy the definitions from your question plus the code in my answer and run them in a new notebook. – bbgodfrey Oct 12 '17 at 1:52