The usual Euler-Lagrange equations can be derived by considering a variation of the action $$S(\mathcal{L})=\int\mathcal{L}(t,q,\dot{q})dt\tag{1}$$ taking a variation $q\mapsto q+\delta q$ and expanding the legrangian to get $$S(\mathcal{L}') = \int\mathcal{L}+\frac{\partial\mathcal{L}}{\partial q}\delta q + \frac{\partial\mathcal{L}}{\partial \dot{q}}\delta \dot{q} + O(\delta q^2) \ dt.\tag{2} $$ The change of the action is then given by $$\Delta S = S(\mathcal{L}')-S(\mathcal{L})=\int\frac{\partial\mathcal{L}}{\partial q}\delta q + \frac{\partial\mathcal{L}}{\partial \dot{q}}\delta \dot{q}dt=\int\left(\frac{\partial\mathcal{L}}{\partial q} + \frac{d}{dt}\frac{\partial\mathcal{L}}{\partial \dot{q}}\right)\delta q \ dt\tag{3}$$ where the last equality follows by integration by parts, and the fact that the variation fixes endpoints ($\delta q=0$ on the boundary). The condition for an extrema is that $\Delta S = 0$ for any choice of $\delta q$, which gives us the Euler-Lagrange equations.
My problem is that I am trying to adapt this type of reasoning to figure out a similar formula for the second variation. My guess for how to do this is to consider the second order terms in $\delta q$ which we ignored above. Working this out gave $$\frac{1}{2}\int\frac{\partial^2\mathcal{L}}{\partial q^2}\delta q^2+2\frac{\partial^2\mathcal{L}}{\partial q \partial \dot{q}}\delta q\delta\dot{q}+\frac{\partial^2\mathcal{L}}{\partial \dot{q}^2}\delta\dot{q}^2 \ dt.\tag{4}$$ From here, my guess is that we want to find some total derivatives which allow us to convert the second two terms into ones whos coefficient is $\delta q^2$. For example, consider $$\frac{d}{dt}\left(\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\delta q^2\right) = \frac{d}{dt}\left(\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\right)\delta q^2 + 2\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\delta q\delta\dot{q}.\tag{5}$$ Integrating, this becomes $$\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\delta q^2\bigg|_{\text{boundary}} = \int\frac{d}{dt}\left(\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\right)\delta q^2 + 2\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\delta q\delta\dot{q}dt\tag{6}$$ or, since $\delta q=0$ on the boundary, $$\int 2\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\delta q\delta\dot{q}dt=-\int\frac{d}{dt}\left(\frac{\partial^2\mathcal{L}}{\partial q\partial\dot{q}}\right)\delta q^2 dt\tag{7}$$ Which gives the desired form of the middle term.
My issue is with the third term. I guessed $$\frac{d}{dt}\left(\frac{\partial^2\mathcal{L}}{\partial\dot{q}^2}\delta q^2-2\frac{\partial^2\mathcal{L}}{\partial\dot{q}^2}\delta q\delta\dot{q}\right)=\frac{d^2}{dt^2}\left(\frac{\partial^2\mathcal{L}}{\partial\dot{q}^2}\right)\delta q^2 - 2\frac{\partial^2\mathcal{L}}{\partial\dot{q}^2}\delta\dot{q}^2-2\frac{\partial^2\mathcal{L}}{\partial\dot{q}^2}\delta q\delta\ddot{q}\tag{8}$$ which gets a little close, but I can't seem to think of how to get rid of the $\delta\ddot{q}$ term.
Have I missed something simple, or is it not possible to get it into the desired form?
I'm too lazy to type out the full solution, but let me give some hints:
If the Lagrangian $L(q,t)$ does not depend on $\dot{q}$, then the 2nd functional derivative of the action $S=\int \!dt~L$ is $$ \frac{\delta^2 S}{\delta q(t)\delta q(t^{\prime})}~=~\frac{\partial^2 L(t)}{\partial q(t)\partial q(t)}\delta(t\!-\!t^{\prime}).$$ It turns out there will be 3 more terms if the Lagrangian $L(q,\dot{q},t)$ depends on $\dot{q}$. It will still be true that $\frac{\delta^2 S}{\delta q(t)\delta q(t^{\prime})}$ has support on the diagonal $t=t^{\prime}$.
So integral $\int \!dt$ in OP's eq. (4) should be replaced by a double-integral $\iint \!dt~dt^{\prime}~\delta(t\!-\!t^{\prime})$.
As usual integrate by parts either wrt. $\frac{d}{dt}$ or $\frac{d}{dt^{\prime}}$ to get rid of derivative(s) of $\delta q$.
For further clues, see my related Phys.SE answer here.