This was previously a comment to space_cadet's answer but became long (down-vote wasn't me though).
I don't understand space_cadet's talk about unstable orbits. Recall that two-body system with Coulomb interaction has an additional $SO(3)$ symmetry and has a conserved Laplace-Runge-Lenz vector which preserves the eccentricity. Because interactions between planets themselves are pretty negligible one needs to look for explanation elsewhere. Namely, in the initial conditions of the Solar system.
One can imagine slowly rotating big ball of dust. This would collapse to the Sun in the center a disk (because of preservation of angular momentum) with circular orbits and proto-planets would form, collecting the dust on their orbits. Initially those planets were quite close and there were interesting scattering processes happening. The last part of the puzzle is mystery though. If there were still large amount of dust present in the Solar system it would damp the orbits to the point of becoming more circular than they are today. The most popular explanation seems to be that the damping of the eccentricity was mediated by smaller bodies (like asteroids). Read more in "Final Stages of Planet Formation" - Peter Goldreich, Yoram Lithwick, Re'em Sari.
Short answer: Use the radial acceleration component $a_\mathrm r$ of the total $\bf a$ by the central force and then deduce the SHM equation using the properties of circular orbital motion (effective potential energy curves can also be used) that $\dot r= 0$ and keeping in mind the fact that SHM would only occur when the circular orbit is stable.
Here the approach is elaborately discussed:
The radial component of the central force $\bf F$ is given by
$$\mathbf F\cdot \mathbf e_\mathrm r~=~ m\left[\ddot r - r(\dot \theta)^2\right]\;.\tag 1 $$
Due to the small perturbation, the radius of the circular orbit changed from $r_0$ to $r(t) ~=~r_0 + \rho(t)~~_; ~~~\rho(0)\ll r_0\;.$
Using $\mathit l= mr^2\dot \theta=\textrm{const}$, $(1)$ becomes
$$\ddot\rho(t)- \frac{l^2}{m^2(r_0+ \rho(t))^3}~=~ \frac{F_\mathrm r(r_0+ \rho(t))}{m}~\tag {1-a} $$
Now, \begin{align}(r_0+ \rho(t))^{-3}&=r_0^{-3}\left(1+ \frac{\rho}{r_0}\right)^{-3}\\ &\approx r_0^{-3}\left(1-3\frac{\rho}{r_0}\right)\\ F_\mathrm r(r_0+\rho)&= F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho+ \ldots \\ &\approx F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho\end{align}
So,
$$\ddot\rho- \frac{l^2}{m^2r_0^3}\left(1-3\frac{\rho}{r_0}\right)= \frac{F_\mathrm r(r_0)}m +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\cdot \frac{\rho}m\tag{1-b} $$
From the orbit-equation$^\ddagger$, we get
$$\mathit l^2=- mr_0^3 F_\mathrm r(r_0)_; $$
So,
$$\ddot \rho+ \left[\frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\right]\rho~=~0\;.$$
This can be written as
$$\ddot \rho + \omega^2 \rho ~=~0\;.\tag 2$$
Now, this is a second-order differential equation for simple harmonic oscillator with frequency $\omega\;,$
where $$\omega^2 ~=~ \frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\;.$$
Unless, the circular orbit is unstable, a simple harmonic radial oscillation about $r=r_0$ would ensue (i.e. SHM would occur if $\omega^2\gt 0\;.$)
$\ddagger$ We would derive the orbit-equation from $(1)$ expressing $r$ and its derivatives in terms of $\theta\;.$
\begin{align}\dot r &= \frac{\mathrm dr}{\mathrm d\theta}\frac{\mathrm d\theta}{\mathrm dt}\\ \ddot r&=\frac{\mathrm d}{\mathrm dt}\left(\frac{\mathrm dr}{\mathrm d\theta}\dot \theta\right)\\&=\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta \end{align}
Let
$$\mathbf H \stackrel{\text{def}}{=}
\mathbf r\times \dot{\mathbf r} =\frac{\mathbf L}m \;.$$
Therefore
\begin{align}\dot \theta &= \frac{|\mathbf H|}{r^2}\\ \implies \ddot \theta &=-\frac{2|\mathbf H|}{r^3} \, \dot r\;.\end{align}
Putting all these in $(1),$ we get
\begin{align}\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta- r(\dot \theta)^2& = \frac{F_\mathrm r}{m} \\ \implies\frac{\mathrm d^2r}{\mathrm d\theta^2}~\left(\frac{|\mathbf H|}{r^2}\right)^2 +\frac{\mathrm dr}{\mathrm d\theta}~\left(-\frac{2|\mathbf H|}{r^3}\dot r\right) - r ~\left(\frac{|\mathbf H|}{r^2}\right)^2 &= \frac{F_\mathrm r}{m}\\ \implies \left(\frac{|\mathbf H|}{r^2}\right)^2\left[\frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r\right]&=\frac{\mathrm F(r)}{m}\\ \implies \frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r&= \frac{r^4 F_\mathrm r}{|\mathbf H|^2 m}\;.\tag{i}\end{align}
$(\rm i)$ is the orbital equation.
Now, applying the facts that for circular orbits, $\dot r,~\ddot r~=~0\,,$ we get the value of $$|\mathbf H|~=~\sqrt{ -\frac{r^3F_\mathrm r}{m}}\;,$$ where $r~=~r_0\;.$
Best Answer
Not every central force admits a circular orbit since only an attractive interaction can balance the repulsive centrifugal term. On the other hand, every attractive force has a circular orbit since by an appropriate choice of angular momentum the centrifugal term, $L^2/2\mu r^2$, can be chosen to cancel the attractive interaction, $U(r)$. Finally, attractive potential is not a sufficient for the existence of stable orbits.
The circular orbit of radius $r_0$ is stable if and only if $U_{\mathrm{eff}}(r_0)$ corresponds to a minimum of the effective potential, $$U_{\mathrm{eff}}(r)=\frac{L^2}{2\mu r^2}+U(r).$$ This implies that the second derivative of $U_{\mathrm{eff}}$ at $r_0$ must be positive. Hence, $$U''(r_0)>-\frac{3L^2}{\mu r_0^4}.\tag1$$
For circular orbits, the radial effective force (which includes interaction and centrifugal forces) vanishes and then $U_{\mathrm{eff}}(r_0)'=0$. Thus $$r_0^3=\frac{L^2}{\mu U'(r_0)}.$$ Plugging this into Eqn (1) we obtain the following condition $$\frac{U''(r_0)}{U'(r_0)}+\frac{3}{r_0}>0,$$ for stable orbits. As an example, assuming an attractive potential with a power law $U=kr^{n}$, the last equation gives us that is has stable orbits only for $n>-2$.
It is worth mentioning that there are different concepts of stability regarding orbital motion. The one assumed here, for which the above result holds, says that a circular orbit is stable if it remains bounded under small perturbations (a bounded orbit is one whose radius is limited by $r_{\mathrm{min}}\leq r\leq r_{\mathrm{max}}$). A different and also common concept is Lyapunov Stability. In this case, among all power law central forces, only the harmonic oscillator gives stable orbits.