我仔细考虑了一下，x、y、z的取值范围不宜取这么大。
这是因为，当x、y、z远大于R时，分母里的θ项基本可以忽略，积分得到的各轴向的值都接近0，而这时分母又比较大，所以与x、y、z小于R的值远远不在一格量级，画出的图基本上只有坐标原点附近才有非零值。而受限于系统资源，点距又不可能取的太密集，这样画出的图没什么实际价值。
所以，我在-R~R之间取值，画出的锥形图如下（第一个是三维图，后面三个是三视图）：
参考代码：
L = R;
I = 0.1; R = 0.1;
u0 = pi*4e-7;
dB = @(x,y,z,t) u0/4*pi * [ R*I*z*cos(t); R*I*z*sin(t); ...
R*I*(sin(t).*(y-R*sin(t))+cos(t).*(x-R*cos(t))) ] / ...
sqrt( (x-R*cos(t)).^2 + (y-R*sin(t)).^2 + z.^2 ).^3;
B = @(x,y,z) quadv(@(t)dB(x,y,z,t),0,2*pi);
[x,y,z] = meshgrid(linspace(-L,L,10));
B = arrayfun(B,x,y,z,'UniformOutput',false);
[m,n,l] = size(B);
[m,n,l] = meshgrid(1:m,1:n,1:l);
Bx = arrayfun(@(i,j,k)B{i,j,k}(1),m,n,l);
By = arrayfun(@(i,j,k)B{i,j,k}(2),m,n,l);
Bz = arrayfun(@(i,j,k)B{i,j,k}(3),m,n,l);
B = sqrt(Bx.^2+By.^2+Bz.^2);
clf reset
colordef(gcf,'black')
[cx cy cz] = meshgrid(linspace(-L,L,10));
h=coneplot(x,y,z,Bx,By,Bz,cx,cy,cz,B,0.8);
set(h,'EdgeColor', 'none');
view(35,25)
xlabel('x');
ylabel('y');
zlabel('z');
box on;
axis tight
camlight left;
camlight right;
lighting phong
