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\;.$
Let's start from
$$
L = \frac{1}{2}(\dot{r}^2 + r^2\dot{\theta}^2) - V(r)
$$
There are two degrees of freedom: $r$ and $\theta$. Start with $\theta$
$\theta$
$$
\frac{{\rm d}}{{\rm d}t}\frac{\partial L}{\partial \dot{\theta}} = 0 = \frac{{\rm d}}{{\rm d}t}(r^2\dot{\theta})~~~\Rightarrow~~~ r^2\dot{\theta} = l = {\rm const} \tag{1}
$$
which means that the number $l = r^2\dot{\theta}$ is a constant (angular momentum!).
$r$
$$
\frac{{\rm d}}{{\rm d}t}\frac{\partial L}{\partial \dot{r}} - \frac{\partial L}{\partial r} = 0 = \ddot{r} - r\dot{\theta}^2 + \frac{{\rm d}V}{{\rm d}r} \tag{2}
$$
Use (1) to write
$$
\frac{{\rm d}}{{\rm d}t} = \frac{l}{r^2}\frac{{\rm d}}{{\rm d}\theta}
$$
and define the variable $u = 1/r$, if you replace both into (2) you will get
$$
\frac{{\rm d}^2u}{{\rm d}t^2} + u = \frac{1}{l^2 u^2}\frac{{\rm d}V(1/u)}{{\rm d}r} \tag{3}
$$
For a Keppler potential
$$
V(r) = -\frac{GM}{r} = -GMu
$$
Replace that in (3) and you get
$$
\frac{{\rm d}^2u}{{\rm d}t^2} + u = \frac{GM}{l^2}
$$
whose solution is
$$
u(\theta) = C\cos(\theta - \theta_0) + \frac{GM}{l^2}
$$
Now define
$$
e = \frac{Cl^2}{GM} ~~~\mbox{and}~~~ a = \frac{l^2}{GM(1 - e^2)}
$$
such that the equation above becomes
$$ \bbox[5px,border:2px solid blue]
{
r(\theta) = \frac{a(1 - e^2)}{1 + e\cos(\theta - \theta_0)}
}
\tag{4}
$$
And there you have it, the solution are conics
Best Answer
Absorb the dimensional GM into the units, to remove excuses for not recognizing the plane-geometry structure. Note the rotational and translational invariance to be used in fixing your coordinate system below. All vectors are then 2-vectors on that plane, $$\ddot{ x }=-{ { x}\over r^3},\qquad \ddot{ y }=-{ { y}\over r^3} ~. $$ It is then self-evident that $$ L=x\dot y -y\dot x $$ is a constant of the motion (in the z-direction, the only non vanishing one, of course), $\dot L=0$.
Moreover, the mass/normalized LRL 2-vector on that plane, $$ \vec e= L \begin{pmatrix} \dot y\\-\dot x\end{pmatrix}-{1\over r} \begin{pmatrix} x\\ y\end{pmatrix} $$ is also conserved, $\dot{\vec e }=0$. Its magnitude will turn out to be the eccentricity.
Dotting by $\vec r$ you have $$ \vec r\cdot \vec e =L^2 -r, \leadsto \\ r= L^2-\vec r\cdot \vec e , \leadsto \\ x^2+y^2= (L^2-\vec r\cdot \vec e )^2. $$ You may use rotational invariance to take $\vec e$ along the x-axis, $\vec e=-\epsilon \hat x$ to prettify your ellipse orientation, and work out the r.h.s. square, as a quadratic polynomial in x, with the obvious constants. Elementary algebra leads you to the ellipse your teacher taught you in terms of the constants ε and L. That is, shifting the origin of xs to $L^2 \epsilon/(1-\epsilon^2)$ and taking $\epsilon ^2= 1-b^2/a^2$, you get $$ \frac{x^2}{a^2}+\frac{y^2}{b^2} = L^4 a^2/b^4. $$