Differential Equation of orbit of a planet Every where I found the solution of the differential equation orbit ( they make the differential equation in the 1st place )  by letting r=1/u  , using reciprocal coordinates.
I was trying  to do away with this substitution but I couldn't form the equation.
I was having a problem in whether (r) is a function of theta which in turn is a function of t or is r directly a function of t. Or is theta the function of r and then r is a function of t or theta directly a function of t. I tried this but produced a mess
I need help in forming the equation without using u , just r and in understanding what is function of what both implicitly and explicitly !
Edit:
After searching on the net I found out that without taking the substitution , I had actually formed the correct equation but that equation is non linear ! I got that correct non linear equation. To solve the equation we use the substitution r=1/u which will linearize it. So it's just linearization !
What I want is that please someone check 
This statement 
"whether (r) is a function of theta which in turn is a function of t or is r directly a function of t. Or is theta the function of r and then r is a function of t or theta directly a function of t"
What is correct in the above statement ? I am getting confused because if theta is function of r then the differential equation would be different . Would it be so that we just replace r and theta in final answer and get the same answer back ? 
 A: Polar coordinates
The movement vector $\vec{OP}=\vec r$ in polar coordinates depends of $\theta$ which in turn depends of $t$, so $r=r(\theta(t))$.
But we will see that if $\theta$ is a convenient angle to describe the trajectory curve it is not straightforward to describe how $P$ evolves in time (or equivalently, to have a formula for $\theta(t)$). Instead we will introduce a new "angle" $E$ named eccentric anomaly and we will see that there are formulas to relate $E\leftrightarrow\theta$ and also to relate $E\leftrightarrow t$. 
Thus giving a way via $E$ intermediate to express both $r$ and $\theta$ as direct functions of time $t$.

Central force movement
gravity law : $m\vec{\ddot r}=\vec F=-\frac{\mathcal GmM}{r^3}\vec r\quad$ where $m$ is planet mass and $M$ is sun mass.
moment law : $\vec\sigma_O=\vec{OP}\wedge m\vec v=\vec r\wedge m\vec{\dot r}$
When derivating we get $\frac{\partial\vec\sigma_O}{\partial t}=\vec{\dot r}\wedge m\vec{\dot r}+\vec r\wedge m\vec{\ddot r}=0+\vec r\wedge\vec F=\vec 0\quad$ since $F$ is a radial force.
Thus $\vec\sigma_O$ is a constant $\vec C$, we name it the area constant.
$\bbox[5px,border:1px solid] {\vec C=\vec r\wedge\vec{\dot r}}$
In the frame $(\vec i,\vec j,\vec k)$ defined by $\vec i=\frac{\vec r}{r},\ \vec j=\frac{\partial\vec i}{\partial\theta},\ \vec k=\vec i\wedge\vec j$.
Note: $\vec i=(\cos\theta,\sin\theta,0)$ so $(\vec i,\vec j,\vec k)$ is a well defined orthonormal frame.
Let's multiply by $\vec v=\vec{\dot r}$ and calculate the kinetic energy : $m\,\vec{\ddot r}.\vec{\dot r}=-\mathcal GmM\,\frac{\vec r.\vec{\dot r}}{r^3}\ $ 
integrating this gives $\frac12mv^2=+\mathcal GmM\frac 1r+cst$ 
We call the arbitrary integration constant $h=\frac{cst}{m}$ and also $\bbox[5px,border:1px solid]{\mu=\mathcal GM}$
$\vec{\dot r}=\frac{\partial(r\vec i)}{\partial t}=\dot r\,\vec i+r\dot\theta\,\vec j\quad$ from this we get $\vec C=r^2\dot\theta\,\vec k\ $ and $\ v^2={\dot r}^2+r^2{\dot\theta}^2$
$\bbox[5px,border:1px solid] {\begin{array}{cc} C=r^2\dot\theta\\ \frac12({\dot r}^2+r^2{\dot\theta}^2)-\frac{\mu}{r}=h \end{array}}$

