I know, from mathematics basis, that a polynomial with one variable can be factor in function of its roots, so I can generate any one-variable polynomial from its zeros. 
But I want know if is possible to generate a two-variable polynomial giving its coordinates where value is zero? 
This problem appears on calculus of *shape functions* in **Finite Elements Method**. For example, in an one-dimensional element with four nodes $(p\_1,p\_2,p\_3,p\_4)$, I have to generate four *shape functions* that have the following requisites: 


The function $N\_i(x)$ has value $1$ at the node $p\_i$ and $0$ on others nodes.


$p\_1 \to x = -1$  

$p\_2 \to x = -1/2$  

$p\_3 \to x = 1/2$  

$p\_4 \to x = 1$ 


$N\_1(p\_1) = 1$  

$N\_1(p\_2) = N\_1(p\_3) = N\_1(p\_4) = 0$ 


From this I can find the general form of the function $N\_1$: 
$N\_1(x) = (x + 1/2)(x - 1/2)(x - 1)(C\_1)$ 


$N\_1(p\_1) = 1 \to C\_1 = -2/3$ 


$N\_1(x) = (-2/3)(x + 1/2)(x - 1/2)(x - 1)$ 


This is the $N\_1$ shape function. 
Now, I want to think this problem to **two-variables** and generate those functions with the same "algorithm". 


A simple square element with coordinates: $p\_1=(-1,-1),p\_2=(1,-1),p\_3=(1,1),p\_4=(-1,1)$.  

And I want the $N\_1(x,y)$ function that has value $1$ at $p\_1$ and $0$ on other three points.


This can be possible?



Generate a two-variable polynomial from its "roots
There a number of very well researched techniques that can be used to solve your problem in a practical sense. Most of these have come from computer graphics and computer vision where taking a set of points on a surface, possibly with noise, and trying to create a surface that either interpolates or approximates those points. For example this [presentation](http://www.cs.cornell.edu/courses/cs4620/2013fa/lectures/17surfaces.pdf) from Cornell gives you some ideas of the techniques used - bezier patches, spline surfaces, subdivision surfaces, implicit surfaces, though of course not all those are directly related to a polynomial interpolation though most relate easily to [piecewise](https://en.wikipedia.org/wiki/Piecewise) polynomial interpolation or approximation.


Generating a continuous, smooth surface that interpolates your points may not really be necessary if you take the local [interpolation](https://en.wikipedia.org/wiki/Multivariate_interpolation) viewpoint. For example if you have a function from $\mathbb{R}^2\rightarrow \mathbb{R}$ which you have sampled at a discrete set of points $(x\_i,y\_i)$ you can define the value of the function at a point $(x,y)$ to be a simple function of the nearby values - the simplest case being [nearest neighbour interpolation](https://en.wikipedia.org/wiki/Nearest-neighbor_interpolation) but going to a slightly more complex [bilinear interpolation](https://en.wikipedia.org/wiki/Bilinear_interpolation) you an achieve high accuracy assuming you have a regular grid of samples - essentially you interpolate linearly along one axis and then along the other. Note that swapping the order does not matter. This would work for your specific example:


\begin{equation}
\begin{aligned}
N\_1(x,y)&=\frac{y+1}{2}(\frac{x+1}{2}N(p\_3)+\frac{1-x}{2}N(p\_4))+\frac{1-y}{2}(\frac{x+1}{2}N(p\_2)+\frac{1-x}{2}N(p\_1))\\&= \frac{(x+1)(y+1)}{4}N(p\_3)+\frac{(1-x)(y+1)}{4}N(p\_4)\\&+\frac{(1-x)(1-y)}{4}N(p\_1)+\frac{(1+x)(1-y)}{4}N(p\_2)\\&=\frac{(1-x)(1-y)}{4}
\end{aligned}
\end{equation}


