codes runs too long

Hi, i wonder if there is somebody who can help. I have this program where it calculates pressure inside a cavity. For normal case, the program just took 30 minutes to give result. I did modification to the program where the concerned codes a shown below, could this codes are the culprit where it took 15 hrs to run and give result. a, b, c, d are bessel functions that i dont bother to write here. Note: this is just partial of the program, i doubt other part is the reason because i only change numeric values for them.

HHvp=1e-10;
HHvp1=1e-10;
HHvp2=1e-10;
HHvp3=1e-10;
HHvp4=1e-10;
HHvp5=1e-10;

     for m1=1:5
         for n1=1:5
             for l1=1:5
              omega_n(m1,n1,l1)=m1*n1*l1 % a lot more to this, to long to write)
             end
         end
     end

     for n1=1:5
       for l1=1:5
          HHvp1= HHvp1 + a*n1*l1
       end
     end

     for n1=1:5
       for l1=1:5
          HHvp2= HHvp2 + b*n1*l1
       end
     end

     for n1=1:5
       for l1=1:5
          HHvp3= HHvp3 + c*n1*l1
       end
     end

     for n1=1:5
       for l1=1:5
          HHvp4= HHvp4 + d*n1*l1
       end
     end        

     HHvp=HHvp1+HHvp2+HHvp3+HHvp4;

     Hvp(mmm2,nnn2,lll2,num2)= HHvp;

[EDITED, Jan, code formatted]

[Azzi Abdelmalek]

What computer are you using? It took 0.02 s with a simple Pentium dual core (2.3Ghz).

[Jan]

@Azzi: Have you seen the comment "% a lot more to this, to long to write)"? Perhaps the author abbreviated the code and expect us to be able to imagine the omitted parts.

@zamri: We cannot help to improve the code, when we see only parts of it. Did you use the profiler already to find the most time-consuming lines? Unfortunately the profiler disables the JIT acceleration, such that its measurements are strongly biased when loops are used. Therefor I recommend some TIC/TOC measurements to determine the runtimes of the different parts.

[Azzi Abdelmalek]

Simon, I did'nt see the comment, maybe because the question was not well formated. But what is the aim to post a part of code that took 0.02 s?

[Jan]

@Azzi: I hope that the OP will explain this.

Answer 1

[Jan]

In reply to your posted code:

When you use the profile to find the most time consuming lines, a standard procedure gets obvious: Avoid repeated calculations inside a loop. Here the Bessel function waste a lot of time and using temporary variables accelerate the code by a factor > 100:

r=1e-10+(lll2-1)*dr;  % After this line insert this:

tmp1 = 7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r);
tmp2 = -22*besselj(0,96.7*r)-18.2*bessely(0,96.7*r);
tmp3 = -3*besselj(0,145*r)-17.5*bessely(0,145*r);
tmp4 = -0.8*besselj(0,193*r)-1.6*bessely(0,193*r);
tmp5 = -1.8*besselj(0,242*r)-0.37*bessely(0,242*r);

% And now replace the correspondig calculations by "tmpX" using the
% Replace dialog of the editor.

Btw., the code is extremely ugly. A substantial cleanup would at first increase the speed and at second allow a successful debugging. Most of all consider the MLint warnings (see the orange bars in the Editor). Then add comments and split the very large lines into readable pieces. Then it will be much easier to identify the parts, which are calculated repeatedly. Finally I guess, that a speedup of factor 1000 is possible without any magic or sophisticated tricks.

• Is sqrt(m1*pi/(R2-R1))^2 really wanted?
• In the calculation of Hvp_0_1_0 you assume that n1 has the last value from the last FOR loop, which uses n1 as counter. Better define n1 and l1 explicitly outside the loops.
• Is "3.142" a crude approximation of pi? If so, don't do it. Never loose accuracy, when there is no benefit.

[zamri]

