Solving $\textbf{r}''(t)=\frac{GM}{(r(t))^3}\textbf{r}(t)$ I have the following problem that I don't know how to solve:

Let $\textbf{r}(t)$ be a curve and let $r(t):= |\textbf{r}(t)|$. Let the curve $\textbf{r}(t)$ be defined by:
$$\textbf{a}(t)=\frac{GM}{(r(t))^3}\textbf{r}(t)$$
Where $\textbf{a}(t) = \textbf{r}''(t)$ is the acceleration of the particle and $G,M \in \mathbb{R}$. Suppose that $\textbf{r}(t)$ and $\textbf{r}'(t)$ are not parallel.
The objective is to determine the trajectory $\textbf{r}(t)$ in polar coordinates.

So basically we have the following equation:
$$\textbf{r}''(t)=\frac{GM}{(r(t))^3}\textbf{r}(t)$$
The thing is that I don't know how I should proceed from now on. If we'd have just normal scalar function this would be a differential equation and I think I would be able to solve it, but like this I have no idea how to approach the problem. How can I solve this?

Edit:
If $$\textbf{r}(t):=(x(t),y(t))$$ then the equation turns into the following system of equations:
$$x''=\frac{GM}{(\sqrt{x^2 + y^2})^3}x$$
$$y''=\frac{GM}{(\sqrt{x^2 + y^2})^3}y$$
I don't know if this helps but I think that this turns thing easier perhaps.
 A: Your equation is wrong (you forgot the minus sign). It should instead be
$$\ddot{\mathbf{r}}=\frac{-GM}{r^3}\mathbf{r}$$
Or,
$$\ddot{\mathbf{r}}=\frac{-GM}{r^2}\hat{\mathbf{r}} ~~~~~(1)$$
With that out of the way, let's get started. It is well known that in polar coordinates,
$$\ddot{\mathbf{r}} =\left(\ddot{r} -r\dot{\theta }^{2}\right)\hat{\mathbf{r}} +( r\ddot{\theta } +2\dot{r}\dot{\theta })\widehat{\mathbf{e}_{\theta }}$$
However, since the gravitational force is entirely radial, in our case this is simply
$$\ddot{\mathbf{r}} =\left(\ddot{r} -r\dot{\theta }^{2}\right)\hat{\mathbf{r}}$$
Then, returning to (1),
$$\ddot{r} -r\dot{\theta }^{2}=\frac{-GM}{r^2}$$
Now here's the clever part: We'll make a change of variable $u \equiv \frac{1}{r}$. Therefore
$$\dot{r}=\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\frac{\mathrm{d}\theta}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}$$
We now recall the fact that $ r\ddot{\theta } +2\dot{r}\dot{\theta }=0$ and also observe that
$$r\ddot{\theta } +2\dot{r}\dot{\theta }=\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})=0$$
Therefore $r^2\dot{\theta}$ is a conserved quantity, and we make the definition $h \equiv r^2\dot{\theta}=\frac{1}{u^2}\dot{\theta}$. Therefore,
$$\dot{r}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\frac{h}{r^2}=-h\frac{\mathrm{d}u}{\mathrm{d}\theta}$$
So,
$$\ddot{r}=\frac{\mathrm{d}}{\mathrm{d}t}\left(-h\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)=-h^2u^2\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}$$
Then returning to (1) again, we can substitute in our new variables to obtain
$$-h^2u^2\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}-\frac{1}{u}\left(hu^2\right)^2=-GMu^2$$
$$\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}+u=\frac{GM}{h^2}$$
This is a second order linear constant coefficient ODE for $u(\theta)$ which has the following solution:
$$u(\theta)=A\sin(\theta)+B\cos(\theta)+\frac{GM}{h^2}=\frac{1}{r(\theta)}$$
And if we let $r(\theta)=r_{\text{max}}$ at $\theta=0$ and $r'(\theta)=0$ at $\theta=0$, then we find that $A=0$ and $B=\frac{1}{r_{\text{max}}}-\frac{GM}{h^2}$
Thus
$$r(\theta)=\frac{1}{B\cos(\theta)+\frac{GM}{h^2}}$$
If we define $C \equiv \frac{h^2B}{GM} = \frac{h^2}{GMr_{\text{max}}}-1$,
$$r(\theta)=\frac{h^2/GM}{1+C\cos(\theta)}$$
More work is needed to express this in terms of the semi minor axis $a$, the semi major axis $b$, and the eccentricity $\epsilon$, but I hope this is satisfactory enough. Finding  $r(\theta)$ instead of $r(t)$ is considerably more useful, so I would recommend to stick to this sort of formulation. Please ask if anything is unclear!

THE PARAMETRIC FORM (if you're interested)
The parametric form of the ellipse is $$\mathbf{r}(t)=a\cos(t)\hat{\mathbf{i}}+b\sin(t)\hat{\mathbf{j}}$$
Which in polar coordinates is
$$\mathbf{r}(t)=\sqrt{a^2\cos^2(t)+b^2\sin^2(t)}\hat{\mathbf{r}}+\operatorname{atan2}\left(b\sin(t),a\cos(t)\right)\widehat{\mathbf{e}_{\theta}}$$
Information on $\operatorname{atan2}$ can be found here
EDIT: A more general analysis:
Suppose our equation is
$$\ddot{\mathbf{r}}=\frac{D}{r^3}\mathbf{r}$$
Where $D \in \mathbb{R}$ is some positive or negative constant. Then by similar analysis as before we reach
$$r(\theta)=\frac{1}{B\cos(\theta)-\frac{D}{h^2}}$$
The solutions behave rather differently based on the sign of $D$. If it is negative, we get elliptical orbits as before. If $D$ is positive, however, the "orbits" (bit of an abuse of terminology, as they aren't periodic) are either parabolic or hyperbolic, depending on $h$, $r(0)$, and $r'(0)$.
