求MATLAB大神帮我看看程序哪有问题，谢谢！

我把步长改了改, 因为不然电脑运行太慢不出结果, 所以改长了, 但是没有其他影响. 下面的代码是没改的.

t=0.1;
R=0.5;
a=2; %D的取值区间上限;
b=0.04; %D的步长, 腊物我取的0.1; 
c1=0; %O取值区间下限;
c2=10; %O取值区间上限;
d=0.1; %O的步长, 我取的0.5;
P21=1;r21=0.1;
syms y1;
g=zeros(((c1-c2)/d+1),(a/b*2+1));
n=0;
for O=c1:d:c2
if (O-fix(O)==0)
clear maplemex;
end
n=n+1;
m=0;
for x=-a:b:a
m=m+1;
Rs=(O^2*(R+2*r21))/(4*(x+0.002*y1)^2+(R+2*r21)^2); 
D=4*Rs^2+R^2+4*Rs*P21+2*P21*R+P21^2; 
s1=0.5*(-(2*Rs+P21+R)+sqrt(D)); 
s2=0.5*(-(2*Rs+P21+R)-sqrt(D));
r11=Rs*(Rs+P21)/sqrt(D)*(1/(Rs+s2)*(exp(-Rs*t)-exp(s2*t))-1/(Rs+s1)*(exp(-Rs*t)-exp(s1*t)))+exp(-Rs*t); %ρ11
r22=Rs/sqrt(D)*(exp(s1*t)-exp(s2*t)); 
r12=0.5*O*(R/2+r21)*(r11-r22)/((x+0.002*y1)^2+(R/2+r21)^2); 
f=2.4*10^(-8)*y1.^2.*exp(-4.83*10^(-6).*y1.^2).*r12; %得到被积函数
f=@(y1)subs(f)+0.*y1;  % 这里加了一项0, 就能运行了, 结果肯定不受影响差局虚.
g(n,m)=quad(f,0,100); 
end
end
y=c1:d:c2;
x=-a:b:a;
Z=g;
surf(x,y,Z);

图像是精简之后的, 不过你画出来应该是类似的虚燃, 仅仅是更精细:

原因是f是你的函数的名字，又是你函数里变量的名字，导致冲突。你把你的函数改个名字就可以了。