Using NDSolve with Manipulate when singularity or stiffness detected

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.

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}]
• what is eqx and eqy in this case? and what is ndsol[x0,y0]?? Is it the same as sol[x0,y0]? – Ayan Nussupbekov Nov 30 '17 at 9:26
• They are intermediate variables to store the values returned by the function ndsol. sol and ndsol are refering to the same function. I corrected my answer to reflect this. – TheYeda Nov 30 '17 at 9:29
• does it mean that if I put sol instead of ndsol I will get the same answer? When I tried it using just sol I didnt get anything, but with ndsol I get desired result? Do you know why it is like this? Also, thanks, everything worked perfectly – Ayan Nussupbekov Nov 30 '17 at 9:33
• It was my typographic error. I meant to write it as ndsol, but ended up writing sol. Its sort of a variable name. You just need to use the same symbol at all places. If you want to use sol, use sol everywhere and likewise for ndsol – TheYeda Nov 30 '17 at 9:35

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.