Conic trajectory
This is where the manipulation is a bit artificial, we express $\dot r=\dot\theta\frac{dr}{d\theta}$ in terms of $\dot\theta=\frac{C}{r^2}$
giving $\frac12(\frac{C^2}{r^4}(\frac{dr}{d\theta})^2+\frac{C^2}{r^2})-\frac{\mu}{r}=h,\quad$and we do the change of variable $\frac1r=u+\frac{\mu}{C^2}$.
I do not detail the calculacy but this gives :
$\bbox[5px,border:1px solid] {(\frac{du}{d\theta})^2+u^2=\frac{\mu^2e^2}{C^4}}\quad$ where $\ \bbox[5px,border:1px solid]{e=\sqrt{1+\frac{2hC^2}{\mu^2}}}\ $ is called the eccentricity.
This integrates into $u=\frac{\mu e}{C^2}\cos(\theta-\theta_0)$ and we will introduce yet some new notations.
$\bbox[5px,border:1px solid] {\begin{array}{c} \nu=\theta-\theta_0 \\ p=\frac{C^2}{\mu} \\ \end{array}}\quad$ giving $\quad\bbox[5px,border:3px solid blue] {r=\frac{p}{1+e\cos \nu}}$
We can remark that :


*

*$e=0\quad$ this is a circle

*$h<0\Rightarrow 0<e<1\quad$ this is an ellipse

*$h=0\Rightarrow e=1\quad$ this is a parabola

*$h>0\Rightarrow e>1\quad$ this is an hyperbola
The case which interest us currently is the ellipse. Let's call $2a$ the length of the major axis. 
We see that $r$ is minimum for $\nu=0$ and $r_{min}=a(1-e)=\frac p{1+e}$
$\bbox[5px,border:1px solid] {p=a(1-e^2)}\quad$ remark that we have also $a=\frac{\mu}{-2h}$, we will use it later.
Note: for Earth $e\simeq0.0167$ and $a\sim1\,UA\simeq149\,598\,000\ km$

Third law of Kepler
Let's introduce our new angles.


*

*$\nu$ is called the true anomaly 

*$E$ is called the eccentric anomaly

*$M=n(t-t0)=2\pi(\frac{t-t0}{T})$ is called the mean anomaly, $n$ the frequency and $T$ the period.
With the new angle $E$ we have
$\begin{cases}
x=r\cos\nu=a(\cos E-e) \\
y=r\sin\nu=b\sin E=a\sqrt{1-e^2}\sin E \\
\end{cases}$
The equation of the ellipse becomes $\bbox[5px,border:3px solid blue] {r=a(1-e\cos E)}$
We deduce also without much trouble that $\bbox[5px,border:3px solid green] {\tan(\frac{\nu}{2})=\sqrt\frac{1+e}{1-e}\tan(\frac E2)}$
The tedious part is establishing the Kepler's equation. I will not detail everything here, but basically you define $dt=\frac{r}{\sqrt{-2h}}d\tilde E$ an elementary area element.
Then by replacing $r,\dot r,\ddot r$ by the appropriate terms $r,r'=\frac{dr}{d\tilde E},r''=\frac{d^2r}{d\tilde E^2}$ in the equation $r\ddot r+{\dot r}^2=2h+\frac{\mu}{r}$ obtained by eliminating $\dot \theta$ from the movement equation.
You eventually arrive to $\bbox[5px,border:1px solid] {r''+r=a}$ which combined to $\sqrt{-2h}(t-t_0)=\int\limits_0^{\tilde E} r(\tilde E)d\tilde E$ proves to be compatible with the $E=\tilde E$ and $r(E)$.
The latter integral relation gives you $\sqrt{-2h}(t-t_0)=a(E-e\sin E)$
We define $\bbox[5px,border:3px solid green] {M=2\pi(\frac{t-t_0}{T})=E-e\sin E}$ and gain $Ma=\frac{2\pi a}{T}=\sqrt{-2h}=\sqrt\frac{\mu}{a}$ 
which is the famous $\bbox[5px,border:1px solid] {\frac{T^2}{a^3}=\frac{4\pi^2}{\mu}}\ $.

As promised we have now the two blue movement trajectories $r(\theta),r(E)$ expressed in two different frames, and the angles $\theta,E$ are linked to the time $t$ thanks to the green equations.
A: Without the substitution $r=1/u$ it will indeed be a mess. $ r, \theta $ as time dependents are involved with elliptic integrals/functions.
The equations of motion can be numerically integrated more easily.
Unfortunately the planetary orbit equation in its simplest form ( major axis along x-axis, focus at origin) combines  coordinates from Cartesian and Polar forms
$$ p =r - e\cdot x $$
which when differentiated with respect to arc I found leads to a simple form as its geometric property:
$$ \dfrac{\cos \psi}{\cos \phi} =e. $$
( Angles which conic tangent makes to radius and x-axis respy.)
