What is wrong in my script?

clc; clear all; close all;
syms m1 m2 s1 s2 j1 j2 
I = 0; 
lambda = 1060*10^-9;
M=0
z=linspace(0.00001,8000)
wo = 0.02;
C = 10^(-7);
k=2*pi/lambda;
po=(0.545*C^2*k^2*z).^(-(3/5));
b=0.1;
T1=0;
T2=0;
T3=0;
T4=0;
T5=0;
T6=0;
delta= ((1i*k)./(2*z))+(1./(wo.^2))+(1./(po.^2));
delta_star= subs(delta, 1i, -1i);
etta = delta_star - (1./(delta.*po.^4));
X=((k./(2.*z)).^2).*((1./(2.*1i.*sqrt(delta))).^M).*(1./(16.*delta.*etta));
alpha = (1i.*k./(2.*z)).*(1./(delta.*po.^2)-1);
beta_p = (b./(2.*wo)).*(1./(delta.*po.^2)+1);
beta_n = (b./(2.*wo)).*(1./(delta.*po.^2)-1);
A = (alpha.^2./etta)-((k.^2)-(4.*z.^2.*delta));
B1_p = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta);
B1_n = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta);
B2_p = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta));
B2_n = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta));
C1_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
C1_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);

e1=(exp(C1_p).*hermiteH(m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(m2+j2,((-1i.*beta_p)./sqrt(etta)))));
e2=(exp(C1_p).*hermiteH(M-m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(M-m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-m2+j2,((-1i.*beta_p)./sqrt(etta)))));
O = zeros(size(z));

for m1 =0:M
    T2=0;
    for m2=0:M
        T3=0;
        for s1=0:m1/2
            T4=0;
            
            for j1=0:m1-(2*s1)
                T5=0;
                
                for s2=0:(M-m1)/2
                    T6=0;
                    
                    for j2=0:(M-m1-(2*s2))
                                          
                                    
                        T6=T6+(factorial(M-m1-(2.*s2))/(factorial(j2)*factorial(M-m1-(2.*s2)-j2))).*((1./(2.*1i.*sqrt(etta))).^(M-m2+j2)).*((1./(po.^2)).^j2).*((b./2.*wo).^(M-m1-2*s2-j2)).*e2;
                    end
                    T5=T5+(((factorial(M-m1).*(-1).^s2)./(factorial(s2).*factorial(M-m1-2.*s2))).*(((2.*1i)./(delta.^(0.5))).^(M-m1-2.*s2))).*T6;
                end
                 if (m1-(2*s2))<0
                    break
                 end
                
                T4=T4+(((factorial(m1-(2.*s1)))/(factorial(j1)*factorial(m1-(2.*s1)-j1))).*((1./(2.*1i.*sqrt(etta))).^(m2+j1)).*((1./(po.^2)).^j1).*((b./2.*wo).^(m1-2.*s1-j1))).*e1.*T5;
            end
            T3=T3+(((factorial(m1).*(-1).^s1)./(factorial(s1).*factorial(m1-(2.*s1)))).*(((2.*1i)./(sqrt(delta))).^(m1-2.*s1))).*T4;
        end
        T2=T2+nchoosek(M,m2)*((-1i)^(M-m2))*T3;
    end
    T1=T1+nchoosek(M,m1)*((1i)^(M-m1))*T2;
    
end
I0 = -real(X.*T1)

I_numerical = double(vpa(I0))
I_o=abs(I_numerical)

writematrix(I_o,"Normalized_intensity_M=0_b0.1_if.xlsx")
writematrix(z,"Propagation_distance.xlsx")

[Jan]

You forgot to mention, why you assume, that anything is wrong.

[Athira T Das]

Not matching with the numerical results

[Jan]

@Athira T Das: Your comment does now allow to reproduce, what you assume as problem.

[Athira T Das]

@Jan:Evaluated this equation for different M values from 0 to 3.

And plotted the normalized I_o.

But my plotted results is not matching with the actual result.

[Torsten]

Could you include a plot of the expected results ?

And these "expected results" come from a reliable source ?

[Athira T Das]

@TorstenYes it is from reliable source

Answer 1

[Walter Roberson]

