之后的运行又出现：Subscript indices must either be real positive inte

问题出在Q(i)上,你的i前面未赋值,被理解为虚数单位,当然不能用作数组下标.

看你的意思,似乎是想取不同的f和r,分别计算该参数条件下的重积分.
这个可以做,涉及到dblquad传递附加参数的问题.但问题又来了：在integrnd函数中,f是直接取值,而r则是由前面的数据算出来的,那么r究竟用哪个数?

帮你大致改了改,方法基本上没问题,但参数及公式需要你自己认真核对.

把以下代码全部保存到一个M文件运行即可.
function zd
[f,r]=meshgrid(500:10:1500, 0.015:0.00020:0.035);
for ii=1:length(f)
for jj=1:length(r)
a1=0.01;%孔半径
ua=1;%空气柱振速
w=2*pi.*f(ii);
c=340;%声速
k=w./c;
%a=0;%位置点的x坐标
b=0.015;%位置点的y坐标
%c1=0.0315;%位置点的z坐标（见笔记本）
o=abs(acos(b./r(jj)));%位置矢量与孔中心线之间的弧度
Q = dblquad(@integrnd, 0, 0.01, 0, 2*pi,[],[],f(ii),r(jj));
H(ii,jj)=abs(j.*k.*a1.*c.*ua.*abs(Q)./(2*pi));
end
end
mesh(f,r,H)
function z = integrnd(x, y, f, r)
% f=1000;
a1=0.01;%孔半径
ua=1;%空气柱振速
w=2*pi*f;
c=340;%声速
k=w/c;
a=0;%位置点的x坐标
b=0.01;%位置点的y坐标
c=0.0315;%位置点的z坐标（见笔记本）
o=atan(abs(c)/sqrt(a^2+b^2));%位置矢量与孔中心线之间的弧度
% r=sqrt(a^2+b^2+c^2);%点距
z = x .* exp(-j*k*(r*sqrt(1-2*x*cos(y)*sin(o)/r+x.^2/r^2))) ./ (r*sqrt(1-2*x*cos(y)*sin(o)/r+x.^2/r^2));