How to solve the inverse square law equation of motion From $$m\boldsymbol{\ddot{r}}=\boldsymbol{\hat{r}} f(r)$$ I can get
$$r''-r \theta '^2=-\frac{k}{m r^2}$$
$$2 r' \theta '+r \theta ''=0$$
Now it seems that all the books tells me the method to solve this equation by eliminating $t$ and obtain the equation of the orbit. But I am wondering if it is possible to solve the two equation of $r$ with respect to $t$ and $\theta$ with respect to $t$.
 A: 
But I am wondering if it is possible to solve the two equation of r with respect to t and θ with respect to t.

Yes, you can. In the case of an elliptical orbit with non-zero angular momentum, you need to introduce the concepts of eccentric anomaly, mean anomaly, and mean motion.
Eccentric anomaly is related to true anomaly (your theta) via
$$\sqrt{1-\varepsilon} \tan \frac \theta 2 = \sqrt{1+\varepsilon} \tan \frac E 2 $$
where $\varepsilon$ is the eccentricity of the orbit.
Eccentric anomaly is related to mean anomaly via Kepler's equation,
$$M = E - \varepsilon \sin E$$
Mean anomaly is a linear function of time, with the mean mean motion $n$ specifying the time constant:
$$M(t) = M(t_0) + n(t-t_0)$$
Finally, the mean mean motion $n$ is a constant determined by the masses of the orbiting bodies and the distance at which they are orbiting one another:
$$n = \sqrt {\frac {\mu_1+\mu_2}{a^3}}$$
where $a$ is the semi-major axis of the orbit, and $\mu_1$ and $\mu_2$ are the gravitational parameters of the orbiting bodies: $\mu_k = GM_k$, where $G$ is Newton's gravitational constant and $M_k$ is the mass of body $k$. (But in practice, one typically uses the gravitational parameter rather than mass and $G$. The gravitational parameter is highly observable. Mass and the gravitational constant are not.)
In the case of a parabolic or hyperbolic trajectory, you will need to use the parabolic or hyperbolic anomaly in lieu of the elliptical anomaly. In the case of zero angular momentum, you will need to solve the radial Kepler problem.

To find $\theta(t)$, you'll need the mean motion $n$, the eccentricity $\varepsilon$, and the value of $\theta$ at some epoch time $t_0$, $\theta_0$. The above equations let's you compute the initial mean anomaly $M_0$ given this initial true anomaly $\theta_0$ and the eccentricity. The mean anomaly at some time $t$ is trivial as mean anomaly is a linear function of time. The tricky part is finding the eccentric anomaly at that time.
Kepler's equation is a transcendental equation. It certainly does have an inverse, but that inverse cannot be expressed in terms of the elementary functions. Many hundreds of papers have been written ways to solve the inverse Kepler problem in the four centuries since Kepler found that planetary orbits are ellipses. The most widely used approach is to use a Newton-Raphson iterator. This works great so long as the eccentricity isn't huge ($\varepsilon < 0.8$). Specialized techniques are needed for highly eccentric orbits.

Alternatively, you can entirely forego the Keplerian approach and integrate the equations of motion numerically.  Keplerian orbits are a nice fiction. They only exist in the case of a two body problem that follows Newton's law of gravitation. Our solar system has a central mass (the Sun), eight planets (Mercury, Venus, ... Neptune), and a pile of lesser bodies (Pluto, Ceres, etc.) Orbits in our solar system are only approximately elliptical. Moreover, Newton's law of gravitation is only approximately correct. General relativity more accurately describes orbital motion than does Newtonian mechanics.
There is no way to solve this simply, or even iteratively. The only way around this mess is to use numerical integration techniques. Now we're about not just hundreds of technical papers but many thousands of papers, some of which are being published to this very day.
A: It is possible to find $t(r)$ and $t(\theta)$ easily, however inverting to find $r(t)$ and $\theta(t)$ is hard to do in general (unless you use Fourier series etc..)
\begin{align*}
t(r) &= \sqrt{\frac{m}{2}} \int \frac{dr}{\sqrt{\frac{k}{r} - \frac{l^2}{2mr^2} + E}} \\
t(\theta) &= \frac{l^3}{mk^2} \int \frac{d\theta}{\left( 1 + e \text{cos }(\theta - \theta')\right)^2}
\end{align*}
You can solve these integrals in terms of elementary functions, but once you do you will see that you cannot invert.
As @David Hammen noted, you can define some auxiliary variables to simplify your expressions, but at the end the relation to invert will remain trancendental and you can only have $r(t)$ implicitly.
