Merry Christmas, I think I can still say it. I am back this time for my problem in a simpler case, without coupling. I have a problematic code now and I want your help. I want to solve numerically the below coupled equations:




the primes are $d/dx$ and μ the mass scale $μ=Μ/Μ_{planck}$

with boundary conditions for $\tilde{\phi}$ to be: $\lim_{x\rightarrow \infty}\tilde{\phi}(\tilde{x})=\tilde{\phi}_{fv}\sim 1,\qquad \frac{d\tilde{\phi}}{d\tilde{x}}\Bigg|_{\tilde{x}=0}=0$

and for $\tilde{\rho}:$

$\tilde{\rho}(0)= 0,\qquad \tilde{\rho}'(0)= 1$

with the following potential

$ \tilde{U}(\tilde{\phi})=\frac18 (\tilde{\phi}^2-1)^2+\frac{\tilde{\epsilon}}{2}(\tilde{\phi}-1)$

I'm following the method of this paper https://arxiv.org/pdf/2205.03140v2.pdf (if you have questions about my coding go to the APPENDIX A.6 and check the method's description in the flat and curved space cases, it is quick 3 pages ), where they use Taylor expansion to transform the BCP to a numerically solvable ICP. I want to solve for small $\epsilon$ lets say 0.11.

Here it is. My amateur, super problematic code. I post it as detailed as I can in order to be clear to you.

ro = r[x]; fi = f[x]; (*ρ(x) and φ(x)*)
rod = D[ro, x];
rodd = D[ro, {x, 2}];
fid = D[fi, x];
fidd = D[fi, {x, 2}];
eps = 0.11 (*ε*); xmin = 0.01; xmax = 12; fv = 0.93; 
f0 = -1.05;(*α in the paper, I put it equal to true vacuum value ~ -1, 
because I want the graph to start from there*)
eq1 = fidd + (3 rod/ro) fid - (1/2) fi (fi^2 - 1) - eps/2 == 0;
eq2 = rodd + (8 Pi/ 3) mrat^2 ro (fid^2 + (1/8) (fi^2 - 1)^2 + (eps/2) (fi - 1)) == 0;
(*initial and boundary conditions*)
incon1 = f[xmin] == f0 + (1/8) ((1/2) f0 (f0^2 - 1) + eps/2) xmin^2;
incon2 = f'[xmin] == (1/4) ((1/2) f0 (f0^2 - 1) + eps/2) xmin;
incon3 = r[xmin] == xmin;
incon4 = r'[xmin] == 1;
bcon1 = f[xmax] == fv;
bcon2 = f'[xmin] == 0;
bcon3 = r[0] == 0;
bcon4 = r'[0] == 1;
solution = NDSolveValue[{eq1, eq2, bcon1, bcon2, bcon3, bcon4}, {fi, ro}, {x,xmin,xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {incon1, incon2, incon3, incon4}}];
fsol = Part[Evaluate[f[x] /. solution], 1];
Plot[fsol, {x, xmin, xmax}, AxesLabel -> {x, φ}];
(* Manipulate Command for "f0" and "xmin-xmax"??*)
  • I think I need Manipulate command (for f0 and xmin-xmax?), but I can't use it in my Mathematica because it is cracked. When I used it shows an error.

  • For a solution I expect an instanton, which a good guess is $\phi(x)\sim \tanh(smthing)$. Check the graph ($\xi=0$ our case and $\eta$ is our $x$) enter image description here

  • When I run the code, stiff errors arise. I observed that in a number of shooting method questions here too, but I have no idea how to overcome these stiff problems.

  • I have already checked the answer in the flat case question here with a little different boundary conditions:How to program efficient undershoot/overshoot , I can solve the flat space problem too. Also, I can solve the curved space problem with another method, I have a code solving it via minimization of the energy functional, but I have to make a specific hypothesis for the function ρ(x) in the metric. So I want to solve it in a more general way via the equations of motion posted here.


Any ideas how to fix the boundary condition phi[xmax]=fv???