Your code defines e1 and e2 in terms of symbolic variables m2 and j2 . Later it does for loops that assign values to m2 and j2. However, assigning a value to m2 and j2 at the MATLAB level does not affect the symbolic variables by the same name. You need to subs() the numeric values in.

format long g
syms m2 s1 s2 j1 j2 
I = 0; 
lambda = 1060*10^-9;
M=0;
z=linspace(0.00001,8000);
wo = 0.02;
C = 10^(-7);
k=2*pi/lambda;
po=(0.545*C^2*k^2*z).^(-(3/5));
b=0.1;
T1=0;
T2=0;
T3=0;
T4=0;
T5=0;
T6=0;
delta= ((1i*k)./(2*z))+(1./(wo.^2))+(1./(po.^2));
delta_star= subs(delta, 1i, -1i);
etta = delta_star - (1./(delta.*po.^4));
X=((k./(2.*z)).^2).*((1./(2.*1i.*sqrt(delta))).^M).*(1./(16.*delta.*etta));
alpha = (1i.*k./(2.*z)).*(1./(delta.*po.^2)-1);
beta_p = (b./(2.*wo)).*(1./(delta.*po.^2)+1);
beta_n = (b./(2.*wo)).*(1./(delta.*po.^2)-1);
A = (alpha.^2./etta)-((k.^2)-(4.*z.^2.*delta));
B1_p = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta);
B1_n = ((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta);
B2_p = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_p)./etta));
B2_n = -(((1i.*k.*b)./(2.*z.*wo.*delta))+((2.*alpha.*beta_n)./etta));
C1_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);
C1_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_p = ((b.^2)./(4.*wo.^2.*delta))+(((beta_n).^2)./etta);
C2_n = ((b.^2)./(4.*wo.^2.*delta))+(((beta_p).^2)./etta);

e1=(exp(C1_p).*hermiteH(m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(m2+j2,((-1i.*beta_p)./sqrt(etta)))));
e2=(exp(C1_p).*hermiteH(M-m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-m2+j2,((-1i.*beta_n)./sqrt(etta))))+((exp(C1_n).*hermiteH(M-m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-m2+j2,((-1i.*beta_p)./sqrt(etta)))));
O = zeros(size(z));

for m1 =0:M
    T2=0;
    for m2=0:M
        T3=0;
        for s1=0:m1/2
            T4=0;
            
            for j1=0:m1-(2*s1)
                T5=0;
                
                for s2=0:(M-m1)/2
                    T6=0;
                    
                    for j2=0:(M-m1-(2*s2))
                                          
                                    
                        T6=T6+(factorial(M-m1-(2.*s2))/(factorial(j2)*factorial(M-m1-(2.*s2)-j2))).*((1./(2.*1i.*sqrt(etta))).^(M-m2+j2)).*((1./(po.^2)).^j2).*((b./2.*wo).^(M-m1-2*s2-j2)).*subs(e2);
                    end
                    T5=T5+(((factorial(M-m1).*(-1).^s2)./(factorial(s2).*factorial(M-m1-2.*s2))).*(((2.*1i)./(delta.^(0.5))).^(M-m1-2.*s2))).*T6;
                end
                 if (m1-(2*s2))<0
                    break
                 end
                
                T4=T4+(((factorial(m1-(2.*s1)))/(factorial(j1)*factorial(m1-(2.*s1)-j1))).*((1./(2.*1i.*sqrt(etta))).^(m2+j1)).*((1./(po.^2)).^j1).*((b./2.*wo).^(m1-2.*s1-j1))).*subs(e1).*T5;
            end
            T3=T3+(((factorial(m1).*(-1).^s1)./(factorial(s1).*factorial(m1-(2.*s1)))).*(((2.*1i)./(sqrt(delta))).^(m1-2.*s1))).*T4;
        end
        T2=T2+nchoosek(M,m2)*((-1i)^(M-m2))*T3;
    end
    T1=T1+nchoosek(M,m1)*((1i)^(M-m1))*T2;
    
end
I0 = -real(X.*T1);

I_numerical = double(vpa(I0));
I_o=abs(I_numerical);

writematrix(I_o,"Normalized_intensity_M=0_b0.1_if.xlsx")
writematrix(z,"Propagation_distance.xlsx")