Calculate Divergence of B field

[tiwwexx]

Hello! I'm trying to compute the divergence of the B field in a simulation to double check that it is 0.
To do this I export out an overall field sensor B(x,y,z,f0) to matlab, read it in, use the cell centers convention
(axis_n(1)+axis_n(2))/2
to reshape the data into an ndgrid

x_pos= (a.Axis0(1:end-1)+a.Axis0(2:end))/2;
y_pos= (a.Axis1(1:end-1)+a.Axis1(2:end))/2;
z_pos= (a.Axis2(1:end-1)+a.Axis2(2:end))/2;
[x,y,z] = ndgrid(x_pos,y_pos,z_pos);

and a 4D array (x,y,z,3) where the 4th dimension is the Bx,By,Bz components.

B_field =  reshape(a.Snapshot0,numel(x_pos),numel(y_pos),numel(z_pos),[]);

I then get the divergence of my b field as

div_b = divergence(y,x,z,B_field(:,:,:,1),B_field(:,:,:,2),B_field(:,:,:,3))

When looking at this field (shown below) it is clearly not 0 everywhere which somewhat confused me. I verified that my B_flux through different square and spherical surfaces was equal to 0 in the analysis tab which I thought should directly translate to div_b ~ 0 everywhere.

Any suggestions on why this didn't work would be much appreciated!