I'm leaving the code in its initial form in the case something its going wrong with the mass scale- eps change.

  • 1
    $\begingroup$ Some more information might help. What is the physical justification for thinking the solution should rise up to .93 for the potential you're using, and for choosing the given form for the potential? What is the significance of eps in your equation? Where does the plot in Fig. 3 above come from? It's not the same paper you linked to. I guess my point is if it's not in the math, maybe there's an error in the equation somewhere or in the physical assumptions. $\endgroup$
    – user87932
    Commented Jan 8, 2023 at 18:57
  • $\begingroup$ The output of the code can be verified by plugging the solutions back into the equations and boundary conditions after supplying a guess for 'a'. The equations should yield interpolating functions which are zero across the range of x, and the bc's should be satisfied. $\endgroup$
    – user87932
    Commented Jan 8, 2023 at 18:58
  • $\begingroup$ @jdp The plot and the potential are from here: sci-hub.se/https://www.worldscientific.com/doi/epdf/10.1142/…. The eps=0.11 as a choice is also from the paper of Lee and Lee (LL), ξ=0 is the case without coupling, the case we are studing in my question, the curved space from A.6 Appendix. I took the potential and the values LL used in their numerical section and tried to apply the insructions of Appendix A.6 from the other paper. Eps is used to be small in our approaches. Check the numerical section of the paper it is 1 page and it does not help me. $\endgroup$ Commented Jan 8, 2023 at 20:35
  • 1
    $\begingroup$ Thanks for the information. I'll try to set up a Manipulate with code to verify that the outputs do indeed properly solve the equations. I don't expect ParametricNDSolve to fail, since was released some years back, but it can't hurt to check. I'll add it to my answer when I get it ready - probably tomorrow. I'll look over the paper you linked. The thing is, there aren't many parameters you can vary; you can try adding eps to the Manipulate and see if that does anything. "Can it make great changes?" Well, the equations are nonlinear, so it's hard to say. $\endgroup$
    – user87932
    Commented Jan 9, 2023 at 2:18
  • 1
    $\begingroup$ Aside from this, the only other suggestions I have for now are to see if there's a way to adapt the code to reproduce another published result you know the answer to, and to try and treat gravity as a small perturbation - decrease the coupling strength to approach the flat solution. I don't know how hard the latter might be, but you did mention you knew what the solution should look like for this case. Maybe I'll come up with something better after looking over the paper. $\endgroup$
    – user87932
    Commented Jan 9, 2023 at 2:21

1 Answer 1


I deleted most of my earlier answer after reading the paper, and deleted or folded earlier comments into this answer.

The authors point out that the problem isn't a Cauchy initial value problem (IVP), but a boundary value problem (BVP), so they're using a shooting method. For flat space, they keep the phi'[0] boundary condition, but replace the phi[infinity] = fv condition (which makes the original problem a BVP instead of an IVP) with phi[0] = a, 'a' being a parameter to be varied. You now have a Cauchy IVP.

Due to the singularity at the origin, they need to shift their initial values over to xmin, which they do via a Taylor expansion. This gives you a new set of IVP's at xmin. You can now solve the resulting problem as a function of parameter 'a', then vary 'a' to match the condition at phi[infinity] = fv to desired precision.

For the curved space case, they follow a similar approach, but also shift the r variable over to xmin using a similar Taylor expansion, which they don't show the derivation for, just the result. This is how they get r[xmin] = xmin, etc. So their "shooting method" is partially implicit in the derivation they performed, except for the dependence on a parameter 'a'. I.e. they've restructured the BVP as a new IVP, and shifted the IVP equations over to xmin instead of 0 to avoid the singularity.

One thing that perhaps wasn't stressed enough was that this shift to xmin needs to be applied to both the boundary and initial value conditions. Three of the four are the same - only the condition on phi at infinity is swapped out to convert from a BVP to an IVP. The shooting method will attempt to test its intermediate results against the supplied boundary conditions, and if they remain singular, the process won't converge.

