There's no reason to think that you'll find a plane that intersects all four points in general, but you can try various other approximations. A few quick ideas:

(1) Some sort of midpoint rule where you assume that the average value of f across the region [j/n,(j+1)/n] x [k/n,(k+1)/n] is simply the average of the four corner values.

(2) Instead of assuming a plane through all four corners, assume a linear interpolation in x along the slices y = k/n and (k+1)/n, then a linear interpolation in y between those two linear fits; this results in an equation of the form a + bx + cy + dxy.

Some Googling for things like "2D numerical integration" will turn up a variety of methods. It's all variations on the idea of estimating the value of f between your data points and then integrating between the data points. Which method you choose will depend on how accurate you want your results to be, and how much work you want to do to compute the integral.

There's no reason to think that you'll find a plane that intersects all four points in general, but you can try various other approximations. A few quick ideas:

1. Some sort of midpoint rule where you assume that the average value of $f $ across the region $[j/n,(j+1)/n] \times [k/n,(k+1)/n]$ is simply the average of the four corner values.

2. Instead of assuming a plane through all four corners, assume a linear interpolation in $x$ along the slices $y = k/n$ and $(k+1)/n$, then a linear interpolation in $y$ between those two linear fits; this results in an equation of the form $a + bx + cy + dxy$.

Some Googling for things like "2D numerical integration" will turn up a variety of methods. It's all variations on the idea of estimating the value of $f$ between your data points and then integrating between the data points. Which method you choose will depend on how accurate you want your results to be, and how much work you want to do to compute the integral.