I would like to be able to numerically integrate an exact differential 2-form. Mathematica can do this symbolically using DSolve (though I don’t know how to insert an initial condition), but I get a “system is overdetermined” error when I try to do this numerically using NDSolve.

For instance, consider the following symbolic example posted by bbgodfrey on this forum:

```
DSolveValue(
{D(f(x, y), x) == 1 + 2 E^(x^2 + y^2) x Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
D(f(x, y), y) == 3 y^2 + 2 E^(x^2 + y^2) y Cos(x + y) - E^(x^2 + y^2) Sin(x + y)},
f(x, y), {x, y}) // FullSimplify
output: x + y^3 + C(1) + E^(x^2 + y^2) Cos(x + y)
```

That works fine… but there are two problems. The first (smaller) problem is that putting in the initial condition f(0,0)=0 (which forces c(1) to be -1) doesn’t work:

```
DSolveValue({D(f(x, y), x) ==
1 + 2 E^(x^2 + y^2) x Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
D(f(x, y), y) ==
3 y^2 + 2 E^(x^2 + y^2) y Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
f(0, 0) == 0}, f(x, y), {x, y}) // FullSimplify
output: DSolveValue({D(f(x, y), x) ==
1 + 2 E^(x^2 + y^2) x Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
D(f(x, y), y) ==
3 y^2 + 2 E^(x^2 + y^2) y Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
f(0, 0) == 0}, f(x, y), {x, y}) // FullSimplify
```

But the second (more important) problem is that I can’t get NDSolve to *numerically* solve the exact PDE (with the initial condition f(0,0)=0):

```
NDSolveValue({D(f(x, y), x) ==
1 + 2 E^(x^2 + y^2) x Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
D(f(x, y), y) ==
3 y^2 + 2 E^(x^2 + y^2) y Cos(x + y) - E^(x^2 + y^2) Sin(x + y),
f(0, 0) == 0}, f, {x, -2, 2}, {y, -2, 2})
output: NDSolveValue::overdet: There are fewer dependent variables, {f(x,y)}, than equations, so the system is overdetermined.
```

Can anyone help?