I have a dynamical system which I am modelling using Lagrangian analytical mechanics. The Euler-Lagrange equations for my $N$ generalized coordinates can be arranged into a (linearized with respect to time) matrix system
$$ \mathbf{A} \ddot{\mathbf{x}} = \mathbf{b} $$
to solve for the generalized accelerations $\ddot{\mathbf{x}}$, which are then readily integrated to produce generalized velocity $\dot{\mathbf{x}}$ and generalized position $\mathbf{x}$. I'm solving this numerically using Python, and it is all hunky dory and beautiful.
However, I'm running into some difficulty figuring out how to configure the equations in matrix form once I add holonomic constraints. Let's call my constraints $f_1(\mathbf{x})$ and $f_2(\mathbf{x})$, which encapsulate the two constraints in the form
$$ f_1(\mathbf{x}) = 0 \\ f_2(\mathbf{x}) = 0 $$
Let's also assume that $f_1(\mathbf{x})$ and $f_2(\mathbf{x})$ are fiendishly difficult to linearize algebraically a priori, because they're complicated and there's no obvious point about which to linearize them (i.e., $\mathbf{x}$ is likely to spend a lot of time nowhere near any such state you guess to linearize about).
I understand that with constraints, my Euler-Lagrange equations will look like this, where $L$ is my Lagrangian and $x_i$ are the scalar elements of vector $\mathbf{x}$:
$$ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x_i}}\right) - \frac{\partial L}{\partial x_i} - \lambda_1 \frac{\partial f_1}{\partial x_i} - \lambda_2 \frac{\partial f_2}{\partial x_i} = 0 $$
$\frac{\partial f_1}{\partial x_i}$ and $\frac{\partial f_2}{\partial x_i}$ here are not difficult to calculate. $\lambda_1$ and $\lambda_2$ here are my Lagrangian multipliers, and represent 2 new unknowns to solve for. These Lagrangian multipliers, as unknowns, have to join the matrix $\mathbf{A}$ on the left hand side (LHS) of the system of equations.
What I can't seem to figure out, in a practical situation like this where I need to solve the system numerically, rather than the analytical convenience of textbook problems, is how to incorporate the equations of constraint into the matrix system. The new system looks something like this:
$$ \begin{bmatrix}\mathbf{A} & -\nabla{f_1} & -\nabla{f_2} \\ 0 & ?? & ?? \\ 0 & ?? & ?? \end{bmatrix} \begin{bmatrix}\ddot{\mathbf{x}} \\ \lambda_1 \\ \lambda_2 \end{bmatrix} = \begin{bmatrix}\mathbf{b} \\ 0 \\ 0 \end{bmatrix} $$
What should the ?? be here? I guessed perhaps they should be $f_1$ and $f_2$ themselves, since the equations of constraint imply that
$$ \lambda_1 f_1(\mathbf{x}) = 0 \\ \lambda_2 f_2(\mathbf{x}) = 0 $$
for general $\lambda_1, \lambda_2 \ne 0$, but I can't find any reference saying so. If that were correct, then the matrix system would look like
$$ \begin{bmatrix}\mathbf{A} & -\nabla{f_1} & -\nabla{f_2} \\ 0 & f_1 & 0 \\ 0 & 0 & f_2 \end{bmatrix} \begin{bmatrix}\ddot{\mathbf{x}} \\ \lambda_1 \\ \lambda_2 \end{bmatrix} = \begin{bmatrix}\mathbf{b} \\ 0 \\ 0 \end{bmatrix} $$
ChatGPT 3.5 told me I have to linearize my equations of constraint about some chosen state of $\mathbf{x}$, but I don't see how that will get them into a usable matrix form since the equations of constraint don't themselves contain $\lambda_1$ and $\lambda_2$.