Now I know that $dy/dx$ is not a fraction etc., but if you look at the works of 17th -- early 19th century mathematicians you'll see they played fast & loose with their infinitesimals and arrived at correct results.
I tried to give that a go on a differential equation, and I think I obtained an incorrect result. My problem is I can't give an argument as to why any of the steps are illegal. Can anyone help me? Here's the derivation:
A point mass $M$ fixed at the origin gravitationally attracts another point mass. The motion of the test mass satisfies:
$$ \ddot{\boldsymbol r} = -\frac{gM\hat{\boldsymbol r}}{r^2} \tag{1}\label{1}$$
The acceleration vector in polar coordinates is given by
$$ \ddot{\boldsymbol r} = (\ddot r - r\dot\varphi^2)\hat{\boldsymbol r} + \frac{1}{r}\frac{d}{dt}(r^2\dot\varphi)\hat{\boldsymbol\varphi} $$
Plugging that in \eqref{1} and equating like components we get that the $\hat{\boldsymbol\varphi}$ component is zero, which implies
$$ r^2\frac{d\varphi}{dt} = k \tag{2}\label{2} $$
a constant. So we only have the radial equation:
$$ \frac{d^2r}{dt^2} - r\left(\frac{d\varphi}{dt}\right)^2 = -\frac{gM}{r^2} $$
Let's try and prove that this is a conic. We'll use \eqref{2}, rewritten as
$$ \frac{r^2\,d\varphi}{k} = dt $$
to get rid of the $dt$'s in the equation. Also, as long as we're playing "$dy/dx$ is a fraction", we may as well go ahead and expand that squared derivative:
$$ \frac{d^2r}{dt^2} - \frac{r\,d\varphi^2}{dt^2} = -\frac{gM}{r^2} $$
Now rewriting $dt^2$ as $r^4\,d\varphi^2/k^2$ in the above,
$$ \frac{k^2\,d^2r}{r^4\,d\varphi^2} - \frac{k^2r\,d\varphi^2}{r^4\,d\varphi^2} = -\frac{gM}{r^2} $$
$$ \frac{d^2r}{r^4\,d\varphi^2} - \frac{r\,d\varphi^2}{r^4\,d\varphi^2} = -\frac{gM}{k^2r^2} $$
$$ \frac{d^2r}{d\varphi^2} - \frac{r\,d\varphi^2}{\,d\varphi^2} = -\frac{gMr^2}{k^2} $$
$$ \frac{d^2r}{d\varphi^2} = -\frac{gMr^2}{k^2} + r $$
Now, this equation already seems to be wrong. It still holds for a circular orbit (don't forget that $k=r^2\dot\varphi$ though), but I don't think it holds for other conics. Moreover, I worked a bit more on this equation (now with more "conservative" steps) and in more than one way arrived at equations that IIRC the general solutions of which are in terms of elliptic functions (and didn't seem to necessarily reduce to something simpler). See:
(1) If the solution is a conic then
$$ r = \frac{ep}{1-e\cos\varphi} $$
and thus
$$ \frac{d^2r}{d\varphi^2} = \frac{2e^3p-e^2p\cos\varphi-e^3p\cos^2\varphi}{(1-e\cos\varphi)^3} $$
but, abbreviating $gM/k^2$ with $a$,
\begin{align} -ar^2 + r &= -\frac{ae^2p^2}{(1-e\cos\varphi)^2} + \frac{ep}{1-e\cos\varphi} \\ &= \frac{-ae^2p^2(1-e\cos\varphi)+ep(1-e\cos\varphi)^2}{(1-e\cos\varphi)^3} \\ &= \frac{-ae^2p^2+ae^3p^2\cos\varphi+ep+e^3p\cos^2\varphi-2e^2p\cos\varphi}{(1-e\cos\varphi)^3} \\ &= \frac{(ep-ae^2p^2)+(ae^3p^2-2e^2p)\cos\varphi+e^3p\cos^2\varphi}{(1-e\cos\varphi)^3} \end{align}
Now we must have these two equal. But the coefficients of the cosine-squared terms can only match if $e=0$ or $p=0$. So general conics don't satisfy the eq.
(2) In $r''=-ar^2 + r$, complete the square on the RHS to arrive at something like $r''=-A\rho^2 + B$. We'll have $d\rho/dr=1$ so the eq. is
$$ \rho'' = A\rho^2 + B $$
Multiply by $\rho'$ on both sides and integrate,
$$ \frac{1}{2}\rho'^2 = 2A\rho^3 + B\rho + C $$
Separate the variables:
$$ \frac{d\rho}{\sqrt{4A\rho^3 + 2B\rho + C}} = d\varphi $$
Integrating again the LHS becomes an elliptic integral in general. So, incorrect.
Where is the illegal step and why is it illegal? Or is the mistake somewhere else?
Special thanks to OP who helps me a lot with this.
Your radical substitution where the differentials are treated as fractions are correct, when $r$ is treated as a constant. You can actually check by more conservative methods. $$ \frac{d^2 r}{dt^2}=\frac{d}{dt}\left(\frac{d r}{dt}\right)=\frac{d}{dt}\left(\frac{d r}{d\varphi}\cdot\frac{d\varphi}{dt} \right)=\frac{d}{dt}\left(\frac{d r}{d\varphi}\cdot\frac{k}{r^2} \right)=\frac{k}{r^2}\frac{d}{d\varphi}\left(\frac{d r}{d\varphi} \right)\cdot \frac{d\varphi}{dt}=\frac{k^2}{r^4}\frac{d^2r}{d\varphi^2}. $$
However, in general, we have
$$ \frac{d^2 r}{dt^2}=\frac{d}{dt}\left(\frac{d r}{dt}\right)=\frac{d}{dt}\left(\frac{d r}{d\varphi}\cdot\frac{d\varphi}{dt} \right)=\frac{d}{dt}\left(\frac{d r}{d\varphi}\cdot\frac{k}{r^2} \right)=\frac{k}{r^2}\frac{d}{d\varphi}\left(\frac{d r}{d\varphi} \right)\cdot \frac{d\varphi}{dt}-\frac{2k}{r^3}\frac{dr}{d\varphi}=\frac{k^2}{r^4}\frac{d^2r}{d\varphi^2}-\frac{2k}{r^3}\frac{dr}{d\varphi}. $$ The second term of this does not show up in your radical method, so the final result is incorrect.
However, $dr/d\varphi$ is zero when the orbit is a circle, so your equation is satisfied by a circle.
Your tangential components of the acceleration seems to be wrong $$ \frac{d^2\mathbf{r}}{dt^2}=(\ddot r - r\dot\varphi^2)\hat{\boldsymbol r}+(r\ddot \theta + 2\dot r \dot\varphi)\hat{\boldsymbol \varphi}. $$ For a proof, follow this Link.
In your equation, you write $\frac{1}{4}\frac{d}{dt}(r^2\dot\varphi)\hat{\boldsymbol\varphi}=\frac{1}{4}r(r\ddot \theta + 2\dot r \dot\varphi)\hat{\boldsymbol \varphi}$ instead, which differs from the first equation by a factor of $(1/4)r$. However, that does not affect your results, since they are equavalent when set to zero, in your situation.