great, i'll try what u suggested, i'm sorry about the program looks ugly, i'm really new to this, i just need to get some result done for study..i've looked everywhere for other sample program that calculate pressure inside a cavity but n/a. The bessel coefficient was solved manually and using maple.

anyway, the term sqrt(m1*pi/(R2-R1))^2.. is necessary..i don't have idea of how to integrate it to the bessel loop.

[zamri]

dear Jan,

thanks so much for your advise, i was able to run it in 5 minutes. I appreciate it so much. Sorry for being like 5 yr old on this.

[Jan]

I've replaced some constants in the loops and abbreviated whatever I found logical. Of course I cannot debug the results, so I leave this up to you. Some details look strange, like the mentioned sqrt(x)^2.

With less redundancies, the code runs in 100 seconds. Perhaps the code is even less readable as before, sorry. When this code is used for serious science, I would suggest to rewrite it from scratch after mathematical simplifications.

clc
%Torus Dimension
R2=0.275;
R1=0.21;
W=0.205;
ht=0.015;
Theta=2*pi;
%Young's Modulus
E=7.5E+7;
% Poisson's Ratio
GAMA=0.40;
% Plate density
Lo=1350;
Lo0=1.3;
%Plate loss factor
Eta1=0.03;
Eta2=0.03;
c=343;
%excitation force
%F=F0*sin(2*pi*f);
F0=1;
% Excitation force location
thetaf=0;
zf=W/2;
%Plate longitudinal wave velocity
cLt=sqrt(E/Lo/(1-GAMA^2));
% Plate Bending Stiffness
B=E*ht^3/(12*(1-GAMA^2));

dr0=(R2-R1)/20
dtheta0=2*pi/20;
dz0=W/20;

df=100/100;

dr=(R2-R1)/2;
dtheta=2*pi/3;
dz=W/2;

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

HFv=zeros(3,3,51);
HFvv=zeros(3,3,51);
Hvp=zeros(3,3,3,51);
CHvp=zeros(3,3,3,51);
omega_mn1=zeros(5,5);
omega_mn=zeros(5,5);
p_omega=zeros(3,3,3,51);
mean_p2=zeros(3,3,3,51);
omega_n=zeros(5,5,5);
Freq=zeros(51,1);

tmp6 = (1-GAMA^2);

h = waitbar(0,'Please wait...');

