Numerically integrating exact differentials – Mathematica Stack Exchange

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?