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