for num2 = 200:250
   disp(num2);

     f=1e-10+(num2-1)*df;
     omega=2*3.1415928*f;
     omega2 = omega ^ 2;
     Freq(num2,1)=f;

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

     coef_Hvp=0.0005i*(Theta/2)*(W/2)*omega*c^2*Lo0;
     coef_HFv=2i*omega*F0/(Lo*ht*pi*W);

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

     coef_1_0_0=0.105i*Theta*W*omega*c^2*Lo0;
     coef_0_1_0=0.041i*(Theta/2)*W*omega*c^2*Lo0;
     coef_0_0_1=0.0002i*Theta*(W/2)*omega*c^2*Lo0;
     coef_1_1_0=2.3i*(Theta/2)*W*omega*c^2*Lo0;
     coef_0_1_1=0.0002i*(Theta/2)*(W/2)*omega*c^2*Lo0;
     coef_1_0_1=0.0006i*Theta*(W/2)*omega*c^2*Lo0;

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

     EEc=1e-10;
     EEp=1e-10;

     for mmm2=1:3
        theta=1e-10+(mmm2-1)*dtheta;

        for nnn2=1:3
           z=1e-10+(nnn2-1)*dz;

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

           HFFVV=1e-10;

           gama=W*sqrt(ht)/((R2)^(3/2));
           q2 = (1/sqrt(tmp6*Lo/E));
           for m=1:5
              q1 = sin(m*pi*z/W);
              q3 = sin(m*pi*zf/W);
              for n=1:5
                 omega_mn1(mmm2,nnn2)=sqrt(tmp6*((1+(m*3.142*R2/(n*W))^2)^-2)*(m*3.142/(n*gama)^4)+ ...
                    (1+(m*3.142*R2/(n*W))^2)^2*(n^4)/12);
                 omega_mn(mmm2,nnn2)=omega_mn1(mmm2,nnn2)*(ht)/R2^2*q2;
                 HFFVV=HFFVV+q3*sin(n*pi*thetaf/Theta)*q1*sin(n*pi*theta/Theta)/ ...
                    (omega_mn(m,n)^2-omega2+1i*Eta1*omega_mn(m,n)^2);
              end
           end

           HFvv(mmm2,nnn2,num2)=coef_HFv*HFFVV;

           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           q3 = (ht)/R2^2*(1/sqrt((tmp6)*Lo/E));
           gama=W*sqrt(ht)/((R2)^(3/2));
           q5 = sin((1:5)*pi*thetaf/Theta) .* sin((1:5)*pi*theta0/Theta);

           for lll2=1:3
              r=1e-10+(lll2-1)*dr;

              tmp1 = 7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r);
              tmp2 = -22*besselj(0,96.7*r)-18.2*bessely(0,96.7*r);
              tmp3 = -3*besselj(0,145*r)-17.5*bessely(0,145*r);
              tmp4 = -0.8*besselj(0,193*r)-1.6*bessely(0,193*r);
              tmp5 = -1.8*besselj(0,242*r)-0.37*bessely(0,242*r);

              pw=1e-10;

              for nn= 1:21
                 theta0=1e-10+(nn-1)*dtheta0;

                 for ll=1:21
                    z0=1e-10+(ll-1)*dz0;

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

                    HHFv=1e-10;
                    q12 = (1:5) .^ 4;
                    q13 = ((1:5) * gama) .^ 4;
                    for m=1:5
                       q1 = sin(m*pi*zf/W);
                       q2 = q1 * sin(m*pi*z0/W);
                       q11 = 1 / q6^2;
                       for n=1:5
                          q6 = 1+(m*3.142*R2/(n*W))^2;
                          omega_mn1(mmm2,nnn2) = sqrt(tmp6 * q11 * ...
                             (m*3.142/q13(n)) + q6^2*q12(n)/12);
                          omega_mn(mmm2,nnn2)  = omega_mn1(mmm2,nnn2)*q3;
                          q4 = omega_mn(m,n)^2;
                          HHFv=HHFv + q2 * q5(n) / (q4-omega2+1i*Eta1*q4);
                       end
                    end

                    HFv(mmm2,nnn2,num2)=coef_HFv*HHFv;

                    HHvp=1e-10;
                    HHvp1=1e-10;
                    HHvp2=1e-10;
                    HHvp3=1e-10;
                    HHvp4=1e-10;
                    HHvp5=1e-10;

                    q7 = ((1:5) * pi / W) .^ 2;
                    for m1=1:5
                       q1 = abs(m1*pi/(R2-R1));
                       for n1=1:5
                          q2 = c * q1 + (2*n1/(R2+R1))^2;
                          for l1=1:5
                             omega_n(m1,n1,l1) = q2 + q7(l1);
                          end
                       end
                    end

                    omega_n2 = omega_n .^ 2;
                    omega_x  = omega_n2 - omega2+1i*Eta2*omega_n2;

                    q7 = cos((1:5)*pi*z0/W) .* cos((1:5)*pi*zf/W);
                    q8 = (2/(R2+R1))*theta0;
                    q9 = (2/(R2+R1))*thetaf;
                    for n1=1:5
                       q1 = cos(n1*q8)* cos(n1*q9);
                       for l1=1:5
                          q10   = q1 * q7(l1) / omega_x(m1,n1,l1);
                          HHvp1 = HHvp1 + tmp1 * q10;
                          HHvp2 = HHvp2 + tmp2 * q10;
                          HHvp3 = HHvp3 + tmp3 * q10;
                          HHvp4 = HHvp4 + tmp4 * q10;
                          HHvp5 = HHvp5 + tmp5 * q10;
                       end
                    end

                    HHvp=HHvp1+HHvp2+HHvp3+HHvp4+HHvp5;

                    Hvp(mmm2,nnn2,lll2,num2)= HHvp;

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    n1_ = 5;
                    l1_ = 5;

                    Hvp_1_0_0=tmp5/((c^2*(pi/(R2-R1))^2*(1+1i*Eta2))-omega2);
                    Hvp_0_1_0=(cos(n1_*q8)*cos(n1_*q9))/((c^2*(2/(R2+R1))^2*(1+1i*Eta2))-omega2);
                    Hvp_0_0_1=(cos(l1_*pi*z0/W)*cos(l1_*pi*zf/W))/((c^2*(pi/W)^2*(1+1i*Eta2))-omega2);
                    Hvp_1_1_0=((tmp1)*(cos(n1_*q8)*cos(n1_*q9)))/((c^2*((pi/(R2-R1))^2+(2/(R2+R1))^2)*(1+1i*Eta2))-omega2);
                    Hvp_0_1_1=(cos(n1_*q8)*cos(n1_*q9)*((cos(l1_*pi*z0/W)*cos(l1_*pi*zf/W))))/((c^2*((2/(R2+R1))^2+(pi/W)^2)*(1+1i*Eta2))-omega2);
                    Hvp_1_0_1=((tmp1)*(cos(l1_*pi*z0/W)*cos(l1_*pi*zf/W)))/((c^2*((pi/(R2-R1))^2+(pi/W)^2)*(1+1i*Eta2))-omega2);

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

                    CHvp(mmm2,nnn2,lll2,num2)=coef_Hvp*HHvp+coef_1_0_0*Hvp_1_0_0+coef_0_1_0*Hvp_0_1_0+coef_0_0_1*Hvp_0_0_1+coef_1_1_0*Hvp_1_1_0+coef_0_1_1*Hvp_0_1_1+coef_1_0_1*Hvp_1_0_1;

                    IntFunction=HFv(mmm2,nnn2,num2)*CHvp(mmm2,nnn2,lll2,num2);
                    pw=pw+dtheta0*dz0*IntFunction;

                 end
              end

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

              p_omega(mmm2,nnn2,lll2,num2)=pw;
              mean_p2(mmm2,nnn2,lll2,num2)=abs(p_omega(mmm2,nnn2,lll2,num2)*conj(p_omega(mmm2,nnn2,lll2,num2))/2);

           end
        end
     end

     waitbar(num2/250,h)

