The simplest way to keep track of all this is to introduce a "dummy" parameter $\epsilon$ which just "counts how small things are".  Thus start with your initial set of full equations
\begin{align*}
(l+a)\ddot\theta + 2\dot\theta\dot a + g\sin\theta &= 0 \tag{1} \\
\ddot a - (l+a)\dot\theta^2 + g(1-\cos\theta) + \frac{k}{M}a &= 0. \tag{2}
\end{align*}
and, under the approximations where $\theta$ and $a$ are small, replace
$\theta\to \epsilon\theta$ and $a\to \epsilon a$.  With this
$\dot \theta \to \epsilon \dot\theta$ etc so your equations become
\begin{align*}
(l+\epsilon a)\epsilon \ddot\theta + 
2\epsilon^2\dot\theta\dot a + g\sin\epsilon\theta &= 0 \tag{3}\\
\epsilon\ddot a - (l+\epsilon a)\epsilon^2 \dot\theta^2 + g(1-\cos\epsilon\theta) + \epsilon\frac{k}{M}a &= 0. \tag{4}
\end{align*}
You can then *linearize* your equations of motion, i.e. expand everything to terms *linear* in $\epsilon$, meaning you dump anything $\epsilon^2$ or above.
This immediately gives
\begin{align}
\epsilon l\ddot{\theta}+g\epsilon\theta&=0\\
\epsilon \ddot{a}+\epsilon\frac{k}{M}a&=0\, .
\end{align}
If you have done the job right, the counter $\epsilon$ just drops out (as it does here).  This makes it clear that terms in your first equation like $a\ddot{\theta}$ and $\dot{\theta}\dot{a}$ are of size $\epsilon^2$ and can be ignored.  Likewise this cleanly kills in your second equation the entire $(l+\epsilon a)\epsilon^2\dot{\theta}^2$ as it contains terms in $\epsilon^2$ and $\epsilon^3$, and also kills the $(1-\cos\epsilon \theta)$ term.

In problems such as the one you have, nothing replaces being systematic and carefully keeping track of what's expected to be small.

-------

Edit 1: resonance requires an external driving force and the response of the system to resonant or near resonant input is not linear: the amplitude of the oscillation is dominated by a non-linear factor of $1/\sqrt{(\omega^2-\omega_0^2)^2+(\omega_0^2\omega^2/Q^2)}$ which is nearly singular when the quality factor $Q$ is small.  As a result the assumption of small amplitude of oscillation near a minimum of the effective potential is not valid and linearizing will produce nonsense (terrible approximate solution).

-----------
Edit 2: If you want "the next order", you need to solve by successive approximation.  I will sketch this only for the $\theta$ degree of freedom.  You assume
\begin{align}
\theta(t)&=A\cos\omega_\theta t + \epsilon F_\theta(t) \tag{5}\\
a(t)&=B\cos\omega_a t +\epsilon G_a(t) \tag{6}
\end{align}
(I trust the notation is obvious.)  Inserting Eq.(5) and Eq.(6) into the equation of motion for $\theta$, Eq.(3), and setting separate powers of $\epsilon$ to $0$; if my algebra is right, you get
\begin{align}
0&=\epsilon\,A \cos(\omega_\theta t)
\left(g - l\omega^2_\theta\right) \tag{7} \\
0&=\epsilon^2\left(-A B \omega^2_\theta \cos(\omega_a t)\cos(\omega_\theta t)+g F_\theta(t)+2 A B \omega_a\omega_\theta \sin(\omega_\theta t)\sin(\omega_a t)+ l \ddot{F_\theta}\right)\, . \tag{8}
\end{align}

Amazingly (and still provided my algebra is right!), this equation is still decoupled in $\epsilon^2$, in the sense it involves only $F_\theta$ but does not involve $G_a$ so you are left with a standard 2nd order differential equation equivalent to a driven oscillator.

The terms in $\cos(\omega_a t)\cos(\omega_\theta t)$ and $\sin(\omega_\theta t)\sin(\omega_a t)$ indicate that you must use
as ansatz 
$$
F_\theta(t)=C \cos\left((\omega_\theta-\omega_a)t\right)+ D \cos\left((\omega_\theta+\omega_a)t\right)\, ,
$$
(there might also be terms in 
$\sin\left((\omega_\theta\pm\omega_a)t\right)$ also depending on your initial conditions). 

You need to watch carefully for some resonance conditions, which would invalidate the assumption that the extra term $\epsilon F_\theta$ is small.  In general, one must also watch for the appearance of so-called *secular* terms. I don't think this happens here, but when this occurs, one must then write the first order frequencies $\omega_a$ and $\omega_\theta$ as a series of the form
$$
\Omega_\theta^2=\omega_\theta^2+\epsilon (\omega^{(1)}_\theta)^2+\ldots
$$
in a procedure called the Lindstedt-Poincare method.