Biot-Savart law application: calculating a magnetic field of a square loop at a random point

Hey! I've created this program that calculates the magnetic field at a distance d_x, d_y, d_z from a point that dists b from the coil plane, but when i substitute it for d_x=d_y=d_z=0 it is incongruent with the biot savart's approximation:

I=3.81; %current

a=0.4; % 1/2 of the side of the square coil

b=0.2; % distance between the reference point and the square coil

u_0=1.2566e-6; %permeability of free space

pi=3.14159265359;

%coil in the yz plane

B_1=[0 0 0];

y=-a;

z=a;

d_x=0; %deviation from the reference point (b; 0; 0) in x

d_y=0; %deviation from the reference point in y

d_z=0; %deviation from the reference point in z

int=[0 0 0] %this parameter will help us integrate dIxr/r^3

cro=[0 0 0]

while y<a %magnetic field due to the upper side of the square coil

r=[b+d_x d_y-y d_z-a];

cro=cross([0 1 0], r);

nor=(norm(r))^3;

int=int+cro/nor;

y=y+0.00001;

end

B=u_0*I*int*N*10000/(4*pi)

B_1=B_1 + B

z=a;

y=a;

int=[0 0 0]

while z>-a %magnetic field due to the right side of the square coil

r=[b+d_x d_y-a d_z-z];

cro=cross([0 0 -1],r);

nor=norm(r)^3;

int=int+cro/nor;

z=z-0.00001;

end

B=u_0*I*N*int*10000/(4*pi)

B_1=B_1 + B %sum of the contributions

z=-a;

y=a;

int=[0 0 0]

while y>-a %magnetic field due to the lower side of the square coil

r=[b+d_x d_y-y d_z+a];

cro=cross([0 -1 0],r);

nor=norm(r)^3;

int=int+cro/nor;

y=y-0.00001;

end

B=u_0*I*int*N*10000/(4*pi)

B_1=B_1 + B %sum of the contributions

z=-a;

y=-a;

int=[0 0 0]

while z<a %magnetic field due to the left side of the square coil

r=[b+d_x d_y+a d_z-z];

cro=cross([0 0 1],r);

nor=norm(r)^3;

int=int+cro/nor;

z=z+0.00001;

end

B=u_0*I*int*N*10000/(4*pi)

B_1=B_1 + B %sum of the contributions

B_n=norm(B_1) %norm of the magnetic field vector

could you please tell me what I'm doing wrong?

Thank you in advance!