end

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

m_p2_CLT1_1_1=zeros(51,1);
m_p2_CLT1_2_1=zeros(51,1);
m_p2_CLT1_3_1=zeros(51,1);
m_p2_CLT2_1_1=zeros(51,1);
m_p2_CLT2_2_1=zeros(51,1);
m_p2_CLT2_3_1=zeros(51,1);
m_p2_CLT3_1_1=zeros(51,1);
m_p2_CLT3_2_1=zeros(51,1);
m_p2_CLT3_3_1=zeros(51,1);

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

for m2=200:250

     m_p2_CLT1_1_1(m2,1)=10*log10(mean_p2(1,1,1,m2));
     m_p2_CLT1_2_1(m2,1)=10*log10(mean_p2(1,2,1,m2));
     m_p2_CLT1_3_1(m2,1)=10*log10(mean_p2(1,3,1,m2));

     m_p2_CLT2_1_1(m2,1)=10*log10(mean_p2(2,1,1,m2));
     m_p2_CLT2_2_1(m2,1)=10*log10(mean_p2(2,2,1,m2));
     m_p2_CLT2_3_1(m2,1)=10*log10(mean_p2(2,3,1,m2));

     m_p2_CLT3_1_1(m2,1)=10*log10(mean_p2(3,1,1,m2));
     m_p2_CLT3_2_1(m2,1)=10*log10(mean_p2(3,2,1,m2));
     m_p2_CLT3_3_1(m2,1)=10*log10(mean_p2(3,3,1,m2));

