matlab符号积分计算很慢,三重积分使用定义用数值计算，for循环太多速度太慢，怎......

对一个函数进行三重积分，函数在附件中。

在球坐标系中对这个式子积分，按照剖分定义不断叠加for循环，因剖分次数太多，所以电脑一直出不来结果。想利用矩阵化代替for，写了以后出现了问题。贴出原程序和矩阵化以后的程序，想请大家帮我看看下一步要怎么进行。clear

clc

bx=0;by=0;

b_x=0.05;b_y=0.05;

Nx=49;Ny=49;

f1=1200;

lamda1=340/f1;

k1=2*pi/lamda1;

f2=1e9;

lamda2=3e8/f2;

k2=2*pi/lamda2;

Pt=1000;%w

Pt0=10*log(Pt);

thetav=0.75*pi/180;%天线波束宽度

ge=40;%dB,增益

A0=1;

Z0=A0^2/(4*pi*Pt*ge);%阻抗

%%%%%%%%%%%%%%%环境因子

T=288.15;  %K

e=10.02;   %hpa

P=1013.25;  %hpa

W=10000;      % w   功率

rho_0=1.29; % Kg/m^3   空气密度

c_0=331.6+0.6*15.15;    % 声音速度  温度15.15摄氏度

r_0=5*10^(-3);          % 单个声源半径

A=5*10^(-4);            % 声波传播路径上的衰减

%%%%%%%%%%%%%%%%

theta2=linspace(2*pi/180,179*pi/180,300);%散射角

% theta2=[179*pi/180,119*pi/180,59*pi/180,9*pi/180];

H0=10;%不均匀体内参考点高度

distance=2*H0./tan(theta2/2);%距离

R0=(distance/2)./cos(theta2/2);%波传播“标准路径”

gama=(distance/2)/(2*R0);

radius=R0*thetav/(4*sin(2*gama));%球域的半径

H=H0+radius;%球中心与坐标原点的距离

Z=radius./H;

fei1=linspace(0,2*pi,200);

%%%%%%%%%%%%%%%%

h=2*pi/200;

F=zeros(1,300);

E01=zeros(1,300);

Q1=[];

Q2=[];

Q3=[];

Q4=[];

G=[];

for m=1:length(theta2)

theta1=linspace(0,asin(Z(m)),200);

g=asin(Z(m))/20