Ok, I guess it's time to solve the problem, but I think you already did a good job!
We start from your equation (3)
\frac{\dot{r}^2}{2}+\frac{h^2}{2 r^2} - \frac{G M}{r}=\text{const}.
Here h=L/m.
First we substitute r=1/u[\theta(t)]. Using the chain rule this gives
\dot{r}=-\frac{u'}{u^2} \dot{\theta},
where the prime refers to differentiation wrt. \theta and the dot to differentiation wrt. time.
Now use the angular-momentum conservation (it's also on your paper):
h=\frac{L}{m}=r^2 \dot{\theta}.
Subsitute this into the equation above, you get
\dot{r}=-\frac{u'}{u^2} \frac{h}{r^2}=-h u'.
Put this into the energy-conservation equation above
\frac{h^2 u'^2}{2} + \frac{h^2}{2} u^2-G M u=\text{const}.
The time derivative of this gives
h^2 u' u'' \dot{\theta}+h^2 u u' \dot{\theta}-G M u' \dot{\theta}=0.
Using again the angular-momentum conservation equation gives after some algebra
u''=\frac{G M}{h^2}-u.
This is the equation of an harmonic oscillator subject to an external constant force. The general solution obviously is
u(\theta)=\frac{G M}{h^2}+A \cos \theta + B \sin \theta.
Now it is the usual convention to choose the polar coordinate system such that \vartheta=0 denotes the perihel, i.e., the point, where the planet comes next to the sun. At this point u=1/r should be maximal, i.e., u'(0)=0 \; \Rightarrow \; B=0 and (assuming A>0) yields
u(\theta)=\frac{G M}{h^2} + A \cos \theta.
This gives
r(\theta)= \frac{1}{G M/h^2 + A \cos \theta}.
Put this a bit in a more familiar shape:
r(\theta)=\frac{p}{1+\epsilon \cos \theta},
where
p=\frac{h^2}{G M}, \quad \epsilon=\frac{h^2 A}{G M}.
One can show that, if the total energy is negative (bound motion), then \epsilon<1. This means that the planet moves along an ellipse with the sun in one of its foci. That's Kepler's first law.