end

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

figure(1)
plot(Freq(:,1), m_p2_CLT1_1_1(:,1), Freq(:,1), m_p2_CLT1_2_1(:,1), Freq(:,1), m_p2_CLT1_3_1(:,1), 'LineWidth',2);
grid on

figure(2)
plot(Freq(:,1), m_p2_CLT2_1_1(:,1), Freq(:,1), m_p2_CLT2_2_1(:,1), Freq(:,1), m_p2_CLT2_3_1(:,1), 'LineWidth',2);
grid on

figure(3)
plot(Freq(:,1), m_p2_CLT3_1_1(:,1), Freq(:,1), m_p2_CLT3_2_1(:,1), Freq(:,1), m_p2_CLT3_3_1(:,1), 'LineWidth',2);
grid on

Answer 2

[Jan]

Look at this piece of code:

for n1=1:5
  for l1=1:5
    HHvp4= HHvp4 + d*n1*l1
  end
end

This can be abbreviated to:

HHvp4 = HHvp4 + 225 * d;

A pre-allocation is strongly recommended, when you fill an array like omega_n iteratively. But your example runs in 0.000032 seconds on my computer, such that I do not think that this is a bottleneck - except if you call this millions of times, and if so, please explain this. Anyhow, I would create this array like this:

omega_n = zeros(5, 5, 5);
omega_n(:, :, 1) = (1:5)' * (1:5);
for ii = 2:5
  omega_n(:, :, ii) = ii * omega_n(:, :, 1);
end

Puh, looks more like Matlab, but it takes 0.000020 also. I do not assume that 0.01 thousandth seconds have been worth to improve this part. But who cares about minutes of programming time, when the execution time is 0.1 seconds shorter?! This is the fastest method I've found under R2009a/64 on a Windows7/Core2Duo:

omega_n = zeros(5, 5, 5);
for l1=1:5
  for n1=1:5
    for m1 = 1:5
       omega_n(m1, n1, l1) = m1 * n1 * l1;
    end
  end
end

0.0000089 seconds. The order of loops has been reversed, such that neighboring elements are written. But as long as 125*8 bytes match into the processor cache, this is not essential. Using a temporary variable c = n1 * l1 in the inner loop is more efficient in theory, but the improvement is not significant.

[zamri]

Thanks for your answer Jan, i am new to matlab, same to programming. I paste all the codes here, just in case u/others can point out what's the problem that it takes 15 hrs to complete (for 50 hz range). If i want t get 0-500 Hz range, i'll probably take 15*10=150 hrs.

clear all
close all
clc
%Torus Dimension
R2=0.275;
R1=0.21;
W=0.205;
ht=0.015;
Theta=2*pi;
%Young's Modulus
E=7.5E+7;
% Poisson's Ratio
GAMA=0.40;
% Plate density
Lo=1350;
Lo0=1.3;
%Plate loss factor
Eta1=0.03;
Eta2=0.03;
c=343;
%excitation force 
%F=F0*sin(2*pi*f);
F0=1;
% Excitation force location
thetaf=0;
zf=W/2;
%Plate longitudinal wave velocity
cLt=sqrt(E/Lo/(1-GAMA^2));
% Plate Bending Stiffness
B=E*ht^3/(12*(1-GAMA^2));

dr0=(R2-R1)/20
dtheta0=2*pi/20;
dz0=W/20;

df=100/100;

