Removing a singularity from a diferential equation I have the following ordinary differential equation for $y(x)\small x\in [0 \ h]$ where $h>0$
\begin{equation}
Ty''-(P_{0}+\gamma x)[1+y'^{2}]^{\frac{3}{2}} = 0; \quad
y(0) = 0;\quad \frac{1}{y'(0)}=0
\end{equation}
The ode comes from a force balance on the surface of a flexible cross-section. You can see the footprints of radius of curvature in the ode.
Anyways, I want a numerical solution of this equation given $h,\gamma,T,P_{0}$, all parameters $\geq0$. My strategy for the numerical solution is to isolate the singularity and then write $y(x) = singularity + z(x)$ or $y(x)=singularity \cdot z(x)$ near $x=0$. Then recast the problem in terms of $z(x)$. Then apply any RK4 scheme.
A further analysis near $x=0$ showed me that for non-zero $P_{0}$, we have
$singularity \sim iCx^{\frac{1}{2}}$ where $i=\sqrt(-1)$. And then I'm stuck with no good ideas. Recasting doesn't help even if one neglects $i$.
Is this approach correct? Is there any other way?
 A: Before trying to answer to the question, I would like to clarify a possible mess in the wording of the problem :
$$Ty''-(P_{0}+\gamma x) [1+y'^{2}]^{\frac{3}{2}} = 0; \quad
y(0) = 0;\quad y'(0)\rightarrow\infty \quad ;\quad 0\leq x \leq h \tag 1$$
The ODE is separable :
$$T\frac{y''}{[1+y'^{2}]^{\frac{3}{2}} }=P_0+\gamma x \tag 2$$
Integration leads to: 
$\quad T\int \frac{d(y')}{[1+y'^{2}]^{\frac{3}{2}} }=\int (P_0+\gamma x)dx \tag 3$
$$T\frac{y'}{[1+y'^{2}]^{\frac{1}{2}} }=P_0x+\frac{\gamma}{2}x^2+c \tag 4$$
$y'(0)\rightarrow\infty \quad\implies\quad \frac{y'}{[1+y'^{2}]^{\frac{1}{2}} }\rightarrow 1 \quad\implies\quad T=c \tag 5$
$$T\frac{y'}{[1+y'^{2}]^{\frac{1}{2}} }=P_0x+\frac{\gamma}{2}x^2+T \tag 6$$
$\text{All terms are positive.} \quad \begin{cases} 
T\frac{y'}{[1+y'^{2}]^{\frac{1}{2}} } <T \\
P_0x+\frac{\gamma}{2}x^2+T >T
\end{cases} \tag 7$
Thus the equality $(6)$ is impossible. Is there a typo in $(1)$ ?
IN ADDITION :
So, the problem of equation $(1)$ and given conditions has no real solution.
Nevertheless, just for curiosity, one can continue for complex solution.
$$T^2\frac{y'^2}{1+y'^2}=\left(P_0x+\frac{\gamma}{2}x^2+T\right)^2 \quad\to\quad  y'^2=\frac{\left(P_0x+\frac{\gamma}{2}x^2+T\right)^2}{T^2-\left(P_0x+\frac{\gamma}{2}x^2+T\right)^2} $$
$$y'=\pm\, i\frac{P_0x+\frac{\gamma}{2}x^2+T}{\sqrt{\left(P_0x+\frac{\gamma}{2}x^2+T\right)^2-T^2} }\qquad\qquad \left(P_0x+\frac{\gamma}{2}x^2+T\right)^2-T^2 \geq 0 $$
$$\boxed{ y(x)=\pm\, i\int_0^x \frac{P_0\xi+\frac{\gamma}{2}\xi^2+T}{\sqrt{\left(P_0\xi+\frac{\gamma}{2}\xi^2+T\right)^2-T^2}}d\xi }\qquad\text{so that}\quad y(0)=0$$
The integral could be expressed on closed form in terms of elliptic integral of second kind (Huge formula). In practice, numerical calculus of the integral is probably the simplest method. 
The integral is convergent for $x\to 0$ because $\frac{P_0\xi+\frac{\gamma}{2}\xi^2+T}{\sqrt{\left(P_0\xi+\frac{\gamma}{2}\xi^2+T\right)^2-T^2}} \simeq \frac{T}{\sqrt{2TP_0\xi}}$
$$y(x)\simeq \pm\, i\sqrt{\frac{2T}{P_0}}\sqrt{x}\qquad x\text{ small}$$
A: My previous post was not an answer to the question about the numerical method of solving the PDE. In fact this was a comment in order to draw attention to a possible typo in the PDE.
To be able to give a convenient answer to the question raised, one have to suppose that the PDE is corrected.
Close to $x=0$, the function $y(x)$ can be approximated to (equivalent) :
$$y(x)\sim c\:x^p \quad\text{with}\quad p>0\quad\text{so that}\quad y(0)=0$$
Hence
$$y'(x)\sim cpx^{p-1} \quad\text{then}\quad y'(0)=\infty \quad\implies \begin{cases} c>0\\p<1\end{cases}$$
All this implies 
$$0<p<1$$
$$y''(x)\sim c\,p\,(p-1)\,x^{p-2} \quad\implies\quad y''(x)<0$$
This is contradictory with the given equation :
$$Ty''=(P_{0}+\gamma x)[1+y'^{2}]^{\frac{3}{2}} \quad \text{where}\quad y''>0$$
So, probably there is a wrong positive/negative sign. 
Supposing that the corrected equation is :
$$Ty''+(P_{0}+\gamma x)[1+y'^{2}]^{\frac{3}{2}}=0$$
Then, your strategy for the numerical solution is correct : Isolate the singularity and then write $y(x) = \text{singularity} + z(x)$ or $y(x)=\text{singularity} \cdot z(x)$ near $x=0$. Then recast the problem in terms of $z(x)$. Then apply any RK4 scheme.
A further analysis near $x=0$ will show you that for non-zero $P_{0}$, we have
singularity $y(x)\simeq \sqrt{\frac{2T}{P_0}}\sqrt{x}$. So nothing prevents further numerical calculus.
Alternatively one can solve analytically the PDE :
$$y(x)=\int_0^x \frac{T-P_0\xi-\frac{\gamma}{2}\xi^2}{\sqrt{T^2-\left(T-P_0\xi-\frac{\gamma}{2}\xi^2\right)^2}}d\xi $$
A: After integration of the separable equation in $y'$, and using your initial condition on $y'$,
$$\frac{y'}{\sqrt{y'^2+1}}=ax^2+bx+1$$ (where $a,b$ are known).
This is enough to show that your equation cannot have a real solution with such a condition !

Now if you want to go ahead in the imaginaries, you have to integrate
$$y'=\frac{ax^2+bx+1}{\sqrt{(ax^2+bx+1)^2-1}}=\frac{ax^2+bx+1}{\sqrt{(ax^2+bx)(ax^2+bx+2)}}$$
($\pm i$ factor omitted !) and you can remove the singularity with
$$\left(y-\,\sqrt{\frac{2x}b}\right)'=\frac{ax^2+bx+1}{\sqrt{(ax^2+bx)(ax^2+bx+2)}}-\frac1{\sqrt{2bx}}.$$
