Bilinear interpolation

Question

I need to find the coefficients $a_{jk}$ given the following equation: $$ f(x,y)=\sum_{j,k=0}^Na_{jk}x^jy^k$$ in which are known $(N \times N)$ values of the function $f(x,y)$. What is the best method to find them? In the case of $f(x)$, we can write: $$f(x)=\displaystyle\sum_{j=0}^Na_jx^j$$ so if we know N values of $f(x)$ we get the matrix equation: $$ \begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \ \\ f(x_N) \end{pmatrix}= \begin{pmatrix} 1 & x_1 & \cdots & x_1^N \\ 1 & x_2 & \cdots & x_2^N \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_N & \cdots & x_N^N \end{pmatrix}\begin{pmatrix} a_0 \\ a_1 \\ \vdots \ \\ a_N \end{pmatrix} $$ and inverting the coefficient matrix we can get the coefficient vector. How to get the coefficient matrix $a_{jk}$ in $2D$ case? Thanks.