How to code following equation?

My code for the above equation is given below.

But looks like it don't work well. Could anyone have a look at whether I implemented the above formula well?

clear all;clc;close all;
syms j l k p q s K
m = 3;
z=linspace(0.0001,5000);
w0=0.05;
Omega = 30;
k_w = 4.05*10^6;
gama=(1./(w0.^2))+((1i.*k_w)./(2.*z));
gama_star = subs(gama, 1i, -1i);
con = 1./(16.*pi.*pi.*z.*z);

total = 0;
for j = 0:m
    J = nchoosek(m,j).*((1i).^(m-j));
    sum = 0;
    for l = 0:m
        L = nchoosek(m,l).*((1i).^(m-l));
        for k = 0:1/2:l/2
            KK = (((-1).^k).*factorial(l))./(gamma(k+1).*factorial(l-(2.*k)));
            for q= 0:(l-2*k)
                 Q = nchoosek((l-(2.*k)),q);
                 for p= 0:1/2:(m-l)/2
                     P = (((2.*1i)./(sqrt(gama))).^(m-(2.*p)-(2.*k)));
                     for s = 0:(m-l-2*p)
                        a = (hermiteH(m-j+s,((Omega.*1i)./sqrt(gama)))).*(hermiteH(j+q,((Omega.*1i)./sqrt(gama))));
                        b = (exp(Omega).*((Omega).^(m-s-q-(2.*p)-(2.*k)))) - (exp(-Omega).*((-Omega).^(m-s-q-(2.*p)-(2.*k))));
                        S = nchoosek((m-l-(2.*p)),s).*a.*b.*(1./((2.*1i.*sqrt(gama_star)).^(m+s+q)));
                     end
                 end
            end
        end
    end
    sum = sum + S.*P.*Q.*KK.*L;
end
total = sum.*con;
F = vpa(real(total))

F = 

[Torsten]

But looks like it don't work well.

Did you check using pencil and paper or how do you know ?

[Athira T Das]

I am trying to plot F w.r.t z for different values of m. And the graph is not good enough.

[Torsten]

"not good enough" means "wrong" ?

[Athira T Das]

@TorstenYes

Answer 1

[Image Analyst]

Many times with complicated equations it's difficult to code and often it's because the parentheses were misplaced. I recommend you break it up into smaller terms, like term1, term2, term3, sum1, sum2, sum3, etc., and then combine them.

[Athira T Das]

How can I split these terms, since running index of each summation is related to the previous summation?

[Image Analyst]

Start with the innermost look and compute the sum for that, like

sumOverS = 0
for s = 1 : 0:(m-l-2*p)
    sumOverS = whatever.....
end

Then after the s loop but before the end of the p loop, use the sumOverS that you just got from the s loop.

[Athira T Das]

I tried this way, but still it is wrong

F = 0;
c=1./(16.*pi.^2.*z.^2)
for l = 0:m
  L = ((1i).^(m-l)).*nchoosek(m,l);
  sum_j = 0;
  for j = 0:m
      J= ((1i).^(m-j)).*nchoosek(m,j);  
      sum1 = 0.0;
      for k = 0:1/2:l/2
            faktor1_K = (((-1).^k).*factorial(l))./(gamma(k+1).*factorial(l-(2.*k)));
            for q = 0:l-2*k
                faktor1_Q =  nchoosek((l-(2.*k)),q);
                for p = 0:1/2:(m-l)/2
                    faktor1_P = (((2.*1i)./(sqrt(Gamma))).^(m-(2.*p)-(2.*k)));
                    sum_inner = 0;
                    for s = 0:m-l-2*p
                        E1 = ((exp(Omega).*(Omega.^(m-s-q-2.*p-2.*k))) - (exp(-Omega).*(-Omega.^(m-s-q-2.*p-2.*k)))).*hermiteH(j+q,((1i.*Omega)./(sqrt(Gamma)))).*hermiteH(m-j+s,((1i.*Omega)./(sqrt(Gamma))));
                        sum_inner = sum_inner + nchoosek((m-l-(2.*p)),s).*(1./(((2.*1i.*sqrt(Gamma_star))).^(m+q+s))).*E1;
                    end
                end
            end
            sum1 = sum1 + faktor1_K.*faktor1_P.*faktor1_Q.*sum_inner;
      end
      sum_j= sum_j + J.*(sum1);
  end
 F = F + (L.*sum_j);
end  
F0 = c.*F  

[Image Analyst]

That's because you have not defined z:

Unrecognized function or variable 'z'.

Error in test8 (line 2)

c=1./(16.*pi.^2.*z.^2)

[Torsten]

