# NDSolve with a 0/0 error

I am trying to solve this set of coupled first order differential equations. The $$\phi'$$ term gives a 0/0 error for when $$\epsilon=0$$. When I enter:

NDSolve[{ϕ'[ϵ] == -1/ϵ^2 ((1 + (8 c ϕ[ϵ])/αm[M]) (Mhat[ϵ] +
c αm[M]^3 ϵ^3 ghat[ϕ[ϵ]/αm[M]]))/(
1 - ((16 c)/αm[M]) (Mhat[ϵ]/ϵ)),
Mhat'[ϵ] == αm[M]^3 ϵ^2 (((ϕ[ϵ]/αm[M])^2 - 1)^(3/2) + 3 c fhat[ϕ[ϵ]/αm[M]]),
Mhat[0] == 0, ϕ[0] == 1, ϕ'[0] == 0}, {Mhat, ϕ}, {ϵ, 0, 1}]


Where:

ghat[y_] := 3 ArcSinh[Sqrt[y^2 - 1]] + y (2 y^2 - 5) Sqrt[y^2 - 1]
fhat[y_] := y (2 y^2 - 1) Sqrt[y^2 - 1] - ArcSinh[Sqrt[y^2 - 1]]


and $$\alpha m[M]$$ is just a constant.

I get

Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
Power::infy: Infinite expression 1/0.^2 encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at \[Epsilon] == 0..


I have tried adding the condition $$\phi'[0]==0$$, but then I get:

NDSolve::icordinit: The initial values for all the dependent variables are not explicitly specified.
NDSolve will attempt to find consistent initial conditions for all the variables.
NDSolve::ndnco: The number of constraints (3) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (2).


I have also tried replacing the $$\phi'$$ expression with a piecewise function that is the same as the above expression for $$\phi'$$ except that it is set to zero at $$\epsilon=0$$. This gives me a solution that is only valid at $$\epsilon=0$$.

I know this can be solved numerically, I just can't figure out how to get mathematica to do it.

• There are several undefined functions. Can't really help beyond Bill's observation that epsilon = 0 is a singular point. Nov 6, 2020 at 4:18
• Unfortunately I still get the same 1/0 type error when I rearrange both sides to remove the explicit division by zero terms. I have edited the above description to include the undefined terms. Nov 6, 2020 at 12:58
• I don't think multiplying by epsilon will work because the first thing NDSolve does is to solve for the derivatives. Multiplying both sides of an equation won't change that. We need your complete code, because dealing with singular points is tricky. Nov 6, 2020 at 15:39
• I agree that these equations can be solved, but not without values for c and αm[M]. Please provide. Nov 6, 2020 at 19:46

This system of equations can be solved as follows. For convenience, define

eq1 = ϕ'[ϵ] == -1/ϵ^2 ((1 + (8 c ϕ[ϵ])/αm[M]) (Mhat[ϵ] + c αm[M]^3 ϵ^3 ghat[ϕ[ϵ]
/αm[M]]))/(1 - ((16 c)/αm[M]) (Mhat[ϵ]/ϵ));
eq2 = Mhat'[ϵ] == αm[M]^3 ϵ^2 (((ϕ[ϵ]/αm[M])^2 - 1)^(3/2) + 3 c fhat[ϕ[ϵ]/αm[M]]);


Because the system of equations is singular at ϵ = 0, expand the right side of eq2 there.

Normal@Series[Last[eq2], {ϵ, 0, 2}] /. ϕ[0] -> 1
(* ϵ^2 αm[M]^3 (((1 - αm[M]^2)/αm[M]^2)^(3/2) + 3 c (-ArcSinh[Sqrt[-1 + 1/αm[M]^2]] +
((-1 + 2/αm[M]^2) Sqrt[(1 - αm[M]^2)/αm[M]^2])/αm[M])) *)


It follows that MHat varies as ϵ^3 near ϵ = 0, so make the substitution, Mhat[ϵ] -> ϵ^3 Nhat[ϵ] and expand the right side of eq1.

Normal@Series[Last[eq1 /. Mhat[ϵ] -> ϵ^3 Nhat[ϵ]], {ϵ, 0, 1}] /. ϕ[0] -> 1
(* -ϵ (1 + (8 c)/αm[M]) (Nhat[0] + c αm[M]^3 (3 ArcSinh[Sqrt[-1 + 1/αm[M]^2]] +
((-5 + 2/αm[M]^2) Sqrt[(1 - αm[M]^2)/αm[M]^2])/αm[M])) *)


which now is well behaved and also consistent with the OP's desire that ϕ'[0] == 0. So, make the Nhat substitution into the two equations.

eq1r = Simplify[eq1 /. Mhat[ϵ] -> ϵ^3 Nhat[ϵ]]
(* ϕ'[ϵ] == 1/(16 c ϵ^2 Nhat[ϵ] - αm[M]) ϵ (αm[M] + 8 c ϕ[ϵ]) (Nhat[ϵ] +
c (3 ArcSinh[Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2]] αm[M]^3 - 5 αm[M]^2 ϕ[ϵ]
Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2] + 2 ϕ[ϵ]^3 Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2])) *)

eq2r = Simplify[((#/ϵ^2) & /@ eq2) /. Mhat -> Function[{ϵ}, ϵ^3 Nhat[ϵ]]]
(* 3 Nhat[ϵ] + 3 c ArcSinh[Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2]] αm[M]^3 +
3 c ϕ[ϵ] (αm[M]^2 - 2 ϕ[ϵ]^2) Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2] +
ϵ Nhat'[ϵ] == αm[M]^3 (-1 + ϕ[ϵ]^2/αm[M]^2)^(3/2) *)


In order that eq2r be well behaved at ϵ = 0, the boundary condition Nhat[0] must be chosen such that eq2r reduces to 0 == 0 there.

Nhat0 = Nhat[0] /. Flatten@Solve[eq2r /. ϵ -> 0 /. ϕ[0] -> 1, Nhat[0]]
(* 1/3 (6 c Sqrt[-1 + 1/αm[M]^2] + Sqrt[-1 + 1/αm[M]^2] αm[M] -
3 c Sqrt[-1 + 1/αm[M]^2] αm[M]^2 - 3 c ArcSinh[Sqrt[-1 + 1/αm[M]^2]] αm[M]^3 -
Sqrt[-1 + 1/αm[M]^2] αm[M]^3) *)


With these transformations the ODEs can be solved by

ϵ0 = 10^-6;
s = NDSolveValue[{eq1r, eq2r, ϕ[ϵ0] == 1, Nhat[ϵ0] == Nhat0} /. {c -> 1, αm[M] -> 1/6},
{ϕ[ϵ], Nhat[ϵ]}, {ϵ, ϵ0, 1}];
Plot[First@s, {ϵ, ϵ0, 1}, PlotRange -> {0, 1}, AxesLabel -> {ϵ, ϕ},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
Plot[ϵ^3 Last@s, {ϵ, ϵ0, 1}, AxesLabel -> {ϵ, Mhat}, ImageSize -> Large,
LabelStyle -> {15, Bold, Black}]


In the absence of definitions in the Question for the two constants, I chose c = 1 arbitrarily, and αm[M] slightly smaller than ϕ[1], so that -1 + ϕ[ϵ]^2/αm[M]^2 > 0`.