NDSolve for coupled differential equations with large coefficient

Question

I have three coupled differential equations to solve. Those equations are

equation1=neta'[t] == -(7.6458*10^-23) neta[t] - (1.49912 neta[t])/t


equation2=rho'[t] == -(3.89898*10^36/t^3) + (3.97887*10^-23) neta[t]


equation3=npsi'[t] == (9 npsi[t])/(78125000000000000000000 \[Pi]) - (1.49912 npsi[t])/t

I have tried to solve these equations by using

soln = NDSolve[{equation1,equation2, equation3, neta[10^-10] == 0.0025,rho[10^-10]== 4*10^6 ,npsi[10^-10] == 0}, {neta,rho,npsi}, {t, 10^-10, 1}, AccuracyGoal -> Autometic, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, 
"StiffnessTest" -> True}, MaxStepSize -> Autometic, PrecisionGoal -> Autometic]

It is showing the error

NDSolve::ndsz: At t == 1.`*^-10, step size is effectively zero;singularity or stiff system suspected.

However, if I just skip the second equation and solve for the other two then it can be solved by the same code. Even if I skip the first part of equation2 then also it is solving. So I think the problem is in the first part of the equation2, a very large number is coming there and it fails to solve.

Can you please help in this regards?

Thank you.

---

Suppose I have a set of equations

eqn1=neta'[t] == -7.6458*10^-23 neta[t] - 3.566*10^-18 neta[t] Sqrt[rho[t]]

eqn2=rho'[t] == 3.97887*10^-23 neta[t] - 4.75467*10^-18 rho[t]^(3/2)

eqn3=npsi'[t] == 3.66693*10^-23 neta[t] - 3.566*10^-18 npsi[t] Sqrt[rho[t]]

Then if I do it in a similar way then

soln = NDSolve[{eqn1,eqn2, eqn3, neta[10^-10] == 0.0025,rho[10^-10]== 4*10^6 ,npsi[10^-10] == 0}, {neta,rho,npsi}, {t, 10^-10, 1}, AccuracyGoal -> Autometic, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, "StiffnessTest" -> True}, MaxStepSize -> Autometic, PrecisionGoal -> Autometic]

It shows same kind of error

NDSolve::ndtol: Tolerances requested by the AccuracyGoal and PrecisionGoal options could not be achieved at t == 1.`*^-10.

Answer 1

We write the system of equations in the form

equation1 = neta'[t] == a11*neta[t] + a12*neta[t]/t;


equation2 = rho'[t] == a21/t^3 + a22* neta[t];


equation3 = npsi'[t] == a31* npsi[t] + a32 npsi[t]/t;

The exact solution allows us to see the source of the problem.

soln = DSolve[{equation1, equation2, equation3}, {neta, rho, npsi}, t]

(*Out[]= {{neta -> Function[{t}, E^(a11 t + a12 Log[t]) C[1]], 
  rho -> Function[{t}, -(a21/(2 t^2)) + C[2] - 
     a22 t^(1 + a12) (-a11 t)^(-1 - a12) C[1] Gamma[1 + a12, -a11 t]],
   npsi -> Function[{t}, E^(a31 t + a32 Log[t]) C[3]]}}*)

 netaC = First[neta[t] /. soln]; rhoC = 
 First[rho[t] /. soln]; npsiC = First[npsi[t] /. soln];
eq1 = netaC /. t -> t0; eq2 = rhoC /. t -> t0; eq3 = npsiC /. t -> t0;
sol = Solve[
   eq1 == neta0 && eq2 == rho0 && eq3 == npsi0, {C[1], C[2], C[3]}];
netaT = netaC /. sol; rhoT = rhoC /. sol; npsiT = npsiC /. sol;


a11 = -(7.6458*10^-23); a12 = -1.49912; a21 = -3.89898*10^36; a22 = 
 3.97887*10^-23; a31 = 
 9./(78125000000000000000000 \[Pi]); a32 = -1.49912; neta0 = 0.0025; \
rho0 = 4.*10^6; npsi0 = 0.; t0 = 1.*10^-10; 
{netaT, rhoT, npsiT}
{{(2.55117*10^-18 E^(7.6458*10^-33 - 7.6458*10^-23 t))/
  t^1.49912}, {(-1.94949*10^56 + 2.44065*10^-51 I) + 1.94949*10^36/
   t^2 - 9.28269*10^-52 Gamma[-0.49912, 7.6458*10^-23 t]}, {0.}}

To get a numerical solution, we need to increase WorkingPrecision -> 100and all coefficients and initial data set with this precision.

a11 = -7.6458`100*10^-23; a12 = -1.49912`100; a21 = \
-3.89898`100*10^36; a22 = 3.97887`100*10^-23; a31 = 
 SetPrecision[9/(78125000000000000000000 \[Pi]), 
  100]; a32 = -1.49912`100; neta0 = 0.0025`100; rho0 = 
 4.0`100*10^6; npsi0 = 0; t0 = 1.0`100*10^-10;
equation1 = neta'[t] == a11*neta[t] + a12*neta[t]/t;
equation2 = rho'[t] == a21/t^3 + a22* neta[t];
equation3 = npsi'[t] == a31* npsi[t] + a32 npsi[t]/t;
soln = NDSolve[{equation1, equation2, equation3, neta[t0] == neta0, 
   rho[t0] == rho0, npsi[t0] == npsi0}, {neta, rho, npsi}, {t, t0, 
   1.`100}, WorkingPrecision -> 100, MaxSteps -> 10^6]


{LogLogPlot[neta[t] /. soln, {t, t0, 1}, PlotRange -> All, 
  PlotLabel -> "neta"], 
 LogLogPlot[Abs[rho[t]] /. soln, {t, t0, 1}, PlotRange -> All, 
  PlotLabel -> "rho"]}