Calculate Divergence of B field

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 ndgridx_pos= (a.Axis0(1:end1)+a.Axis0(2:end))/2; y_pos= (a.Axis1(1:end1)+a.Axis1(2:end))/2; z_pos= (a.Axis2(1:end1)+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!