dr=(R2-R1)/2;
dtheta=2*pi/3; 
dz=W/2;

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

       HFv=zeros(3,3,51);
       HFvv=zeros(3,3,51);
       Hvp=zeros(3,3,3,51);
       CHvp=zeros(3,3,3,51);
       omega_mn1=zeros(5,5);
       omega_mn=zeros(5,5);
       p_omega=zeros(3,3,3,51);
       mean_p2=zeros(3,3,3,51);
       omega_n=zeros(5,5,5);
       Freq=zeros(51,1); 

       h = waitbar(0,'Please wait...'); 

for num2 = 200:250

f=1e-10+(num2-1)*df; 
omega=2*3.1415928*f;
Freq(num2,1)=f;

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

coef_Hvp=0.0005*(Theta/2)*(W/2)*i*omega*c^2*Lo0;
coef_HFv=2*i*omega*F0/(Lo*ht*pi*W);

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

coef_1_0_0=0.105*Theta*W*i*omega*c^2*Lo0;
coef_0_1_0=0.041*(Theta/2)*W*i*omega*c^2*Lo0;
coef_0_0_1=0.0002*Theta*(W/2)*i*omega*c^2*Lo0;
coef_1_1_0=2.3*(Theta/2)*W*i*omega*c^2*Lo0;
coef_0_1_1=0.0002*(Theta/2)*(W/2)*i*omega*c^2*Lo0;
coef_1_0_1=0.0006*Theta*(W/2)*i*omega*c^2*Lo0;

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

            EEc=1e-10;  
            EEp=1e-10;

           for mmm2=1:3
                   theta=1e-10+(mmm2-1)*dtheta;

           for nnn2=1:3
                   z=1e-10+(nnn2-1)*dz;

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

HFFVV=1e-10;          

    for m=1:5
       for n=1:5
           gama=W*sqrt(ht)/((R2)^(3/2));
           omega_mn1(mmm2,nnn2)=sqrt((1-GAMA^2)*((1+(m*3.142*R2/(n*W))^2)^-2)*(m*3.142/(n*gama)^4)+(1+(m*3.142*R2/(n*W))^2)^2*(n^4)/12);
           omega_mn(mmm2,nnn2)=omega_mn1(mmm2,nnn2)*(ht)/R2^2*(1/sqrt((1-GAMA^2)*Lo/E));
          HFFVV=HFFVV+sin(m*pi*zf/W)*sin(n*pi*thetaf/Theta)*sin(m*pi*z/W)*sin(n*pi*theta/Theta)/(omega_mn(m,n)^2-omega^2+i*Eta1*omega_mn(m,n)^2);
             end
         end

       HFvv(mmm2,nnn2,num2)=coef_HFv*HFFVV;  

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

           for lll2=1:3
                   r=1e-10+(lll2-1)*dr;    

pw=1e-10;                 

for nn= 1:21
theta0=1e-10+(nn-1)*dtheta0;

