麻烦各位帮忙修改一下程序啦

f=125;
c0=344;
r0=0.2;
pa=2;
k=2*pi*f./c0;
p0=0.00002;
m=1;
for r=eps:0.01:0.5  
        n=1;
    for theta=eps:0.01:pi
         pp=0;
           ps=0;
            for l=1:100   % l取0到100                           
                 Pl=legendre(l,cos(theta)); %勒让德函数
                   y1=sqrt(pi./2./(k.*r+eps)).*besselj(l+0.5,(k.*r+eps));%第l阶球贝塞尔函数
                   pp(n,m)=pp+pa.*(2*l+1).*(-1i).^l.*y1.*Pl; %计算入射声压
                    y2=sqrt(pi./2./(k*r+eps)).*besselj(l+0.5,k*r+eps)-1i.*sqrt(pi./2./(k*r+eps)).*bessely(l+0.5,(k.*r+eps));%第二类球汉克尔函数
                    dy1dx=diff(y1);
                     dy2dx=diff(y2);
                     Bl=(-pa).*(-1i).^(l+1).*(2*l+1)./(1i).*subs(dy1dx,'r','r0')./subs(dy2dx,'r','r0'); %求Bl,即分别对第l阶球贝塞尔函数和第二类球汉克尔函数进行求导后做商
                     ps(n,m)=y2.*Bl.*Pl;%计算散射声压
             end
            p(n,m)=pp+ps;%计算总声压
           n=n+1;
    end
   m=m+1;
end  
SPL=20*log10(abs(p)/p0);
imagesc(x,y,SPL);
colorbar
title('头散射后的平面波')
xlabel('theta')
ylabel('r(m)')