Question about surface integration

Hi, I'm facing a problem in calculating surface integration. I defined the function below which calculates surface integration for a discrete dataset using the trapezoidal rule.

.m

function out = trapz2D(x,y,f)

    x = x(:);
    y = y(:);

    nx = length(x);
    ny = length(y);

    dx = diff(x);
    dy = diff(y);
    dS = dy*dx.';

    df = (f(1:ny-1,1:nx - 1) + f(2:ny,1:nx-1) + f(1:ny-1,2:nx) + f(2:ny,2:nx))/4;

    out = sum(sum(dS.*df));

end

I calculated surface integration using this function and the dataset attached, and then compare it to the value obtained from "trapz" which is a MATLAB built-in function.

surfInteg.m

clc
close all
clear
format long

load('dataset.mat')

I1 = trapz(y, trapz(x,Pz3,2));
I2 = trapz2D(x,y,Pz3);

disp([I1; I2])

>>>1.0e-12 *

   0.307618158054522 - 0.000000000004792i
   0.307618158054522 - 0.000000000004792i

I1 = trapz(y, trapz(x,Pz1,2));
I2 = trapz2D(x,y,Pz1);

disp([I1; I2])

>>>1.0e-27 *

   -0.135645057561047 + 0.013248760976931i
   -0.129284381233370 + 0.013497757971006i

The results are similar for "Pz3", but not for "Pz1". Could someone explain why this is happening?

[Torsten]

Scale Pz1 by a factor of 1e27 and see what happens.

Best wishes

Torsten.