Assume we have a function $L(\vec q,\vec{\dot q}, t)$ (where $\vec q$ is a vector in some given number of dimensions, similarly for $\vec{\dot q}$) and we wish to find a condition for which a chosen $q(t)$ would maximize
$$ \int_{t_1}^{t_2} L \left( q(t),\dot q(t), t \right) \, {\rm d} t $$
If we vary the path by $\vec{\delta q}$ (with $\vec{\delta q}(t_1)=\vec{\delta q}(t_2)=0$ the change in the value this perturbation introduces is
$$ \int_{t_1}^{t_2}\frac{\partial L}{\partial \vec q}\cdot \vec{\delta q}+\frac{\partial L}{\partial \vec{\dot q}}\cdot \vec{\dot{\delta q}} \, {\rm d} t = \int_{t_1}^{t_2} \left( \frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}} \right)\cdot \vec{{\delta q}} \, {\rm d} t $$
If $\vec q(t)$ is an extrema the integral we originally introduced, this must be 0. In any ordinary case, we could easily figure out some tricks to say that $(\frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}})_i=0$ for any index $i$, giving us the standard EL equations for a system of multiple particles.
Now, lets assume that we have a constraint on our system and the only admissible choices of $\vec q(t)$ are those such that $h(\vec q(t))=0$. Now we introduce a small perturbation $\vec{\delta q}$, and see to fit the constraints that $\nabla h_{q(t)} \cdot \vec{\delta q}=0$ We again get $\int_{t_1}^{t_2}(\frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}})\cdot \vec{{\delta q}}dt=0$ for when $\vec q(t)$ is an extrema of the original integral, but this time we have less freedom with the $\vec{{\delta q(t)}}$ we can choose--this leads me to a dead end. I dont see how we can pull out the integrand in this case. If we were able to prove a similar result to the one above and say that $(\frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}})\cdot \hat{{\delta q(t)}}=0$, we could use Lagrange multipliers to rewrite this as $(\frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}})=\lambda (t) \nabla h$--the way I see this is that $(\frac{\partial L}{\partial \vec q}-\frac{d}{dt}\frac{\partial L}{\partial \vec{\dot q}})$ will always be perpendicular to the tangent plane if we can pull the ingrained above out, which shows its parallel to $\nabla h$, giving us the existence of such a function of time to scale both values to equality. I could just be oblivious to something but with the equations trapped inside the integral, I am unsure how we get to utilize Lagrange multipliers to find a solution to such problem. I appreciate any help.