m = 3;
z=linspace(0.0001,5000);
w0=0.05;
Omega = 30;
k_w = 4.05*10^6;
Gamma=(1./(w0.^2))+((1i.*k_w)./(2.*z));
Gamma_star = (1./(w0.^2))+((-1i.*k_w)./(2.*z));
F = 0;
c=1./(16.*pi.^2.*z.^2);
for l = 0:m
    L = ((1i).^(m-l)).*nchoosek(m,l);
    sum_j = 0;
    for j = 0:m
        J= ((1i).^(m-j)).*nchoosek(m,j);  
        sum1 = 0.0;
        for k = 0:1/2:l/2
            faktor1_K = (((-1).^k).*factorial(l))./(gamma(k+1).*factorial(l-(2.*k)));
            for q = 0:l-2*k
                faktor1_Q =  nchoosek((l-(2.*k)),q);
                for p = 0:1/2:(m-l)/2
                    faktor1_P = (((2.*1i)./(sqrt(Gamma))).^(m-(2.*p)-(2.*k)));
                    sum_inner = 0;
                    for s = 0:m-l-2*p
                        E1 = ((exp(Omega).*(Omega.^(m-s-q-2.*p-2.*k))) - (exp(-Omega).*(-Omega.^(m-s-q-2.*p-2.*k)))).*hermiteH(j+q,((1i.*Omega)./(sqrt(Gamma)))).*hermiteH(m-j+s,((1i.*Omega)./(sqrt(Gamma))));
                        sum_inner = sum_inner + nchoosek((m-l-(2.*p)),s).*(1./(((2.*1i.*sqrt(Gamma_star))).^(m+q+s))).*E1;
                    end
                end
            end
            sum1 = sum1 + faktor1_K.*faktor1_P.*faktor1_Q.*sum_inner;
        end
        sum_j= sum_j + J.*(sum1);
    end
    F = F + (L.*sum_j);
end  
F0 = c.*F  

F0 = 
   0.0000 + 0.0000i   0.3312 + 0.3312i   0.6623 + 0.6622i   0.9928 + 0.9927i   1.3225 + 1.3224i   1.6512 + 1.6509i   1.9784 + 1.9780i   2.3040 + 2.3034i   2.6276 + 2.6267i   2.9489 + 2.9478i   3.2678 + 3.2663i   3.5839 + 3.5820i   3.8969 + 3.8946i   4.2066 + 4.2039i   4.5129 + 4.5096i   4.8153 + 4.8115i   5.1138 + 5.1093i   5.4081 + 5.4029i   5.6979 + 5.6920i   5.9831 + 5.9764i   6.2635 + 6.2560i   6.5390 + 6.5306i   6.8093 + 6.8000i   7.0742 + 7.0640i   7.3338 + 7.3226i   7.5877 + 7.5755i   7.8359 + 7.8226i   8.0783 + 8.0639i   8.3148 + 8.2993i   8.5452 + 8.5286i   8.7696 + 8.7517i   8.9877 + 8.9687i   9.1997 + 9.1794i   9.4054 + 9.3839i   9.6047 + 9.5819i   9.7977 + 9.7737i   9.9844 + 9.9590i  10.1646 +10.1380i  10.3385 +10.3106i  10.5060 +10.4768i  10.6672 +10.6367i  10.8220 +10.7902i  10.9706 +10.9374i  11.1129 +11.0784i  11.2490 +11.2132i  11.3789 +11.3419i  11.5028 +11.4645i  11.6206 +11.5810i  11.7325 +11.6916i  11.8385 +11.7964i  11.9387 +11.8954i  12.0331 +11.9886i  12.1220 +12.0763i  12.2053 +12.1585i  12.2832 +12.2353i  12.3558 +12.3067i  12.4232 +12.3730i  12.4854 +12.4342i  12.5427 +12.4903i  12.5950 +12.5417i  12.6425 +12.5882i  12.6854 +12.6301i  12.7237 +12.6675i  12.7575 +12.7004i  12.7870 +12.7291i  12.8123 +12.7535i  12.8335 +12.7739i  12.8507 +12.7903i  12.8641 +12.8029i  12.8736 +12.8117i  12.8795 +12.8170i  12.8819 +12.8187i  12.8809 +12.8170i  12.8766 +12.8121i  12.8690 +12.8040i  12.8584 +12.7928i  12.8448 +12.7787i  12.8283 +12.7617i  12.8091 +12.7420i  12.7871 +12.7196i  12.7626 +12.6947i  12.7357 +12.6674i  12.7063 +12.6377i  12.6747 +12.6057i  12.6409 +12.5716i  12.6050 +12.5354i  12.5670 +12.4972i  12.5272 +12.4572i  12.4855 +12.4153i  12.4421 +12.3717i  12.3970 +12.3264i  12.3503 +12.2796i  12.3021 +12.2313i  12.2524 +12.1815i  12.2014 +12.1304i  12.1491 +12.0780i  12.0956 +12.0245i  12.0408 +11.9697i  11.9850 +11.9139i  11.9282 +11.8571i