How to integrate PDEs back to the original function, e.g. $x e^{-x^2-y^2}$? (Generic or numeric approaches are welcome.)

The notation I'll use is $x = (x^1, \ldots, x^n)$ and $y = (y^1, \ldots, y^n).$ Suppose we are given a system,

$$\frac{\partial f}{\partial x^1} = g_1(x),\quad\ldots,\quad\frac{\partial f}{\partial x^n} = g_n(x)$$

of $n$ decoupled partial differential equations defined on some open rectangle[1] containing some point $x_0 \in \mathbb{R}^n.$ All this means is that there are no singularities inside the rectangle. The problem of finding $f$ amounts to asking whether there exists a smooth function $f$ whose differential satisfies,

$$df = g_1(x)\,dx^1 + \cdots + g_n(x)\,dx^n.$$

Viewing this object as a smooth one-form, a necessary and sufficient condition for this to be integrated to construct $f$ is that the exterior derivative vanishes identically on the open rectangle (see the Poincaré Lemma). That is,

$$0 = dg_1 \wedge dx^1 + \cdots + dg_n \wedge dx^n.$$

I will not expand this out completely in the general case (it is not anymore enlightening), but I will show what happens in the 2 variable case. If this condition holds, you can perform a line integral to compute $f.$ That is, for any fixed point $y$ the value of $f$ at $y$ is given by,

$$f(y) = f(x_0) + \int_{x_0}^{y} g_1(x)\,dx^1 + \cdots + g_n(x)\,dx^n,$$

where any path can be taken for the integral as long as it stays within the open rectangle where the condition is satisfied. Note, that if the conditions do not hold, there is no such $f$ that solves the equations.

Two Variable Example

Consider the smooth one-form

$$\omega = g_1(x^1, x^2)\,dx^1 + g_2(x^1, x^2)\,dx^2,$$

describing the PDE

$$\frac{\partial f}{\partial x^1} = g_1(x^1, x^2),\quad\frac{\partial f}{\partial x^2} = g_2(x^1, x^2),$$

over some open rectangle containing $x_0.$ For some intuition, you should see that this is equivalent to asking that the mixed partial derivatives of the function $f,$ if it exists, should be equal independent of the order they are taken.

The exterior derivative of $\omega$ is,

$$d\omega = \left[\frac{\partial g_2}{\partial x^1} - \frac{\partial g_1}{\partial x^2}\right] \,dx^1\wedge dx^2.$$

We need $d\omega = 0$ on an open rectangle containing $x_0,$ so we need

$$\frac{\partial g_2}{\partial x^1} = \frac{\partial g_1}{\partial x^2}$$

on that open rectangle. We can verify that your example does satisfy this. Let us check this. Let

$$g_1(x) = (1 - 2\,(x^1)^2)\,\exp(-((x^1)^2 + (x^2)^2)),\quad g_2(x) = -2\,x^1\,x^2\,\exp(-((x^1)^2 + (x^2)^2)).$$

The required partial derivatives are,

$$\begin{aligned}\frac{\partial g_1}{\partial x^2} &= -2\,x^2\,(1 - 2\,(x^1)^2)\,\exp(-((x^1)^2 + (x^2)^2)),\\ \frac{\partial g_2}{\partial x^1} &= -2\,x^2\,\exp(-((x^1)^2 + (x^2)^2)) + 4\,(x^1)^2\,x^2\,\exp(-((x^1)^2 + (x^2)^2)).\end{aligned}$$

Clearly $\frac{\partial g_1}{\partial x^2} - \frac{\partial g_2}{\partial x^1} = 0.$ As a result, we can directly compute the smooth function $f$ whose differential is $\omega.$ Explicitly, the value of $f$ at some arbitrary point $y$ in an open rectangle containing $x_0$ is,

$$f(y) = f(x_0) + \int_{x_0}^{y} (1 - 2\,(x^1)^2)\,\exp(-((x^1)^2 + (x^2)^2))\,dx^1 -2\,x^1\,x^2\,\exp(-((x^1)^2 + (x^2)^2))\,dx^2,$$

where $f(x_0) \in \mathbb{R}$ is an arbitrary initial value. For this integral, we can take any path from $x_0$ to $y$ that stays within the open rectangle where the condition is satisfied. As a result, let us consider taking the easiest path: first travel from $(x_0^1, x_0^2)$ to $(y^1, x_0^2)$ and then to $(y^1, y^2).$ The notation is a bit confusing here, but I mean perform a line integral in the $x^1$ direction and then along the $x^2$ direction. This becomes,

$$f(x) = f(x_0) + \left(\int_{x_0^1}^{y^1} (1 - 2\,(x^1)^2)\,\exp(-((x^1)^2 + (x_0^2)^2))\,dx^1\right) + \left(\int_{x_0^2}^{y^2} -2\,y^1\,x^2\,\exp(-((y^1)^2 + (x^2)^2))\,dx^2\right).$$

Note how the $x^2$ value is fixed in the first integral with the initial position's second coordinate $x_0^2$ while the $x^1$ value is fixed in the second integral with the final position's first coordinate $y^1.$ This is because of the order of the line integral I have decided to take. This method can be done both analytically (where tractable) and numerically as you did.

In this case, we can solve the problem analytically. Let us fix $x_0 = (0,0)$ with initial value of $f(x_0) = 0.$ Factor out the constants,

$$f(y) = f(x_0) + \exp(-(x_0^2)^2)\,\left(\int_{x_0^1}^{y^1} (1 - 2\,(x^1)^2)\,\exp(-(x^1)^2)\,dx^1\right) + y^1\,\exp(-(y^1)^2)\,\left(\int_{x_0^2}^{y^2} -2\,x^2\,\exp(-(x^2)^2)\,dx^2\right).$$

Now we have fairly easy single-variable integrals, so proceed with the integration,

$$f(y) = f(x_0) + \exp(-(x_0^2)^2)\,\left[x^1\,\exp(-(x^1)^2)\right]_{x_0^1}^{y^1} + y^1\,\exp(-(y^1)^2)\,\left[\exp(-(x^2)^2)\right]_{x_0^2}^{y^2}.$$

Substitute our values for $x_0$ and $f(x_0)$ to find,

$$f(y) = \left[x^1\,\exp(-(x^1)^2)\right]_{0}^{y^1} + y^1\,\exp(-(y^1)^2)\,\left[\exp(-(x^2)^2)\right]_{0}^{y^2}.$$

Simplify and find,

$$f(y) = y^1\,\exp(-(y^1)^2) + y^1\,\exp(-(y^1)^2)\left(\exp(-(y^2)^2) - 1\right) = y^1 \,\exp(-(y^1)^2 - (y^2)^2),$$

as required.

[1]: Any open, convex set works.