Once these steps are performed, you now have a Cauchy problem which can be solved using the equations and the initial values they derived at xmin, using ParametricNDSolve without any shooting method selected, and parameterized by 'a'. The result can be passed to FindRoot with an initial guess for 'a', and solved to match the boundary condition at "infinity" for phi, i.e. match to fv.

You could also try using the shooting method option for NDSolve instead, using the shifted initial and boundary conditions. The shooting tutorial mentions that FindRoot is used internally to seek convergence. There is an advantage to using ParametricNDSolve/FindRoot separately however; if you run into problems, it's easier to see what's going wrong. As long as ParametricNDSolve runs successfully, it's possible to use Plot/Table/etc. to obtain information about what your solution looks like even if FindRoot fails, as will be seen below.

After running an earlier version of this code, it was found that no solution was possible with the specified parameters. After further discussion/investigation with the OP, we found that the variable 'eps' was being used in two different ways. It controls the shape of the potential function the OP uses, but the authors of the referenced appendix in the link use it as the ratio of the scalar field mass to the Planck mass. I've changed the name of this parameter to massRatio in the code below.

It also became apparent that it would be useful to be able to try out different equations for the potential function. To simply this process, and avoid having to re-code equations by hand, I modified the earlier code to allow defining potential functions externally and passing them into the body of the solver. The necessary derivatives are computed internally and used to generate the appropriate equations for each potential.

This is the original potential function:

u["baseline"][eps_][phi_[x_]] := 
  1/8 (phi[x]^2 - 1)^2 + eps/2 (phi[x] - 1);

The [phi[x]] term at the end allow various derivatives to be performed. The [eps_] is for the constant parameter used in the equation; additional parameters may be added here if needed. The u["baseline"] is a way to label trial potentials. Any short descriptive string can serve as a label. It's also possible to omit the label and just use u0, u1, etc.

Here's the main function:

solve[phi_, r_, x_, xmin_, xmax_, a_, a0_, fv_, mrat_, u_] :=
  Module[{dua, eq},
   (* Compute the derivative of the potential wrt phi evaluated at \
phi=a *)
   dua = D[u, phi[x]] /. phi[x] -> a;
   eq["phi"] = phi''[x] + 3 r'[x]/r[x] phi'[x] - D[u, phi[x]] == 0; 
   eq["r"] = r''[x] + 8 Pi/3 mrat^2 r[x] (phi'[x]^2 + u) == 0;
   eq["ic"] = {phi[xmin] == a + 1/8 dua xmin^2, 
     phi'[xmin] == 1/4 dua xmin, r[xmin] == xmin, r'[xmin] == 1};
    Flatten[{eq["phi"], eq["r"], eq["ic"]}], {phi, r}, {x, xmin, 
     xmax}, {a}]

Call it with

 {xmin,xmax,a0,fv,mrat,eps} = {.01,12,-1.05,.93,.25,.11};
    pnds = solve[phi, r, x, xmin, xmax, a, a0,fv,mrat,

Then solve for the optimal 'a' value:

a = (a /. FindRoot[phi[a][xmax] == fv /. pnds, {a, a0}])

and plot the result (I scaled r down to better fit both curves on the plot and added the Epilog to the plot at 'fv' to serve as a reference marker.):

Plot[Evaluate[{phi[a][x], r[a][x]/10} /. pnds], {x, xmin, xmax}, 
         PlotLegends -> {"phi[x]", "r[x]/10"},Epilog ->  Line[{{0,.93},{xmax,.93}}]]

There was also a request for a way to run this code with Manipulate for various inputs. In order to get acceptable output whether FindRoot converges or not, I wrapped the call with Enclose/ConfirmQuiet, which will throw a Failure box to the surrounding enclose if messages are generated.

 Module[{pnds, phi, r, x, a}, 
  pnds = solve[phi, r, x, xmin, xmax, a, a0, fv, mrat, 
  {Plot[Evaluate[{phi[a0][x], r[a0][x]/10} /. pnds], {x, xmin, xmax}, 
    PlotLegends -> {"phi[x]", "r[x]/10"}, 
    Epilog -> Line[{{0, fv}, {xmax, fv}}]], 
     a = (a /. FindRoot[
         phi[a][xmax] == fv /. pnds, {a, a0}])]]}], {{xmin, .01}, .01,
   1}, {{xmax, 12}, 1, 31}, {{a0, -1.05}, -1.1, 0}, {{fv, 0.93}, 
  0.093, .93},
 {{eps, .11}, -.2, .2},
 {{mrat, .25}, .0001, 1.}]

