How to execute double for loop?

I have an equation

.

This particular equation is for an z value.

But I have an array of z and I need to evaluate get C as an array. But currently I have an error in the final step. Can anyone help me?

clc;close all;clear all;
syms r ph phi 
mu=1 
w0=0.01;
lambda=532*10^-9;
R=0.035;
omega=.1;  
zr=pi.*w0.^2./lambda;
k=2*pi/lambda;
z=[1.037	2.014	4.91	6.638	8.097	10.179	12.804	14.491	16.594	18.136	20.468	22.842	26.091	28.695	31.132	32.486	34.735	36.173	38.068	40.256	42.089	44.672	46.358	48.066	50.045	52.816	54.523	56.191	58.919	61.127	62.085	64.106	66.127	67.856	70.084	72.272	74.189	76.063	78.064	80.522	82.876	85.856	86.147	88.189	90.606	92.606	94.584	96.23	98.336	99.794  100									];
T1=[0.261864786659278	0.0279450575402681	0.0167207441613120	0.0585405683053739	0.0467736536041428	0.0721143056065063	0.0204055105411778	0.0414970149223034	0.0370574720556092	0.0323354183816648	0.0557187803983150	0.0209702231071337	0.0523815841575279	0.0478073675740971	0.0452486568524416	0.0454216694742521	0.0455017342270441	0.0433366921457154	0.0416287144440869	0.0325563861379645	0.0453942497662185	0.0308978986984992	0.0388337989311703	0.0395470194181856	0.0412953384586318	0.0359798229891002	0.0343224345520140	0.0388190479592265	0.0357132945620167	0.0338578554357400	0.0294508197964149	0.0346296299864533	0.0328074293196396	0.0327552539365064	0.0332008509577915	0.0266013611176686	0.0351599194601709	0.0325287304943849	0.0320666044429463	0.0321364098469872	0.0321092261908771	0.0104350167431096	0.0155659791192790	0.00845361754599249	0.0304932945576984	0.0286130247196421	0.0302806888132458	0.0290252186616402	0.0323193263219266	0.0321277263184264	0.0320992548815438]  ;
L=3;
N=2*L+1;
m0=1;
m=1 ;
M=[-2 -1 0 1 2];
M0=linspace(-L,L,N);

%to calculate P(m/m0)
wz=w0.*sqrt(1+(z./zr).^2) ;
Con1 = (1i./(2.*lambda.*z)).*exp((-1i.*k.*r.*r)./(2.*z));
for l2=0:m 
    E =-1./(wz.^2) + (1i.*k)./(2.*z) ;
    E2 =  exp((((omega/2)-(r.*cos(ph))).^2)./E).* exp((((omega/2)-(r.*sin(ph))).^2)./E).*hermiteH(l2,1i*(omega./2 - r.*cos(ph))./sqrt(E)).*hermiteH(m-l2,1i*(omega./2 - r.*sin(ph))./sqrt(E))  +  exp((((-omega/2)-(r.*cos(ph))).^2)./E).* exp((((-omega/2)-(r.*sin(ph))).^2)./E).*hermiteH(l2,1i*(-omega./2 - r.*cos(ph))./sqrt(E)).*hermiteH(m-l2,1i*(-omega./2 - r.*sin(ph))./sqrt(E));
    I2 =  ((1i).^(m-l2)).*nchoosek(m,l2).*E2;
end
Is = (Con1.*conj(Con1)).*I2.*conj(I2);


%to calculate Sum over P(m/m0)
parfor j= 1:length(z)
    for m0=1:length(M0)
        IQ1 =Is(j).*r.*exp((-2.*r.*r)./(T1(j)^2)).*exp((-1i.*(m-m0).*(phi-ph)) + (cos(ph-phi).*2.*r.*r)./(T1(j).^2));
        fun = matlabFunction(IQ1,'Vars',[r,ph,phi]);
        result_2 =  (1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)) ;
        Result2(j) = result_2
    end
end
Result2

parfor j= 1:length(z)  
    for m=1:length(M)
        for m0=1:length(M0)
             
            Is = (Con1.*conj(Con1)).*I2.*conj(I2);
            IQ1 =Is(j).*r.*exp((-2.*r.*r)./(T1(j)^2)).*exp((-1i.*(m-m0).*(phi-ph)) + (cos(ph-phi).*2.*r.*r)./(T1(j).^2));
            fun = matlabFunction(IQ1,'Vars',[r,ph,phi]);
            result =log2(N) + (1/N).*(  (1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)) ).*log2(((1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)))./(Result2))  
            C(j)= result         
        end   
    end    
end 
C

Answer 1

[Walter Roberson]

You have nested parfor loops, but each time you are overwriting the (j) entry

            result =log2(N) + (1/N).*(  (1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)) ).*log2(((1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)))./(Result2))  
            

You are dividing by the vector Result2 so you get a vector result.

Answer 2

[Athira T Das]

I changed my code like

fun = matlabFunction(IQ1,'Vars',[r,ph,phi]);
result_3 =  (1/(4.*pi.*pi)).*abs(integral3(@(r,ph,phi)fun(r,ph,phi),0, R/2,0,2*pi,0,2*pi)) ;
result_4 = result_3.*log2(result_3./Result2);
C(j) = log2(N) + (1/N).*result_4

Having error

[Walter Roberson]

result_4 = result_3.*log2(result_3./Result2);

That is a matrix operation involving all of Result2, so it gives a matrix result.

Suggestion: divide by Result2(j)