Infinite expression 1/0.^(1/4) encountered in NDSolve

I am trying to solve very simple system of coupled differential equations using the NDSolve function, but I am getting "Infinite expression 1/0.^(1/4) encountered."

Any idea how to solve the problem? I was looking for solutions where people suggested to shift the "Initial boundary point" or by using the Method -> {"Shooting", "StartingInitialConditions" ->...... Could not manage to resolve it.

The code is as follows:

alphaphi = 2*10^-13;
MX = 1.15*10^-12;
Mphi = 10^13;
alphaX = 10^-2;
gstar = 200;
HubbleI = Mphi;
Mpl = 1.22*10^19;

phiI = (3/(8*Pi))*Mpl^2/Mphi^2 *HubbleI^2/Mphi^2;

T[x] = (30/(gstar*Pi^2))^(1/4)*Mphi*R[x]^(1/4)/x;
Xeq[x] = MX^3/Mphi^3 *(2*Pi)^(-3/2)*x^3*(T[x]/MX)^(3/2)*
Exp[-MX/T[x]];
c1 = Sqrt[3/(8*Pi)]*Mpl*alphaphi/Mphi;
c2 = c1*Mphi*alphaX/(MX*alphaphi);
c3 = c2*Mphi/MX;

sol = NDSolve[{D[phi[x], x] == -c1*x*phi[x]/Sqrt[phi[x]*x + R[x]],
D[R[x], x] ==
c1*x^2*phi[x]/Sqrt[phi[x]*x + R[x]] +
c2*x^-1*(X[x]^2 - Xeq[x]^2)/Sqrt[phi[x]*x + R[x]],
D[X[x], x] == -c3*x^-2*(X[x]^2 - Xeq[x]^2)/Sqrt[phi[x]*x + R[x]],
phi[1] == phiI, R[1] == 0, X[1] == 0}, {phi, R, X}, {x, 1, 10^5}];

Plot[Evaluate[{phi[x], R[x], X[x]} /. sol], {x, 1, 10^5},
PlotRange -> All]

-

I suppose you know that you're dividing by R[x]^(1/4), since it's your equations. And therefore I assume the question is how to integrate the system. Since the division occurs inside an exponential and the argument of the exponential approaches negative infinity at R[1] == 0, the term stays bounded and approaches zero. One solution is to rewrite the Xeq function so that it is zero when R is zero. We also have to rewrite the T function.

Clear[Xeq, T];
T[x_, r_] := (30/(gstar*Pi^2))^(1/4)*Mphi*r^(1/4)/x;
Xeq[_, r_ /; r == 0] = 0.;
Xeq[x_?NumericQ, r_] :=
MX^3/Mphi^3*(2*Pi)^(-3/2)*x^3*(T[x, r]/MX)^(3/2)*Exp[-MX/T[x, r]];


To get past x == 1, we need to set the option Method -> "StiffnessSwitching".

alphaphi = 2*10^-13;
MX = 1.15*10^-12;
Mphi = 10^13;
alphaX = 10^-2;
gstar = 200;
HubbleI = Mphi;
Mpl = 1.22*10^19;

phiI = (3/(8*Pi))*Mpl^2/Mphi^2*HubbleI^2/Mphi^2;

Clear[Xeq, T];
T[x_, r_] := (30/(gstar*Pi^2))^(1/4)*Mphi*r^(1/4)/x;
Xeq[_, r_ /; r == 0] = 0.;
Xeq[x_?NumericQ, r_] :=
MX^3/Mphi^3*(2*Pi)^(-3/2)*x^3*(T[x, r]/MX)^(3/2)*Exp[-MX/T[x, r]];

c1 = Sqrt[3/(8*Pi)]*Mpl*alphaphi/Mphi;
c2 = c1*Mphi*alphaX/(MX*alphaphi);
c3 = c2*Mphi/MX;

eqns = {D[phi[x], x] == -c1*x*phi[x]/Sqrt[phi[x]*x + R[x]],
D[R[x], x] ==
c1*x^2*phi[x]/Sqrt[phi[x]*x + R[x]] +
c2*x^-1*(X[x]^2 - Xeq[x, R[x]]^2)/Sqrt[phi[x]*x + R[x]],
D[X[x], x] == -c3*
x^-2*(X[x]^2 - Xeq[x, R[x]]^2)/Sqrt[phi[x]*x + R[x]],
phi[1] == phiI, R[1] == 0, X[1] == 0};

sol = NDSolve[eqns, {phi, R, X}, {x, 1, 10^5},
Method -> "StiffnessSwitching"];


Plots, scaled to similar ranges:

Plot[
{phi[x]/phi[10^5], R[x]/R[10^5], X[x]/X[0.99 * 10^5]} /. sol // Evaluate,
{x, 1, 10^5}]


-