To begin, the reason the gradient of a scalar-valued function is a vector field is:
$$ \nabla \phi = \langle \frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y},\frac{\partial \phi}{\partial z} \rangle $$
which is the formula for a vector field. For each point $p$ in some region of $\mathbb{R}^3$ we assign $\nabla \phi (p)$. This makes $\nabla \phi$ a vector field.
Now, if $\vec{F} = \nabla \phi$ then $\nabla \times \vec{F} = \nabla \times \nabla \phi = 0$. However, the converse is not necessarily true. It is possible to have $\nabla \times \vec{F} =0$ throughout some domain $U \subset \mathbb{R}^3$ and yet there does not exist $\phi$ on all of $U$ such that $\nabla \phi = \vec{F}$. If you want the curl vanishing to imply the existence of the potential function $\phi$ then you also need to have a topological condition on $U$. In particular, if $U$ is simply connected then we can apply Stokes' Theorem to arbitrary surfaces in $U$ and for each surface $S$:
$$ \int_{\partial S} \vec{F} \cdot d\vec{r} = \iint_S (\nabla \times \vec{F}) \cdot d\vec{S} = 0.$$
Hence integrals of $\vec{F}$ around loops in $U$ vanish hence $\vec{F}$ is path-independent and then we can prove $\displaystyle \phi(\vec{r}) = \int_{p}^{\vec{r}} \vec{F} \cdot d\vec{r}$ is a valid formula for the construction of a potential in $U$.
In any event, if you are currently taking multivariate calculus and you were just introduced to conservative vector fields then relax. In a week or two these things should all gel together. There are a few moving pieces, but, once you see how they all fit it's really pretty. We have the following equivalence:
Suppose $U$ is an open connected subset of $\mathbb{R}^n$ then the following are equivalent
- $\vec{F}$ is conservative; $\vec{F}=\nabla \phi$ on all of $U$
- $\vec{F}$ is path-independent on $U$
- $\oint_C \vec{F} \cdot d\vec{r} =0$ for all closed curves $C$ in $U$
- (add condition $U$ be simply connected) $\nabla \times \vec{F}=0$ on $U$.
You don't have to find the integration constant immediately. Keep proceeding as follows.
After you determined that $f(x,y,z) = xy+g(y,z)$, differentiate with respect to $y$.
This gives $\frac{\partial f}{\partial y}=x+\frac{\partial g}{\partial y}=F_y=x$.
Thus, $\frac{\partial g}{\partial y}=0$, which implies that $g$ is a function of $z$ only. In turn, this means that $f(x,y,z)=xy+h(z)$.
Next, differentiate $f$ with respect to $z$.
This gives $\frac{\partial f}{\partial z}=h'(z)=F_z=z^2$.
Thus, $h(z)=\frac13z^3+C$.
Finally, $f(x,y,z)=xy+h(z)=xy+\frac13z^3+C$.
To check this, we have
$$\vec F=\nabla f(x,y,z)$$
$$=\hat x\frac{\partial f}{\partial x}+\hat y\frac{\partial f}{\partial y}+\hat z\frac{\partial f}{\partial z}$$
$$=\hat xy+\hat yx+\hat zz^2$$which completes the task!
Best Answer
Let $G$ be a rectangular domain $[0, a] \times [0, b]$ with equally spaced grid $N+1 \times M+1$. Also let $h_x = \frac{a}{N}, h_y = \frac{b}{M}$. Also let the known graient be $\mathbf{f}_{i,j}, \; i = 0,\dots,N,\;j = 0,\dots,M$. We hope that $\mathbf{f} = (f_x,f_y) $ really is a conservative field.
Approach 1. Direct gradient approximation
For every point $(x_i, y_j)$ let's write $$ \frac{\phi_{i+1,j}-\phi_{i,j}}{h_x} = \frac{(f_x)_{i,j}+(f_x)_{i+1,j}}{2},\quad i = 0,\dots,N-1,\;j = 0,\dots,M\\ \frac{\phi_{i,j+1}-\phi_{i,j}}{h_y} = \frac{(f_y)_{i,j}+(f_y)_{i,j+1}}{2},\quad i = 0,\dots,N,\;j = 0,\dots,M-1\\ $$ That's trapezoidal rule applied to $$ \frac{\partial \phi}{\partial x} = f_x\\ \frac{\partial \phi}{\partial y} = f_y $$ on every segment of your grid. There are too many equations ($2NM+N+M$) for $(N+1)(M+1)$ unknowns. That's because they are not independent and are consistent only if some condition imposed on $(f_x, f_y)$ is true. That condition is some difference approximation for conservativeness condition imposed on $\mathbf f$. Numerical integration with trapezoidal rule over any path consisting of grid segments should have the same value. If so, just drop extra equations and compute $\phi_{i,j}$ in following steps: Since $\phi$ is defined up to an additive constant, fix $\phi$ at some point: $$ \phi_{0,0} = 0. $$ Compute $\phi_{i,0}$ for $i=1,\dots,N$ using the first difference equation running from $i = 1$ to $i = N$. For each $i$ now compute $\phi_{i,j}$ running from $j = 1$ to $j = M$ using the second difference equation.
The only drawback of this method is that you would obtain different values of $\phi_{i,j}$ if you choose another order of evaluation in the case when $\mathbf f$ fails to satisfy the difference conservative condition.
Approach 2. Poisson's equation
Taking divergence of $\nabla \phi = \mathbf f$ one obtain the Poisson's equation $$ \Delta \phi = \operatorname{div} \mathbf f \text{ in } G. $$ To make a well-possessed problem for Poisson's equation one need to impose some boundary conditions on $\partial G$. The natural conditions for that would be Neumann type conditions $$ \frac{\partial \phi}{\partial n} = \mathbf n \mathbf f \text{ on } \partial G. $$ This could be solved in many ways, but applying Galerkin numerical method gives us the
Approach 3. Weak formulation
Instead of solving $\nabla \phi = \mathbf f$ let's reformulate it in a weak way: find such $\phi$ that $$ \int_G \nabla \phi \nabla v\, dx dy = \int_G \mathbf f \nabla v\, dx dy $$ for every smooth function $v$ defined on $G$.
Let's take $\phi(x,y)$ as a continuous function, bilinear on every cell of the grid (that is $[x_i, x_{i+1}]\times [y_j, y_{j+1}]$). The values $\phi_{i,j}$ uniquelly define such function ($(N+1) \times (M+1)$ unknowns). Take $v(x,y)$ in exactly the same form. There are $(N+1)\times(M+1)$ linearly independend functions of this kind, thus giving us $(N+1)\times(M+1)$ equations to determine $\phi_{i,j}$. Let $\psi_{i,j}(x,y)$ be the function equal to $1$ at $(x_i, y_j)$ and equal to $0$ at every other grid point (but not at the every point, since $\psi_{i,j}(x,y)$ is continious. Note that $v_{i,j}$ support is $[x_{i-1},x_{i+1}]\times[y_{j-1},y_{j+1}]$, so integrating can be performed only on $[x_{i-1},x_{i+1}]\times[y_{j-1},y_{j+1}]$ instead of G. $\psi_{i,j}(x,y)$ is a finite element basis for $\phi(x,y)$ and $v(x,y)$, particulary $$ \phi(x,y) = \sum_{i=0}^N\sum_{j=0}^M \phi_{i,j} \psi_{i,j}(x,y) $$ Taking $\psi_{l,m}$ as different $v$ for weak formulation we obtain $$ \sum_{i=0}^N\sum_{j=0}^M \phi_{i,j} \int_G \nabla\psi_{i,j} \nabla\psi_{l,m} dxdy = \int_G \mathbf f \nabla\psi_{l,m} dxdy $$
This is a huge system of linear equations for $\phi_{i,j}$. The matrix of the system, given by $$ A_{(i,j),(l,m)} = \int_G \nabla\psi_{i,j} \nabla\psi_{l,m} dxdy $$ has nonzero elements only if $(i,j)$ is not far from $(l,m)$, precisely $|i-l| \leq 1$ and $|j - m|\leq 1$. Thus the matrix is sparse. Elements of the matrix can be evaluated by hand or using some simple quadrature rule, precise for quadratic functions. Same is true for the right hand side.
The last to say is that the matrix is singular. That's because there is no unique solution to the problem, there are infinitely many. To resolve the problem one should fix some $\phi_{i,j}$, for example $\phi_{0,0} = 0$. That removes one unknown from the system and one equation at $(0,0)$. The remaining system is nonsingular.