clear all
syms s1 s2 zz lamda alpha c a w;%variance and parameter
C=10^-14;
Z=50:10:300;
k=2*pi/lamda;
rho=(0.545*c^2*k^2*zz)^(-3/5);
f1=k^2/(zz^2)*s1*s2;
f2=exp(j*k/(2*zz)*(s1^2-s2^2))*exp(-1/rho^2*(s1^2+s2^2))*besseli(0,2*s1*s2/rho^2);
%上面f1,f2为被积函数中与源无关部分，下面为sgb
sgb=besselj(0,alpha*s1)*besselj(0,alpha*s2)*exp(-s1^8/w^8)*exp(-s2^8/w^8);
f=f1*f2*sgb;%integrand

f=subs(f,{zz,lamda,alpha,c,w},{Z,1550*10^-9,600,C,0.05});%这是一个矢量；
disp(f);
intValue=ones(size(Z));
for i=1:length(Z)
    disp(i);
    ff=inline(f(i));      
    intValue(i)=quad2dg(ff,0,0.1,0,0.1,10^-20);
end;
plot(Z',abs(intValue))

总是显示有奇点，请教大侠，如何处理？问题1，如何确定奇点的位置？
问题2 为了完成积分确定了位置然后怎样处理呢？
注：quad2dg可以改为dblquad，谢谢！请高手费心！