磁场的方程

function [Bx,By,Bz]=B_(x,y,z)

%计算磁场的三个分量

L1=0.00509;theta0=pi/6;r=0.170;R=0.8;u0=4*pi*10^-7;

I=100;

t0=0;t1=theta0;a=1;b=a*2*pi*(R-r)/L1;

syms t;

x1=r*sin(a*t);

y1=R*sin(b*t)+r*cos(a*t)*sin(b*t);

z1=R*cos(b*t)+r*cos(a*t)*cos(b*t); %以上三行为线圈的参数方程

dx=diff(x1,t);dy=diff(y1,t);dz=diff(z1,t);

r_=[x-x1,y-y1,z-z1];

dl_=[dx,dy,dz]; %将线元参数化

dB_=u0*I/(4*pi).*cross(dl_,r_)/(dot(r_,r_))^1.5; %求磁场关于t的导数

dBx=dB_(1); %取符号函数的第一项

dBx=vectorize(dBx); %变点乘除

Bx=quad(dBx,t0,t1);

dBy=dB_(2); %取符号函数的第二项

dBy=vectorize(dBy); %变点乘除

By=quad(dBy,t0,t1);

dBz=dB_(3) ; %取符号函数的第三项

dBz=vectorize(dBz); %变点乘除

Bz=quad(dBz,t0,t1);

请附带给出画一定范围的磁场所用的时间。谢谢!