Metodo trapezi e cicli for

Ciao a tutti.
Dopo 10 anni di inutilizzo ho rispolverato Matlab per la simulazione di un campo magnetico di un particolare induttore. Cerco di spiegare brevemente.

Questa è una spira percorsa da corrente, in un sistema (x,y,z).

Queste sono le equazioni che determinano il campo magnetico in un generico punto:

Le voglio integrare con il metodo dei trapezi per avere ciascuna componente e poi sommarle vettorialmente per avere una matrice 3D con i valori calcolati nei vari punti. E fin qui ci siamo.

Però non ho una spira sola, ne ho più di una, quindi dovrei ripetere il tutto più volte e sommare il contributo di ognuna, dando in ingresso i parametri che mi servono (numero di spire, raggio di ognuna e corrente circolante).

Ho messo gli input che mi servono però alla fine non riesco a sommare il contributo di ogni spira. Dovrei mettere un ulteriore ciclo per arrivare in 4D e poi sommare il contributo di ogni riga di questa matrice 4D?

Vi copio il codice che sto provando. Oltretutto, il surf in alcuni punti viene una schifezza e non mi rimane nemmeno simmetrica, cosa che mi sarei aspettato.

clear all
clc
clf

n=10;
u=10;
x=linspace(-0.1,0.1,u);  
y=linspace(-0.1,0.1,u);
z=linspace(-0.1,0.1,u);  % tre vettori spaziali

u0=4*pi*10^(-7); % costante
Spire=input('Inserire il numero di spire: '); % numero di spire come input 
I=input('Inserire corrente circolante [A]: '); % corrente come input
for q=1:Spire
    R(q)=input('Inserire Raggio [mm]: '); % ciclo per farmi chiedere tante volte il raggio a seconda di quante sono le spire
    r(q)=R(q)/1000;
    angolo=[0:2*pi/n:2*pi];
Bx_Tot=zeros(length(x),length(y),length(z)); 
By_Tot=zeros(length(x),length(y),length(z));
Bz_Tot=zeros(length(x),length(y),length(z)); % inizializzo le matrici 3D di ogni componente
B=zeros(length(x),length(y),length(z));
    for j=1:length(x)
      for w=1:length(y)
        for t=1:length(z)
          for i=1:(length(angolo)-1)   % ciclo sui tre assi e su ogni settore dell'integrale
Bx=u0*I*r(q)*z(t)/4/pi*(2*pi/(2*n))*((cos(angolo(i))/(x(j)^2+y(w)^2+z(t)^2+r(q)^2-2*r(q)*(x(j)*cos(angolo(i)+y(w)*cos(angolo(i)))))^(3/2))); % metodo dei trapezi
Bx_Tot(j,w,t)=Bx_Tot(j,w,t)+Bx; % sommo il contributo di ogni ciclo per avere l'area totale (ossia l'integrale completo)
By=u0*I*r(q)*z(t)/4/pi*(2*pi/(2*n))*((sin(angolo(i))/(x(j)^2+y(w)^2+z(t)^2+r(q)^2-2*r(q)*(x(j)*cos(angolo(i)+y(w)*cos(angolo(i)))))^(3/2))); % metodo dei trapezi
By_Tot(j,w,t)=By_Tot(j,w,t)+By; % sommo il contributo di ogni ciclo per avere l'area totale (ossia l'integrale completo)
Bz=u0*I*r(q)/4/pi*(2*pi/(2*n))*((r(q)-y(w)*sin(angolo(i))-x(j)*cos(angolo(i)))/(x(j)^2+y(w)^2+z(t)^2+r(q)^2-2*r(q)*(x(j)*cos(angolo(i))+y(w)*cos(angolo(i))))^(3/2)); % metodo dei trapezi
Bz_Tot(j,w,t)=Bz_Tot(j,w,t)+Bz; % sommo il contributo di ogni ciclo per avere l'area totale (ossia l'integrale completo)

B_tot(j,w,t)=sqrt(Bx_Tot(q).^2+By_Tot(q).^2+Bz_Tot(q).^2); % somma vettoriale dei tre contributi per avere una matrice 3D con il valore assoluto ottenuto
% questo valore qui sopra va bene con una spira. Devo però sommare il
% contributo di ogni spira (il valore l'ho passato in input)
          end
        end
    end
    end


B=zeros(q,length(x),length(y),length(z));  % tentativo di fare un ciclo per ottenere una matrice 4D con le q righe che sono il numero di spire
B(q) = (sqrt(Bx_Tot(q).^2+By_Tot(q).^2+Bz_Tot(q).^2)); % credo di ottenere però solo il contributo dell'ultimo ciclo e non sommare il resto

   
end
% %  
for f=1:length(z)
    % surf(x,y,B_Tot)
    figure(f)
surf(x,y,B_tot(:,:,f))
figure(f+length(z))
plot(x,B_tot(:,1,f))
grid on
end

Facendo il programma per una singola spira (trascurando cioè la parte del numero di spire), ottengo due cose strane. La prima sono numeri complessi nella matrice 3D che non mi giustifico (anche se con parte immaginaria sempre =0). La seconda risultati sballati dell'equazione.
Ad esempio, su (x,y) fissando z=1 ottengo un SURF di questo tipo:

Che è perfetto direi ( anche se stranamente non è simmetrico)
Lo stesso ottengo con Z ai valori più alti (ossia inizio e fine vettore z).
Per valori centrali di Z, tipo (x,y,z) = (:,:,6) ottengo una schifezza del genere:

Qualcuno sa illuminarmi?
Grazie mille