In your answer you wrote $\int \limits_{t_1}^{t_2} \frac{\partial L}{\partial q} \delta q dt = [ \frac{\partial L}{\partial q} \delta q(t) ]_{t_1}^{t_2}$ and even said this was zero which is a deep misunderstanding that should carefully be explained.
Lets take the Harmonic oscillator Lagrangian:
\begin{align}
S(q) = \int_{t_1}^{t_2} (\frac{1}{2} m \dot{q}^2 - \frac{1}{2} k q^2)dt
\end{align}
The Euler-Lagrange equations are
\begin{align}
0 &= \frac{d}{dt} \frac{\partial L}{\partial \dot{q}} - \frac{\partial L}{\partial q} \\
&= \frac{d}{dt} [\frac{\partial }{\partial \dot{q}}(\frac{1}{2} m \dot{q}^2 - \frac{1}{2} k q^2)] - \frac{\partial }{\partial q}(\frac{1}{2} m \dot{q}^2 - \frac{1}{2} k q^2) \\
&= \frac{d}{dt} [m \dot{q}] + kq \\
&= m \ddot{q} + k q
\end{align}
which is the same result as Newton's law for the Harmonic oscillator. This is the result we want to end up with from first principles by following the same steps we use to derive the Euler-Lagrange equations, but for the specific Lagrangian above rather than an abstract one.
Recall that in general $S(q)$ is just a function of functions $q$ (which are a function of $t$) and in general we can plug any random $q$ into $S$, so lets take an example, say the $q$ for which $q(t) = t^3$.
\begin{align}
S(q) &= \int (\frac{1}{2} m [\frac{d}{dt} t^3]^2 - \frac{1}{2} k [t^3]^2)dt = \int (\frac{1}{2} m [3 t^2]^2 - \frac{1}{2} k t^6)dt = \int_{t_1}^{t_2} (\frac{9}{2} m t^4 - \frac{1}{2} k t^6)dt \\
&= (\frac{9}{2} m \frac{t^5}{5} - \frac{1}{2} k \frac{t^7}{7})|^{t_2}_{t_1}.
\end{align}
This is now just some number, in other words we just did an integral of some function of $t$, for the special case of $q(t) = t^3$, but I could have chosen a different $q(t)$, e.g. $q(t) = 2t^4$, and gotten a different number as my answer. The point is I am integrating a function of functions and I cannot actually perform the integral until I specify what my function is.
Recall
\begin{align}
\delta q &= (q + \delta q) - q, \\
\delta q^2 &= (q + \delta q)^2 - q^2 = q^2 + 2 q \delta q + (\delta q)^2 - q^2 = 2 q \delta q, \\
\delta \dot{q}^2 &= (\dot{q} + \delta \dot{q})^2 - \dot{q}^2 = \dot{q}^2 + 2 \dot{q} \delta \dot{q} + (\delta \dot{q})^2 - \dot{q}^2 = 2 \dot{q} \delta \dot{q},
\end{align}
where I keep only terms to first order, and notice in $\delta \dot{q}^2$ we just treat $\dot{q}$ as a variable in and of itself, ignoring the dot (time derivative). But what does $q + \delta q$ mean? $q + \delta q$ is just a new function, a function of $t$, $(q + \delta q)(t)$, and so for $q(t) = t^3$ it means I am shifting the function $t^3$ into some new function $(q + \delta q)(t)$, which is potentially any function in general, but it's supposed to mean a function 'close' to $t^3$. For example,
\begin{align}
(q + \delta q)(t) &= (t + \varepsilon)^3,
\end{align}
is close to $q(t) = t^3$, where $\varepsilon$ is some small real number,
\begin{align}
q(t) &= t^3, \\
(q + \delta q)(t) &= (t + \varepsilon)^3, \\
\delta q(t) &= [q + \delta q](t) - q(t) = (t + \varepsilon)^3 - t^3 = 3t^2 \varepsilon + 3 t \varepsilon^2 + \varepsilon^3 = 3t^2 \varepsilon
\end{align}
where in the last equality I keep terms only to first order in $\varepsilon$ because my ultimate goal is to minimize $S(q)$ which involves the first order change in $S$ caused by shifting $q(t) = t^3$ into $[q + \delta q](t) = (t + \varepsilon)^3$ i.e. first order contributions of $\varepsilon$ and not higher order contributions like $\varepsilon^2$. In other words, in a Taylor expansion of $S(q+\delta q) = S[(t + \varepsilon)^3]$ with $\varepsilon$ I expect the first order terms to be like first derivatives of $S$ as a function of $\varepsilon$ which by basic calculus lets us minimize the single-variable function $\tilde{S}(\varepsilon) = S[(t + \varepsilon)^3]$ using basic calculus.
Since we are always treating $(q + \delta q)(t)$ as a first order correction to $q(t)$ we can say
\begin{align}
(q + \delta q)(t) = q(t) + \delta q(t)
\end{align}
which for our example reads as
\begin{align}
(q + \delta q)(t) &= (t + \varepsilon)^3 \\
&= t^3 + 3 t^2 \varepsilon + 3 t \varepsilon^2 + \varepsilon^3 \\
&= t^3 + 3 t^2 \varepsilon \\
&= q(t) + \delta q(t).
\end{align}
In general, how do you determine whether one function is close to another? We need to use metric/normed space theory and define a distance function on our space of functions and then $q(t)$ is any function for which this distance we define carefully is small. This is something studied carefully in the 'calculus of variations', but in physics one ignores this, so lets use $(t + \varepsilon)^3$ as our small shift away from $t^3$. In other words, for $q(t)$ lets always set
\begin{align}
\delta q(t) = q(t + \varepsilon) - q(t).
\end{align}
Notice also that
\begin{align}
\dot{q}(t) &= 3 t^2, \\
\dot{q}(t) + \delta \dot{q}(t) &= 3(t + \varepsilon)^2, \\
\delta \dot{q}(t) &= [\dot{q}(t) + \delta \dot{q}(t)] - \dot{q}(t) = 3(t + \varepsilon)^2 - 3 t^2 = 6 t \varepsilon .
\end{align}
Also notice also that
\begin{align}
\dot{q}(t) &= 3 t^2 = \frac{d}{dt} t^3 = \frac{d}{dt} q(t), \\
\dot{q}(t) + \delta \dot{q}(t) &= 3(t + \varepsilon)^2 = \frac{d}{dt} (t + \varepsilon)^3 = \frac{d}{dt} [q(t) + \delta q], \\
\delta \dot{q} &= \frac{d}{dt} \delta q
\end{align}
and so more generally
\begin{align}
\delta q &= q(t + \varepsilon) - q(t), \\
\delta \dot{q} &= \dot{q}(t + \varepsilon) - \dot{q}(t) = \frac{d}{dt}[q(t + \varepsilon) - q(t)] = \frac{d}{dt} \delta q.
\end{align}
So now lets do the steps up to $\int [\frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \delta \dot{q}] dt$:
\begin{align}
\delta S(q) &= \delta \int (\frac{1}{2} m \dot{q}^2 - \frac{1}{2} m q^2)dt = \int (\frac{1}{2} m \delta \dot{q}^2 - \frac{1}{2} m \delta q^2)dt \\
&= \int (\frac{1}{2} m 2 \dot{q} \delta \dot{q} - \frac{1}{2} m 2 q \delta q)dt \\
&= \int (- m q(t) \delta q(t) + m \dot{q}(t) \delta \dot{q}(t) )dt \\
&= \int [\frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \delta \dot{q}] dt.
\end{align}
Notice in the 2nd last inequality I just re-ordered the integrand to be all $q$ terms first and then all $\dot{q}$ terms (the Lagrangian was written with all $\dot{q}$ terms first then all $q$ terms because of the ugly minus sign). Notice also in the 2nd last inequality I made explicit that everything is a function of $t$, e.g. $\delta q(t)$.
We can re-do these steps for the special example of $q(t) = t^3$ shifted into $q(t) + \delta q(t) = (t + \varepsilon)^3$:
\begin{align}
\delta S(q) &= \delta \int (\frac{1}{2} m \dot{q}^2 - \frac{1}{2} m q^2)dt \\
&= \int (\frac{1}{2} m \delta \dot{q}^2 - \frac{1}{2} m \delta q^2)dt \\
&= \int (\frac{1}{2} m 2 \dot{q} \delta \dot{q} - \frac{1}{2} m 2 q \delta q)dt \\
&= \int (- m q(t) \delta q(t) + m \dot{q}(t) \delta \dot{q}(t) )dt \\
&= \int [- m t^3 (3 t^2 \varepsilon ) + m (3 t^2) (6 t \varepsilon) ]dt \\
&= \int [- 3 m t^5 + 18 m t^3 ] \varepsilon dt
\end{align}
The next thing you want to see is $\int \limits_{t_1}^{t_2} \frac{\partial L}{\partial q} \delta q dt \neq [ \frac{\partial L}{\partial q} \delta q(t) ]_{t_1}^{t_2}$, which in our example is:
\begin{align}
\int \limits_{t_1}^{t_2} \frac{\partial L}{\partial q} \delta q dt &\neq [ \frac{\partial L}{\partial q} \delta q(t) ]_{t_1}^{t_2} \\
\int \limits_{t_1}^{t_2} (- k q) \delta q dt &\neq [ (- k q) \delta q(t) ]_{t_1}^{t_2} \\
\int \limits_{t_1}^{t_2} (- k t^3) \delta q dt &\neq [ (- k t^3) \delta q ]_{t_1}^{t_2} \\
\int \limits_{t_1}^{t_2} (- k t^3) 3 t^2 \varepsilon dt &\neq [ (- k t^3) 3 t^2 \varepsilon ]_{t_1}^{t_2} \\
\int \limits_{t_1}^{t_2} (- 3 k t^5 \varepsilon ) dt &\neq [-3 k t^5 \varepsilon ]_{t_1}^{t_2} \\
(- 3 k \frac{t^6}{6} \varepsilon )|_{t_1}^{t_2} &\neq [-3 k t^5 \varepsilon ]_{t_1}^{t_2}
\end{align}
so hopefully you see whis this was so fundamentally wrong to do. The right thing to do is
to use
\begin{align}
\delta \dot{q}(t) = (\dot{q}(t) + \delta \dot{q}) - \dot{q}(t) = \frac{d}{dt} q(t+\varepsilon) - \frac{d}{dt} q(t) = \frac{d}{dt}[q(t+\varepsilon) - q(t)] = \frac{d}{dt} \delta q
\end{align}
and then do
\begin{align}
\delta S(q) &= \delta \int (\frac{1}{2} m \dot{q}^2 - \frac{1}{2} m q^2)dt = \int (\frac{1}{2} m \delta \dot{q}^2 - \frac{1}{2} m \delta q^2)dt \\
&= \int (\frac{1}{2} m 2 \dot{q} \delta \dot{q} - \frac{1}{2} m 2 q \delta q)dt \\
&= \int (- m q(t) \delta q(t) + m \dot{q}(t) \delta \dot{q}(t) )dt \\
&= \int [\frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \delta \dot{q}] dt \\
&= \int [\frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \frac{d}{dt} \delta q] dt \\
&= \int [\frac{\partial L}{\partial q} \delta q + \frac{d}{dt} (\frac{\partial L}{\partial \dot{q}} \delta q) - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} ) \delta q ] dt \\
&= \int [\frac{d}{dt} (\frac{\partial L}{\partial \dot{q}} \delta q)] dt + \int [\frac{\partial L}{\partial q} \delta q - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} ) \delta q ] dt \\
&= (\frac{\partial L}{\partial \dot{q}} \delta q) |^{t_2}_{t_1} + \int [\frac{\partial L}{\partial q} \delta q - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} ) \delta q ] dt \\
&= (\frac{\partial L}{\partial \dot{q}} \delta q) |^{t_2}_{t_1} + \int [\frac{\partial L}{\partial q} - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} )] \delta q dt
\end{align}
If we want the $q(t)$ that we began from and then varied into $q(t) + \delta q(t)$ to be the $q(t)$ which minimizes the action $S(q)$, we see that it's first order change should be zero, analogous to how in a Taylor expansion $\delta f(x) = f(x) - f(x_0) = f'(x_0)(x - x_0) + \frac{1}{2}f''(x_0)(x - x_0)^2 + \dots$ we see $f'(x_0) = 0$ implies the function $f$ is extremised at $x_0$, and recall it doesn't matter if $f''(x_0)$ is non-zero, this tells us about the curvature around the extremum of $f$. Similarly for $S$ we could figure out higher order changes by starting from $\delta S = S(q + \delta q) - S(q)$ and Taylor expanding in $\delta q$, i.e.
\begin{align}
\delta S = S(q + \delta q) - S(q) = \int_{t_1}^{t_2} [L(t,q + \delta q,\dot{q} + \delta \dot{q}) - L(t,q,\dot{q})]dt = \int_{t_1}^{t_2} [L(t,q ,\dot{q}) + \partial_q L \delta q + \partial_{\dot{q}} L \delta \dot{q} + \frac{1}{2} \frac{\partial ^2 L}{\partial q^2} \delta q^2 + \frac{\partial^2 L}{\partial q \partial \dot{q}} \delta q \delta \dot{q} + \dots - L(t,q,\dot{q})]dt
\end{align}
but we only keep terms to first order in $\delta q$ and $\delta \dot{q}$. So if we are only keeping terms to first order and we want $S$ to be minimized at $q(t)$ we need $\delta S$ above to vanish, which means we need the first term to vanish for any variation $\delta q$, which implies we need $\delta q$ to be zero for any changes in $q(t)$, which means we can only plug in functions which are constant at the end-points $t_1$ and $t_2$, so that $\delta q(t_1) = \delta q(t_2) = 0$ no matter what. This is a huge restriction, it means the $t^3$ example we've been working with is not even allowed for our specific Lagrangian as we'll see below, though in other examples it would be fine. But once we restrict ourselves to functions which vanish at the endpoints, we also need the $q(t)$ which minimizes $S(q)$ for any $\delta q$ to satisfy the Euler-Lagrange equations, otherwise when we plug the $q(t)$ into these equations we will get something non-zero which means this integral above will give a non-zero number and so $\delta S$ is not zero and so $q(t)$ was never a minimum. We can thus view the EL equations as the conditions which a function must satisfy in order to minimize $S$ (as long as it also satisfies $\delta q(t_1) = \delta q(t_2) = 0$. These are strict conditions, and as you can see at the beginning of this post these conditions lead to Newton's law for the Harmonic oscillator starting from that Lagrangian.
As an exercise you can re-do these steps plugging in the explicit values as I did before, what you should get is
\begin{align}
\delta S(q) &= (\frac{\partial L}{\partial \dot{q}} \delta q) |^{t_2}_{t_1} + \int [\frac{\partial L}{\partial q} - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} )] \delta q dt \\
&= ( m q \delta q) |^{t_2}_{t_1} + \int [\frac{\partial L}{\partial q} - (\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} )] \delta q dt \\
&= ( m q \delta q) |^{t_2}_{t_1} + \int [- k q - (\frac{d}{dt} m \dot{q} )] \delta q dt \\
&= ( m q \delta q) |^{t_2}_{t_1} + \int [- k q - m \ddot{q} ] \delta q dt \\
&= ( m q 3 t^2 \varepsilon) |^{t_2}_{t_1} + \int [- k q - m \ddot{q} ] 3 t^2 \varepsilon dt \\
&= ( m t^3 3 t^2 \varepsilon) |^{t_2}_{t_1} + \int [- k t^3 - 6 m t ] 3 t^2 \varepsilon dt \\
&= ( 3 m t^5 \varepsilon) |^{t_2}_{t_1} + \int [- k t^3 - 6 m t ] 3 t^2 \varepsilon dt .
\end{align}
Notice the first term $( 3 m t^5 \varepsilon ) |^{t_2}_{t_1}$ is a non-zero number. In mechanics we are only allowed to use $q(t)$'s which vanish at the endpoints, so we are only allowed to use $q(t)$'s for which this vanishes, and since $t^2$ only vanishes at $t = 0$ we are not even allowed to plug this into the problem physically, though mathematically we did nothing wrong in simply plugging it in and seeing what happens, in other words it is an assumption to restrict the space of functions we are only allowed to use to functions whose variations vanish at the endpoints, in other words, a boundary condition. Lets look at the other term, lets even pretend $q(t) = t^3$ did satisfy the boundary conditions of vanishing at the endpoints for a moment (i.e. lets just ignore that it doesn't vanish at the endpoints). What is going on with
\begin{align}
\int [- k q - m \ddot{q} ] \delta q dt = \int [- k q - m \ddot{q} ] 3 t^2 \varepsilon dt = \int [- k t^3 - 6 m t ] 3 t^2 \varepsilon dt ?
\end{align}
You can literally compute this integral as well and get some number, and that's fine mathematically. If we want $q(t)$ to be a minimum however (i.e. our physics assumption) we need this to vanish, i.e. to be zero, and even more strictly, it has to vanish for all choices of $\varepsilon$, because we want to minimize $S$. This means that the $[- k q - m \ddot{q} ]$ part has to be equal to zero, and this means that we are forced to only work with the $q(t)$ that satisfies the differential equation $ k q + m \ddot{q} = 0$, but this is the Euler-Lagrange equations applied to the Harmonic oscillator Lagrangian we started with. The fact that $k q + m \ddot{q} = - k t^3 - 6 m t$ for $q(t) = t^3$ means this choice of $q$ is not a solution to the Harmonic oscillator problem.
