I have familiarity with lagrange interpolation polynomials in 1-D for a scalar variable of interest. I am currently interested in a 3-D interpolation of a variable that is a vector $\vec{u} = [u, v, w]$.
So essentially, I think what I need is 3 separate polynomials, one for each scalar of vector $\vec{u}$. But I am unsure how to generate a multi-variate lagrange polynomial. Any suggestions?
So let's assume you have a set of data for the unknown model of $\vec{u}$, namely $\mathcal{D} = \lbrace (x_i,y_i,z_i,\vec{u}_i) \rbrace_{i=1}^{n}$. We can then approximate the model using Lagrange Interpolation in the following way:
\begin{align} \vec{u}(x,y,z) \approx \sum_{i=1}^{n} \vec{u}_i \mathcal{L}_{i}(x,y,z) \end{align}
where $\mathcal{L}_i(\cdot)$ represents the $i^{th}$ Lagrange Polynomial and $\mathcal{L}_{i}(x_j,y_j,z_j) = \delta_{i,j}$ where $\delta_{ij}$ is the Kronecker Delta. We can then define the $i^{th}$ Lagrange Polynomial to be the following such that the Kronecker Delta property is satisfied:
\begin{align} \mathcal{L}_{i}(x,y,z) = \prod_{k=1,k \neq i}^{n} \frac{\left(x - x_{k}\right)\left(y - y_{k}\right)\left(z - z_{k}\right)}{\left(x_i - x_{k}\right)\left(y_i - y_{k}\right)\left(z_i - z_{k}\right)} \end{align}
The above formulation should allow you to tackle your problem. Note, as written in the comments, that multidimensional Lagrange Interpolation isn't the safest approach as it can be badly behaved if any of the components of the data are close to each other. This is a negative property of Lagrange Interpolation in general, even in 1D, that risks happening more often when dealing with datasets in higher dimensions.
For higher dimensional interpolation, it would be better to use a different approach.