for ll=1:21
z0=1e-10+(ll-1)*dz0;

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

   HHFv=1e-10;          

    for m=1:5
       for n=1:5  

           gama=W*sqrt(ht)/((R2)^(3/2));
           omega_mn1(mmm2,nnn2)=sqrt((1-GAMA^2)*((1+(m*3.142*R2/(n*W))^2)^-2)*(m*3.142/(n*gama)^4)+(1+(m*3.142*R2/(n*W))^2)^2*(n^4)/12);
           omega_mn(mmm2,nnn2)=omega_mn1(mmm2,nnn2)*(ht)/R2^2*(1/sqrt((1-GAMA^2)*Lo/E));
          HHFv=HHFv+sin(m*pi*zf/W)*sin(n*pi*thetaf/Theta)*sin(m*pi*z0/W)*sin(n*pi*theta0/Theta)/(omega_mn(m,n)^2-omega^2+i*Eta1*omega_mn(m,n)^2);
       end
    end

       HFv(mmm2,nnn2,num2)=coef_HFv*HHFv;  

   HHvp=1e-10;
   HHvp1=1e-10;
   HHvp2=1e-10;
   HHvp3=1e-10;
   HHvp4=1e-10;
   HHvp5=1e-10;

       for m1=1:5
           for n1=1:5
               for l1=1:5
                omega_n(m1,n1,l1)=c*sqrt(m1*pi/(R2-R1))^2+(2*n1/(R2+R1))^2+(l1*pi/W)^2;
               end
           end
       end

       for n1=1:5
         for l1=1:5
            HHvp1= HHvp1+(7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r))*cos(n1*(2/(R2+R1))*theta0)*cos(l1*pi*z0/W)*cos(n1*(2/(R2+R1))*thetaf)*cos(l1*pi*zf/W)/(omega_n(m1,n1,l1)^2-omega^2+i*Eta2*omega_n(m1,n1,l1)^2);
         end
       end

       for n1=1:5
         for l1=1:5
            HHvp2= HHvp2+(-22*besselj(0,96.7*r)-18.2*bessely(0,96.7*r))*cos(n1*(2/(R2+R1))*theta0)*cos(l1*pi*z0/W)*cos(n1*(2/(R2+R1))*thetaf)*cos(l1*pi*zf/W)/(omega_n(m1,n1,l1)^2-omega^2+i*Eta2*omega_n(m1,n1,l1)^2);
         end
       end

       for n1=1:5
         for l1=1:5
            HHvp3= HHvp3+(-3*besselj(0,145*r)-17.5*bessely(0,145*r))*cos(n1*(2/(R2+R1))*theta0)*cos(l1*pi*z0/W)*cos(n1*(2/(R2+R1))*thetaf)*cos(l1*pi*zf/W)/(omega_n(m1,n1,l1)^2-omega^2+i*Eta2*omega_n(m1,n1,l1)^2);
         end
       end

       for n1=1:5
         for l1=1:5
            HHvp4= HHvp4+(-0.8*besselj(0,193*r)-1.6*bessely(0,193*r))*cos(n1*(2/(R2+R1))*theta0)*cos(l1*pi*z0/W)*cos(n1*(2/(R2+R1))*thetaf)*cos(l1*pi*zf/W)/(omega_n(m1,n1,l1)^2-omega^2+i*Eta2*omega_n(m1,n1,l1)^2);
         end
       end        

       for n1=1:5
         for l1=1:5
            HHvp5= HHvp5+(-1.8*besselj(0,242*r)-0.37*bessely(0,242*r))*cos(n1*(2/(R2+R1))*theta0)*cos(l1*pi*z0/W)*cos(n1*(2/(R2+R1))*thetaf)*cos(l1*pi*zf/W)/(omega_n(m1,n1,l1)^2-omega^2+i*Eta2*omega_n(m1,n1,l1)^2);
         end
       end 

       HHvp=HHvp1+HHvp2+HHvp3+HHvp4+HHvp5;

       Hvp(mmm2,nnn2,lll2,num2)= HHvp;    

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

                   Hvp_1_0_0=(7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r))/((c^2*(pi/(R2-R1))^2*(1+i*Eta2))-omega^2);
                   Hvp_0_1_0=(cos(n1*(2/(R2+R1))*theta0)*cos(n1*(2/(R2+R1))*thetaf))/((c^2*(2/(R2+R1))^2*(1+i*Eta2))-omega^2);
                   Hvp_0_0_1=(cos(l1*pi*z0/W)*cos(l1*pi*zf/W))/((c^2*(pi/W)^2*(1+i*Eta2))-omega^2);
                   Hvp_1_1_0=((7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r))*(cos(n1*(2/(R2+R1))*theta0)*cos(n1*(2/(R2+R1))*thetaf)))/((c^2*((pi/(R2-R1))^2+(2/(R2+R1))^2)*(1+i*Eta2))-omega^2);
                   Hvp_0_1_1=(cos(n1*(2/(R2+R1))*theta0)*cos(n1*(2/(R2+R1))*thetaf)*((cos(l1*pi*z0/W)*cos(l1*pi*zf/W))))/((c^2*((2/(R2+R1))^2+(pi/W)^2)*(1+i*Eta2))-omega^2);
                   Hvp_1_0_1=((7.7*besselj(0,48.3*r)-0.32*bessely(0,48.3*r))*(cos(l1*pi*z0/W)*cos(l1*pi*zf/W)))/((c^2*((pi/(R2-R1))^2+(pi/W)^2)*(1+i*Eta2))-omega^2);

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

                   CHvp(mmm2,nnn2,lll2,num2)=coef_Hvp*HHvp+coef_1_0_0*Hvp_1_0_0+coef_0_1_0*Hvp_0_1_0+coef_0_0_1*Hvp_0_0_1+coef_1_1_0*Hvp_1_1_0+coef_0_1_1*Hvp_0_1_1+coef_1_0_1*Hvp_1_0_1; 

