Using NDSolve with Manipulate when singularity or stiffness detected

Question

I want to solve a little bit complicated system of ODEs, where I have to vary initial conditions. My attempt to solution is like this (similar to what I found here Using ndsolve within manipulate):

ClearAll[w, x, y, t, k, x0, x1, y0, y1];
k = 2.5;
x0 = 5.5828*10^(-5);
x1 = 5.5828;
y0 = 1.05371*10^(-18);
y1 = 1.05371*10^(-19);
t0 = 6.7*10^(-5);
p = 3000;
sol[w_, z_] := sol[w, z] =
  NDSolve[{x'[t] ==-((x[t]/p)^(1/k) + x[t])(y[t] + 4 Pi t^3 x[t])/(t^2 - y[t] t), 
      y'[t] == 4 Pi t^2(x[t]/p)^(1/k), 
      x[t0] == w, 
      y[t0] == z}, {x, y}, {t, t0, 10}, MaxSteps -> 1000];
  Manipulate[Plot[({x[t] /. sol[w, z], y[t] /. sol[w, z]}), {t, t0, 10}],
  {w, x0,x1}, {z, y0, y1}]

First, when the initial conditions set to some constant numbers equation solves perfectly, although mathematica shows me that at some point stiffness or singularity detected (it is not shown in this picture because I changed all the units to natural one to make it easier to understand, though there is another error shown instead). The plot shown in the picture is exactly what I predicted Then when I tried to vary initial conditions using Manipulate, some crazy stuff comes out ( Error notifications shown by mathematica). I also tried to solve it using WhenEvent to take into account point of stiffness or singularity, but it doesn't lead to anywhere.

So now I am in thoughts how to solve my problem. I think error comes out, because of the stiffness or singularity detection, but I'm not sure.

P.s. I am relatively new to mathematica, so any help would be appreciated.

EDIT I also need that points of stiffness for further investigation, that is why I cannot use smaller limit in the case above. Actually, after overcoming this problem, I wanted to extract the points where singularity/stiffness detected for different initial conditions and put them in a List using WhenEvent. Although, I cannot even solve problem described above.

Answer 1

In addition to what Ulrich suggested, here is what works for me. The answers are similar in approach:

k = 2.5;
x0 = 5.5828*10^(-5);
x1 = 5.5828;
y0 = 1.05371*10^(-18);
y1 = 1.05371*10^(-19);
t0 = 6.7*10^(-5);
p = 3000;
ndsol[x0_, y0_] :=
  First@NDSolve[
    {x'[t] == -((x[t]/p)^(1/k) + x[t]) (y[t] + 4 Pi t^3 x[t])/(t^2 - 
          y[t] t),
     y'[t] == 4 Pi t^2 (x[t]/p)^(1/k),
     x[t0] == x0,
     y[t0] == y0}, {x, y}, {t, t0, 10}, MaxSteps -> 1000];
GetEqns[x0_, y0_] := {x[t], y[t]} /. ndsol[x0, y0]
Manipulate[{eqx, eqy} = Quiet@GetEqns[w, z]; 
 Plot[{eqx, eqy}, {t, t0, 10}], {w, x0, x1}, {z, y0, y1}]

You are facing problems because when you use Plot directly, it assigns constant values to the variables in equations, making it infeasible to be solved as diffrential equations. GetEqns is an intermediate function that gets all equations in symbolic form, which are easier to interpret for Plot.

Edit:

Here's a better solution:

ndsol[x0_, y0_] := Quiet@NDSolveValue[
    {x'[t] == -((x[t]/p)^(1/k) + x[t]) (y[t] + 4 Pi t^3 x[t])/(t^2 - 
          y[t] t),
     y'[t] == 4 Pi t^2 (x[t]/p)^(1/k),
     x[t0] == x0,
     y[t0] == y0}, {x, y}, {t, t0, 10}, MaxSteps -> 1000];
Manipulate[
 Plot[Evaluate[#[t] & /@ ndsol[x0, y0]], {t, t0, 10}], {w, x0, 
  x1}, {z, y0, y1}]

Answer 2

You should use the function ParametricNDSolve[] instead of NDSolve!

sol = ParametricNDSolve[{x'[t] == -((x[t]/p)^(1/k) + x[t]) (y[t] + 4 Pi t^3 x[t])/(t^2 - y[t] t), y'[t] == 4 Pi t^2 (x[t]/p)^(1/k), x[t0] == w, y[t0] == z},{x, y}, {t, t0, 5}, {w, z}]

solution sol={x,y} depends on the two parameters w,z

Plot[{x[x1, y0][t], y[x1, y0][t]} /. sol, {t, t0, 5}]

and might be further used in Manipulate.