I have a system of differential equations, given by
h0 = 69529/1000; og = 664/1000; c = 3*h0^2 (og); zfinal = 236/100; sigmas = 1;
C1 = -10;
V[z_] := c*(y[z])^C1 ;
Vp[z_] := c C1 y[z]^(-1 + C1);
Phi0 = 1;
Phi0p = 0;
n = 32; m = 64;
Klein = NDSolve[{H[z]^2 == h0 ^2 (1 - og) (1 + z)^3 + 1/3*( V[z] + 1/2*(1 + z)^2 H[z]^2*y'[z]^2 - H[z]^2 (1 + z)^2 y'[z]^2 ( 1/2*H'[z]^2 (1 + z)^2-Vp[z])) - H[z]^4 (1 + z)^3 y'[z]^3 , -1/ 3*(1 + z) D[ (V[z] + 1/2 H[z]^2 (1 + z)^2 y'[z]^2 - H[z]^2 (1 + z)^2 y'[z]^2*( 1/2*H'[z]^2 (1 + z)^2-Vp[z] )) - 3*H[z]^4 (1 + z)^3*y'[z]^3 , z] + H[z]^2 (1 + z)^2*
y'[z]^2 ( 1 - 2*( 1/2*H'[z]^2 (1 + z)^2-Vp[z] )) - 4 H[z]^2 (1 + z) y'[z] - ( -H[z]^2 (1 + z) y'[z] + H[z] (1 + z)*(y''[z]*H[z]*(1 + z) + y'[z] ( H'[z] (1 + z) + H[z]))) == 0, H[0] == h0, y[0] == Phi0, y'[0] == Phi0p}, {H, y}, {z, 0, zfinal}, PrecisionGoal -> n, AccuracyGoal -> n, WorkingPrecision -> m, Method -> "StiffnessSwitching"];
Plot[Evaluate[H[z] /. Klein], {z, 0, zfinal}, PlotRange -> All]
If I evaluate the code I get the following warnings:
NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.
NDSolve::nodae: The method NDSolve`StiffnessSwitching is not currently implemented to solve differential-algebraic equations. Use Method -> Automatic instead.
InterpolatingFunction::dmval: Input value {0.0000482114} lies outside the range of data in the interpolating function. Extrapolation will be used.
Is this normal? I'm doing something wrong?