'Power::infy' and 'Infinity::indet:' error when solving systems of ODE

Question

I am new to mathematica and I am trying to solve the systems of ODE as shown below

Clear["Global`*"]
(*constants*)
phi1 = 1;
s = 0.5;
P = 1;
Cref = 1;
Q = 1;
lambda = 1;

(*the system of ode*)
ode = {x*a''[x] - 2*a'[x] - phi1^2*(1 + s)*phi[x]*a[x] *x == 0,
    x*b''[x] - 2*b'[x] - phi1^2/P*phi[x]*a[x]*x == 0,
   x*phi''[x] - 2*phi'[x] - s*phi1^2*Cref/Q/lambda*phi[x]*a[x]*x == 0};
bcs = {a'[0] == 0, a[1] == 0.9, b'[0] == 0, b[1] == 0.95, 
   phi'[0] == 0, phi[1] == (lambda - 0.6*Cref)/lambda};

(*ndsolve*)
{asol, bsol, phisol} = 
  NDSolveValue[{ode, bcs}, {a, b, phi}, {x, 0, 1}];
Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"FiniteElement"}};
Plot[{asol[x], bsol[x], phisol[x]}, {x, 0, 1}, 
 PlotLegends -> "Expressions"]

I would get error 'Power::infy' and 'Infinity::indet' when i am solving the system from x = 0 to 1. However, this can be overcome when i replace 0 with a close to 0 value. May i ask if there is other better solution for this?

Power::infy: Infinite expression 1/0. encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Answer 1

When this happens means you have issue with starting at zero, where first step gives 1/0 (since you are starting from x=0).

Without looking more at it to find why, and which term, a quick workaround is to start from little after x=0. Using $MachineEpsilon actually worked.

Clear["Global`*"]
(*constants*)
phi1 = 1;
s = 0.5;
P = 1;
Cref = 1;
Q = 1;
lambda = 1;
del = $MachineEpsilon; (*added this *)

(*the system of ode*)
ode1 = x*a''[x] - 2*a'[x] - phi1^2*(1 + s)*phi[x]*a[x]*x == 0
ode2 = x*b''[x] - 2*b'[x] - phi1^2/P*phi[x]*a[x]*x == 0
ode3 = x*phi''[x] - 2*phi'[x] - s*phi1^2*Cref/Q/lambda*phi[x]*a[x]*x ==
   0
bcs = {a'[del] == 0, a[1] == 0.9, b'[del] == 0, b[1] == 0.95, 
  phi'[del] == 0, phi[1] == (lambda - 0.6*Cref)/lambda}

(*ndsolve*)
NDSolveValue[{ode1, ode2, ode3, bcs}, {a, b, phi}, {x, del, 1}]

No errors now.

Again, the issue is that there is a term or derivative in your equations which gives 1/x when x=0 this gives a problem. So more investigation is need to find it.