Here's a sample output:

[![enter image description here][1]][1]

It's also possible to verify that the function does in fact give valid solutions. Write the equations/boundary/initial conditions as lhs == 0, choose a value for the parameter a, generate interpolating functions, and plug these into the lhs of each equation. The equations should give zero for all x in the range, and the initial/boundary conditions at their endpoints. Here's code to test the result:

 Module[{pnds, phi, r, x, a, eq1, eq2, ic, ifn, subst, sol},
  pnds = 
   solve[phi, r, x, xmin, xmax, a, a0, fv, mrat, 
  eq1 = phi''[
     x] + (3 r'[x]/r[x]) phi'[x] - (1/2) phi[x] (phi[x]^2 - 1) - eps/2;
  eq2 = r''[
     x] + (8 Pi/3) eps^2 r[
      x] (phi'[x]^2 + (1/8) (phi[x]^2 - 1)^2 + (eps/2) (phi[x] - 1));
  ic = {phi[xmin] - (a + (1/8) ((1/2) a (a^2 - 1) + eps/2) xmin^2), 
    phi'[xmin] - ((1/4) ((1/2) a (a^2 - 1) + eps/2) xmin), 
    r[xmin] - xmin, r'[xmin] - 1,
    {phi[xmax] - fv}};
  ifn = {phi[a0], r[a0]} /. pnds;
  subst = Thread[{phi, r} -> ifn];
  sol = {{eq1, eq2}, ic} /. a -> a0 /. subst;
  Column[{Plot[sol[[1, 1]], {x, xmin, xmax}], 
    Plot[sol[[1, 2]], {x, xmin, xmax}], sol[[2]]}]
 , {{xmin, .01}, .01, 1}, {{xmax, 12}, 1, 31}, {{a0, -1.05}, -1.1, 
  0}, {{fv, 0.93}, 0.093, .93},
 {{eps, .11}, -.2, .2},
 {{mrat, .25}, .0001, 1.}]

I'm seeing numerical round off errors on the order of 5x10^-6 overall, so the solutions appear valid.

  • $\begingroup$ Thanks a lot, I will run it and I will comment the results you posted in the next days, in general the algorithm i try to follow in my code is in the papers appendix ( without code there), ill check it in the next days because I am out of town, thanks again a lot $\endgroup$ Commented Dec 28, 2022 at 16:18
  • $\begingroup$ My φ and ρ depend only to x (η in the plot), which is Euclidean time here and the problem has an o(4) symmetry $\endgroup$ Commented Dec 28, 2022 at 16:23
  • $\begingroup$ I will check it in a few days and I will comment on your results. I also refer the manipulate command on my notes in the end, on xmin Xmax and maybe a, but I can't use it. I will check the 15day mathematica free trial, and I will back with more info. Thanks a lot again for your help, I will give you some credits in my acknowledgement section $\endgroup$ Commented Dec 31, 2022 at 10:14
  • $\begingroup$ @Jennifer Derleth I'll see if I can add a Manipulate (or Dynamic) later. It shouldn't be too difficult. (Famous last words.) $\endgroup$
    – user87932
    Commented Dec 31, 2022 at 16:11
  • 1
    $\begingroup$ @Jennifer Derleth I added the Manipulate and cleaned up a couple typos. In addition to the 15 day free trial, you might be able to get a free online account. You get charged cloud credits for using curated data and Manipulate, but as far as I know, straight numerical computation (i.e. NDSolve, Plot, etc.) isn't charged. You can check the WRI site for details. $\endgroup$
    – user87932
    Commented Dec 31, 2022 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.