Help with finding magnetic field at distance 4cm and 10cm from two parallel, current carrying wires

I am doing a project which requires us to calculate and plot the magnetic field around 2 parallel wires with equal but opposite currents. Me and my partner have managed to calculate and graph the results but I am not sure how to find the magnetic field at a distance of 4cm and 10cm above the origin. We modified an existing code for one wire and weren't totally sure what each variable did which is why this is especially difficult. Thanks for any help!

clc
close all
clear all

I=.001; %Current through each wire(mA)
u0=1;   %Permeability (Copper = 1)
Nx=13;   %Number of points along x direction
Ny=13;   %Number of points along y direction
Nz=13;   %Number of points along z direction  
N=150;  %Resolution of grid
dist=2; %Distance between the wires

Xw1=(floor(-N/2):floor(N/2));   %X coordinates of wires
Yw1=linspace(dist/2,dist/2,N);  %Y coordinates of wires
Yw2=linspace(-dist/2,-dist/2,N);%Y coordinates of wire 2
Zw1=zeros(N+1,1);               %Z coordinates of wires

xp(1:Nx)=-(Nx-1)/2:(Nx-1)/2;  %X components of the vector space
yp(1:Ny)=-(Ny-1)/2:(Ny-1)/2;  %Y components of the vector space  
zp(1:Nz)=-(Nz-1)/2:(Nz-1)/2;  %Z components of the vector space

X(1:Nx,1:Ny,1:Nz)=0;    %Initialzing the vector space for X
Y(1:Nx,1:Ny,1:Nz)=0;    %Initialzing the vector space for Y
Z(1:Nx,1:Ny,1:Nz)=0;    %Initialzing the vector space for Z

for i=1:Nx
    X(i,:,:)=xp(i);     %Creating the vector space for X
end
for i=1:Ny
    Y(:,i,:)=yp(i);     %Creating the vector space for Y
end
for i=1:Nz
    Z(:,:,i)=zp(i);     %Creating the vector space for Z
end

%*************************************************************************%
%               Calculation of field due to wire 1                        %
%*************************************************************************%
for a=1:Nx      %Loop through each iteration on X axis
for b=1:Ny      %Loop through each iteration on Y axis
for c=1:Nz      %Loop through each iteration on Z axis

for i=1:N-1     %Calculate vector components for each segment of wire
Rx(i)=X(a,b,c)-(0.5*(Xw1(i)+Xw1(i+1)));     %Radius of point from X axis
Ry(i)=(Y(a,b,c)-(0.5*(Yw1(i)+Yw1(i+1))));   %Radius of point from Y axis
Rz(i)=(Z(a,b,c)-(0.5*(Zw1(i)+Zw1(i+1))));   %Radius of point from Z axis
dlx(i)=Xw1(i+1)-Xw1(i);     %Producing differential X
dly(i)=Yw1(i+1)-Yw1(i);     %Producing differential Y
dlz(i)=Zw1(i+1)-Zw1(i);     %Producing differential Z
end

%Calculating vector components at end of wire
Rx(N)=(X(a,b,c)-0.5*(Xw1(N)+1));    %Radius of point from X axis
Ry(N)=(Y(a,b,c)-(0.5*(Yw1(N)+1)));  %Radius of point from Y axis
Rz(N)=(Z(a,b,c)-(0.5*(Zw1(N)+1)));  %Radius of point from Z axis
dlx(N)=-Xw1(N)+1;       %Producing differential X
dly(N)=-Yw1(N)+1;       %Producing differential Y
dlz(N)=-Zw1(N)+1;       %Producing differential Z

for i=1:N
Xcross(i)=(dly(i).*Rz(i))-(dlz(i).*Ry(i));  %Calculating dL X Rx
Ycross(i)=(dlz(i).*Rx(i))-(dlx(i).*Rz(i));  %Calculating dL X Ry   
Zcross(i)=(dlx(i).*Ry(i))-(dly(i).*Rx(i));  %Calculating dL X Rz
R(i)=sqrt(Rx(i).^2+Ry(i).^2+Rz(i).^2);  %Calculating radius from line       
end

