<s>[INST] Is the choice for substitution correct for nonlinear second order differential equation? I am trying to understand why substitution destroys the information in original problem.  I have the following second order nonlinear differential equation,
\begin{equation}
\frac{d^2 Y}{d x^2} +  Y\sqrt{1 - Y^2} =0
 \label{eq:xdef}
\end{equation}
for Y(0) = C and Y(0)' = 0
When I plot the above equation using numerical solver, I get a cosine function.  I am trying to get rid of squareroot function in nonlinear equation.  I am trying to avoid using trigonometric functions or exponential function in substitution.  Therefore my attempt in substitution is to use:
\begin{equation}
z=\sqrt{1 - Y^2} 
 \label{}
\end{equation}
then the derivatives of the above substitution is,
\begin{equation}
\frac{dY}{dx}=-z\frac{dz}{dx}\frac{1}{Y}
 \label{}
\end{equation}
\begin{equation}
\frac{d^2 Y}{d x^2}=(-(\frac{dY}{dx})^2-(\frac{dz}{dx})^2-z\frac{d^2 z}{d x^2})\frac{1}{Y}
 \label{}
\end{equation}
then the nonlinear ODE  simplifies to
\begin{equation}
\frac{d^2 z}{d x^2}+(\frac{dz}{dx})^2\frac{1}{z(1-z^2)}-(1-z^2)=0
 \label{}
\end{equation}
\begin{equation}
z(0)=\sqrt{1 - C^2} 
 \label{}
\end{equation}
and z(0)' = 0
I think that the above results are correct.   However when I try to plot it using numerical solver, instead a cosine curve shifted upward by constant 1 I get either error message (singularity) or some other curve (nonoscillating).   Does that means that my substitution is invalid because information is lost? Or is my approach to plot it is incorrect?  Is there a better substitute that will eliminate squareroot without destroying the information?
 [/INST] You seem to have an error in your substitution formulas. Your first and second derivative formulas are correct. But the final equation in variable $z$ is not right. The right equation after substitution is
$$\frac{d^2z}{dx^2} + \frac{1}{z(1-z^2)}\left(\frac{dz}{dx}\right)^2 + z^2-1 = 0.$$
Is this substitution important or you need to find information on how to solve the equation? Because the original equation is a Hamiltonian system with one degree of freedom and as such is integrable, i.e. it is (more or less) explicitly solvable. One can write the system as a planar system (a system of two equations but involving only first derivatives)
$$\begin{align}
\frac{dY}{dx} &= V\\
\frac{dV}{dx} &= -Y\sqrt{1-Y^2}
\end{align}$$
Now, one way to go about this is to write the system as 
$$\frac{dY}{V} = dx = - \frac{dV}{Y\sqrt{1-Y^2}}$$
$$\frac{dY}{V} = - \frac{dV}{Y\sqrt{1-Y^2}}$$ which after cross-multiplying leads to 
$$Y\sqrt{1-Y^2}dY = -VdV$$
$$VdV + Y\sqrt{1-Y^2}dY = 0.$$
Integrating the last identity on both sides leads to 
$$H(Y,V) = \frac{V^2}{2} + \int Y\sqrt{1-Y^2} dY = \frac{V^2}{2} - \frac{1}{3}(1-Y^2)^{\frac{3}{2}}=E_0$$
where $E_0$ is a constant (energy level). This means that a solution $(Y(x),V(x))$ to the system leaves the function $H$ invariant, i.e. $H(Y(x),V(x)) = E_0$ for all $x$. $H$ is the total energy of the system, where $\frac{V^2}{2}$ is the kinetic energy, while the therm $U(Y) = -\frac{1}{3}(1-Y^2)^{3/2}$ is the potential energy. To understand the trajectories of the system in the $(Y,V)$ plane, one can draw a graph of the potential $U(Y)$ and based on it one can reconstruct the dynamics of the system. On the other hand, since $V = \frac{dY}{dx}$, then the Hamiltonian gives us the equation
$$\left(\frac{dY}{dx}\right)^2 = 2E_0 + \frac{2}{3}(1-Y^2)^{\frac{3}{2}}$$ which turns into
$$\frac{dY}{dx} = \pm \sqrt{2E_0 + \frac{2}{3}(1-Y^2)^{\frac{3}{2}}}$$ and so you can write it as
$$\frac{dY}{ \sqrt{2E_0 + \frac{2}{3}(1-Y^2)^{\frac{3}{2}}}} = \pm dx$$ and by integrating both sides, one gets
$$\int_{Y_0}^{Y(x)} \frac{d\tilde{Y}}{ \sqrt{2E_0 + \frac{2}{3}(1-\tilde{Y}^2)^{\frac{3}{2}}}} = \pm(x-x_0),$$
so basically you can try to solve the integral 
$$\int \frac{d{Y}}{ \sqrt{2E_0 + \frac{2}{3}(1-{Y}^2)^{\frac{3}{2}}}},$$ although the solutions are probably restricted to the real line holomorphic functions that uniformize an algebraic curve of genus 2 (a hyper-elliptic Riemann surface of genus 2), which is a bit tough.
Alternativly, you can try to run some simulations of the equation
$$\frac{dY}{dx} = \pm \sqrt{2E_0 + \frac{2}{3}(1-Y^2)^{\frac{3}{2}}},$$ having in mind that whenever the solution's derivative $Y'(x)$ becomes $0$, the sign of the equation switches (the sign $\pm$ in front of the square root). The behavior, for certain energy levels is oscillatory, kind of like sine and cosine, but definitely not exactly like them. And do not forget $-1 \leq Y \leq 1$. 
Now if one performs the change $z=\sqrt{1-Y^2}$, then this is a type of canonical transformation of the Hamilton system, so the Hamiltonian of the old system, in terms of $(Y,V)$ turns into the Hamiltonian of the new system in terms of $(z,w)$ (which I wrote first), where the 2D transformation is 
\begin{align}
z &= \sqrt{1-Y^2}\\
w &= - \frac{Y}{\sqrt{1-Y^2}}V \left(= \frac{dz}{dx}\right) 
\end{align}
and the inverse is
\begin{align}
Y &= \sqrt{1-z^2}\\
V &= - \frac{z}{\sqrt{1-z^2}}w \left(= - \frac{z}{\sqrt{1-z^2}} \frac{dz}{dx}\right) 
\end{align} 
Consequently,
$$\frac{dY}{dx} = - \frac{z}{\sqrt{1-z^2}} \frac{dz}{dx}$$ and when we substitute it in the first-order equation (which would be equivalent to differentiating the substitution one more time, substituting it in the original equation and getting the equation I wrote first):
$$-\frac{z}{\sqrt{1-z^2}}\frac{dz}{dx} = \pm \sqrt{2E_0 + \frac{2}{3} z^3}.$$  If we square on both sides we get 
$$\left(z\frac{dz}{dx}\right)^2 = \frac{2}{3}(1-z)(1+z)(3E_0 + z^3)$$ where $z=z(x), w = \frac{dz}{dx}(x)$ is a parametrization (after allowing $x \in \mathbb{C}$) of the hyper-elliptic Riemann surface
$$z^2 w^2 = \frac{2}{3}(1-z)(1+z)(3E_0 + z^3).$$
The general solution to the equation in terms of $z$ is then
$$ \int_{z_0}^{z(x)} \frac{\tilde{z} d\tilde{z}}{\sqrt{(1-\tilde{z})(1+\tilde{z})(3E_0 + \tilde{z}^3)}} = \pm \sqrt{\frac{2}{3}}(x-x_0).$$ 
</s>