quiver3 not plotting my arrays

Hi everyone, I have written some code which I hope will give me nice quiver3 plot. However when I plot the data it only displays two arrows in the figure, despite the fact almost all of the elements of the array are non zero. The code is quite long but here it is:

if true
format long e
nele = 5;
M = 4.6.*(10.^25) ; n = 3.*(10.^6) ; mp = 1.67.*(10.^-27) ; v = 400.*(10.^3); mu_0 = 4.*pi.*10.^(-7);
r_0 = (mu_0.*(M.^2)/(2.*n.*mp.*(v.^2))).^(1./6);
phi=linspace(-pi/2,pi/2,nele); theta=linspace(-pi/2,pi/2,nele);
[phi, theta]=meshgrid(phi,theta);
[X,Y,Z] = sph2cart(phi, theta, r_0)
dx=-gradient(X); dy=-gradient(Y); dz=-gradient(Z);
d=[dx,dy,dz]; %gradient vectors
d = reshape(d,nele,nele,3);
[Nx,Ny,Nz]=surfnorm(X,Y,Z);
N=[Nx,Ny,Nz];
j1 = cross(Nreshape,d);
Nreshape = reshape(N,nele,nele,3);
i = j1(:,:,1); j = j1(:,:,2); k = j1(:,:,3);
j1 = [i,j,k];
Jx = ((Ny.*Bgz)-(Nz.*Bgy)).*cos(theta).*cos(phi); %x component of current
Jy = -((Nx.*Bgz)-(Nz.*Bgx)).*cos(theta).*cos(phi); %y component of current
Jz = ((Nx.*Bgy)-(Ny.*Bgx)).*cos(theta).*cos(phi); %z component of current
Xc = ones(size(X));
Xv = kron(X,Xc) - kron(Xc,X);
Yc = ones(size(Y));
Yv = kron(Y,Yc) - kron(Yc,Y);
Zc = ones(size(Z));
Zv = kron(Z,Zc) - kron(Zc,Z);
rPrime = [Xv, Yv, Zv]; %625x3
rPrime = reshape(rPrime, nele.*nele, nele.*nele, 3); %nele.^2 x nele.^2
Jxv = repmat(Jx,nele);
Jyv = repmat(Jy,nele);
Jzv = repmat(Jz,nele);
J = [Jxv, Jyv, Jzv];
% we reshape J so we can use it in the cross function later
J = reshape(J,nele.*nele,nele.*nele,3);
% We now cross J with rPrime  
JXrPrime = cross(J,rPrime)
ii = JXrPrime(:,:,1);
jj = JXrPrime(:,:,2);
kk = JXrPrime(:,:,3);
rMagPriCubed = (((Xv.^2)+(Yv.^2)+(Zv.^2)).^1.5); % |r'|^3
b_ci = ii./rMagPriCubed; % x array of magnetic field elements
b_cj = jj./rMagPriCubed; % y array of magnetic field elements
b_ck = kk./rMagPriCubed;
b_ci(isnan(b_ci))=0; b_cj(isnan(b_cj))=0; b_ck(isnan(b_ck))=0;
sum_b_ci =mat2cell(b_ci, ones(1, size(b_ci,1)./nele).*nele, ones(1, size(b_ci,2)./nele).*nele);
Bcx = cellfun(@(c) sum(c(:)), sum_b_ci); % x array of summed magnetic field elements
sum_b_cj =mat2cell(b_cj, ones(1, size(b_cj,1)./nele).*nele, ones(1, size(b_cj,2)./nele).*nele);
Bcy = cellfun(@(c) sum(c(:)), sum_b_cj); % y array of summed magnetic field elements
sum_b_ck =mat2cell(b_ck, ones(1, size(b_ck,1)./nele).*nele, ones(1, size(b_ck,2)./nele).*nele);
Bcz = cellfun(@(c) sum(c(:)), sum_b_ck);
quiver3(X,Y,Z,Bcx,Bcy,Bcz);
 % code
end

As you can see when this is plotted, there are only two arrows at the 'poles' of the hemisphere. However if you type in Bcx, Bcy or Bcz into the command window you will see arrays which are almost completely non-zero. At first i thought this might be due to the small size of some of the elements in the array so I added a scaling factor. First 10.^20 and then 10.^30 however this still left me with the same plot.

Any ideas or suggestions would be most welcome. If I have made anything unclear do not hesitate to ask me to clear it up. Thanks Again.

[Chad Greene]

That a mighty lot of code, yet it doesn't provide the ability to replicate your result. Can you separate the wheat from the chaff and just give us the X,Y,Z,Bcx,Bcy,Bcz values that can be used to replicate the problem?