Solve an ODE involving an Interpolating function with Mathematica

Question

Using Mathematica, I want to solve an ODE of an Interpolating Function. First I solve the ODE obtaining the function y[s] (or Y[s])

aa = .4; sol =  NDSolve[{-4 + 4/y[s]^2 - y[s]^2 + y[s]^4 + y'[s]^2 ==  0, y[0] == 1.207}, y[s], {s, 0, aa}];
ParametricPlot [Evaluate[{s, y[s]} /. sol[[2]]], {s, 0, aa},  PlotRange -> All]
Y[s_] = y[s] /. sol[[2]];

(with the function Y[s] I can plot the solution y[s])

My problem is solving the ODE w'[s]= Sqrt[y[s]^2 + 2] with initial condition w[0] = 0. My first attempt has been

w1 = NDSolve[{w'[s] == Sqrt[Y[s]^2 + 2], w[0] == 0}, w[s], {s, 0, aa}]

Next, I put W[s_] = w[s] /. w1[[1]] and I plot the function W[s] by means of Plot[W[s], {s, 0, aa}]. However, I think that it is wrong because the plot of W[s] seems to be a straight-line. In fact, I hope that Mathematica solves the ODE for w[s] by taking only the first value for y[s], namely, y[0]=1.207 and then solves the ODE w'[s]=Sqrt[1.207^2+2], w[0]=0 (whose solution is a linear function). I have computed the last ODE:

p1 = NDSolve[{p'[s] == Sqrt[1.207^2 + 2], p[0] == 0},   p[s], {s, 0, aa}]

and the plots of both functions p[s] and w[s] are 'very similar'. Please, help.

Answer 1

The first equation has two branches of the solution, so you need to select a branch, define a function, and then solve the second equation. The solutions of the second equation can be compared for each branch with a linear function. It is obvious that the solutions differ among themselves and differ from the linear function.

aa = .417527; sol = 
 NDSolve[{-4 + 4/y[s]^2 - y[s]^2 + y[s]^4 + y'[s]^2 == 0, 
   y[0] == 1.207}, y, {s, 0, aa}];
Y1[s_] := y[s] /. First[sol]
Y2[s_] := y[s] /. Last[sol]
sol1 = NDSolveValue[{w'[s] == Sqrt[2 + Y1[s]^2], w[0] == 0}, 
   w, {s, 0, aa}];
sol2 = NDSolveValue[{w'[s] == Sqrt[2 + Y2[s]^2], w[0] == 0}, 
   w, {s, 0, aa}];
sol3 = NDSolveValue[{w'[s] == Sqrt[2 + 1.207^2], w[0] == 0}, 
   w, {s, 0, aa}];
{Plot[{Y1[x], Y2[x]}, {x, 0, aa}, PlotLegends -> Automatic, 
  PlotRange -> All], 
 Plot[{sol1[x], sol2[x], sol3[x]}, {x, 0, aa}, 
  PlotLegends -> Automatic, PlotRange -> All]}