IntFunction=HFv(mmm2,nnn2,num2)*CHvp(mmm2,nnn2,lll2,num2);
pw=pw+dtheta0*dz0*IntFunction;

end
end

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

p_omega(mmm2,nnn2,lll2,num2)=pw;
mean_p2(mmm2,nnn2,lll2,num2)=abs(p_omega(mmm2,nnn2,lll2,num2)*conj(p_omega(mmm2,nnn2,lll2,num2))/2);

end
end
end

waitbar(num2/250,h)

end

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

m_p2_CLT1_1_1=zeros(51,1);
m_p2_CLT1_2_1=zeros(51,1);
m_p2_CLT1_3_1=zeros(51,1);
m_p2_CLT2_1_1=zeros(51,1);
m_p2_CLT2_2_1=zeros(51,1);
m_p2_CLT2_3_1=zeros(51,1);
m_p2_CLT3_1_1=zeros(51,1);
m_p2_CLT3_2_1=zeros(51,1);
m_p2_CLT3_3_1=zeros(51,1);

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

for m2=200:250

      m_p2_CLT1_1_1(m2,1)=10*log10(mean_p2(1,1,1,m2))
      m_p2_CLT1_2_1(m2,1)=10*log10(mean_p2(1,2,1,m2))
      m_p2_CLT1_3_1(m2,1)=10*log10(mean_p2(1,3,1,m2))

      m_p2_CLT2_1_1(m2,1)=10*log10(mean_p2(2,1,1,m2))
      m_p2_CLT2_2_1(m2,1)=10*log10(mean_p2(2,2,1,m2))
      m_p2_CLT2_3_1(m2,1)=10*log10(mean_p2(2,3,1,m2))

      m_p2_CLT3_1_1(m2,1)=10*log10(mean_p2(3,1,1,m2))
      m_p2_CLT3_2_1(m2,1)=10*log10(mean_p2(3,2,1,m2))
      m_p2_CLT3_3_1(m2,1)=10*log10(mean_p2(3,3,1,m2))

    end

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

figure(1)
plot(Freq(:,1), m_p2_CLT1_1_1(:,1), Freq(:,1), m_p2_CLT1_2_1(:,1), Freq(:,1), m_p2_CLT1_3_1(:,1), 'LineWidth',2);
grid on

figure(2)
plot(Freq(:,1), m_p2_CLT2_1_1(:,1), Freq(:,1), m_p2_CLT2_2_1(:,1), Freq(:,1), m_p2_CLT2_3_1(:,1), 'LineWidth',2);
grid on

figure(3)
plot(Freq(:,1), m_p2_CLT3_1_1(:,1), Freq(:,1), m_p2_CLT3_2_1(:,1), Freq(:,1), m_p2_CLT3_3_1(:,1), 'LineWidth',2);
grid on

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[Azzi Abdelmalek]

Zamri, there are some missing data, which allows to test your code (like F0). Also there are parts of your code which are not well formated.

[Jan]

@Zamri: I've formatted your code again. Please follow the linked instructions. Thanks.