Bx1=(I*u0./(4*pi*(R.^3))).*Xcross;  %Calculating X component of field 
By1=(I*u0./(4*pi*(R.^3))).*Ycross;  %Calculating Y component of field 
Bz1=(I*u0./(4*pi*(R.^3))).*Zcross;  %Calculating Z component of field 

BX1(a,b,c)=0;  %Clear array holding field data before storing next X value      
BY1(a,b,c)=0;  %Clear array holding field data before storing next Y value
BZ1(a,b,c)=0;  %Clear array holding field data before storing next Z value
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     norm(BX1(a,b,c));
%     norm(BY1(a,b,c));
%     norm(BZ1(a,b,c));
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
for i=1:N       
    BX1(a,b,c)=BX1(a,b,c)+Bx1(i); %Saving X component of field 
    BY1(a,b,c)=BY1(a,b,c)+By1(i); %Saving Y component of field 
    BZ1(a,b,c)=BZ1(a,b,c)+Bz1(i); %Saving Z component of field 
end
end
end
end

%*************************************************************************%
%               Calculation of field due to wire 2                        %
%*************************************************************************%
for a=1:Nx      %Loop through each iteration on X axis
for b=1:Ny      %Loop through each iteration on Y axis
for c=1:Nz      %Loop through each iteration on Z axis

for i=1:N-1     %Calculate vector components for each segment of wire
Rx(i)=X(a,b,c)-(0.5*(Xw1(i)+Xw1(i+1)));     %Radius of point from X axis
Ry(i)=(Y(a,b,c)-(0.5*(Yw2(i)+Yw2(i+1))));   %Radius of point from Y axis
Rz(i)=(Z(a,b,c)-(0.5*(Zw1(i)+Zw1(i+1))));   %Radius of point from Z axis
dlx(i)=Xw1(i+1)-Xw1(i);     %Producing differential X
dly(i)=Yw2(i+1)-Yw2(i);     %Producing differential Y
dlz(i)=Zw1(i+1)-Zw1(i);     %Producing differential Z
end

%Calculating vector components at end of wire
Rx(N)=(X(a,b,c)-0.5*(Xw1(N)+1));    %Radius of point from X axis
Ry(N)=(Y(a,b,c)-(0.5*(Yw2(N)+1)));  %Radius of point from Y axis
Rz(N)=(Z(a,b,c)-(0.5*(Zw1(N)+1)));  %Radius of point from Z axis
dlx(N)=-Xw1(N)+1;       %Producing differential X
dly(N)=-Yw2(N)+1;       %Producing differential Y
dlz(N)=-Zw1(N)+1;       %Producing differential Z

for i=1:N
Xcross(i)=(dly(i).*Rz(i))-(dlz(i).*Ry(i));  %Calculating dL X Rx
Ycross(i)=(dlz(i).*Rx(i))-(dlx(i).*Rz(i));  %Calculating dL X Ry
Zcross(i)=(dlx(i).*Ry(i))-(dly(i).*Rx(i));  %Calculating dL X Rz
R(i)=sqrt(Rx(i).^2+Ry(i).^2+Rz(i).^2);      %Calculating radius from line    
end

Bx2=(-I*u0./(4*pi*(R.^3))).*Xcross;  %Calculating X component of field 
By2=(-I*u0./(4*pi*(R.^3))).*Ycross;  %Calculating Y component of field 
Bz2=(-I*u0./(4*pi*(R.^3))).*Zcross;  %Calculating Z component of field 

BX2(a,b,c)=0;  %Clear array holding field data before storing next X value      
BY2(a,b,c)=0;  %Clear array holding field data before storing next Y value
BZ2(a,b,c)=0;  %Clear array holding field data before storing next Z value

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% norm(BX2(a,b,c));
% norm(BY2(a,b,c));
% norm(BZ2(a,b,c));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:N         
    BX2(a,b,c)=BX2(a,b,c)+Bx2(i); %Saving X component of field 
    BY2(a,b,c)=BY2(a,b,c)+By2(i); %Saving Y component of field 
    BZ2(a,b,c)=BZ2(a,b,c)+Bz2(i); %Saving Z component of field

end
end
end
end

%Calculate magnetic field due to both wires
BX = BX1 + BX2;  %X components
BY = BY1 + BY2;  %Y components
BZ = BZ1 + BZ2;  %Z components

BZ(:,:,4)
BZ(:,:,10)