Error in multiplication of arrays

clc; clear all; close all;
syms m1 m2 s1 s2 j1 j2 
I = 0; 
lambda = 1060*10^-9;
M=2

M = 2

z=linspace(0.0001,10000);
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);



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))
                        e1=(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(m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(m2+j2,((-1i.*beta_n)./sqrt(etta))));
                        e3=(exp(C1_p).*hermiteH(M-m2+j2,((1i.*beta_p)./sqrt(etta))))+(exp(C2_n).*hermiteH(M-m2+j2,((-1i.*beta_n)./sqrt(etta))));
                        e4=(exp(C1_n).*hermiteH(M-m2+j2,((1i.*beta_n)./sqrt(etta))))+(exp(C2_p).*hermiteH(M-m2+j2,((-1i.*beta_p)./sqrt(etta))));
                        T6=T6+nchoosek(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)).*e1.*e2.*e3.*e4;
                    end
                    T5=T5+T6.*(((factorial(M-m1).*(-1).^s2)./(factorial(s2).*factorial(M-m1-2.*s2))).*(((2.*1i)./(delta.^(0.5))).^(M-m1-2.*s2)));
                end
                T4=T4+T5.*nchoosek(m1-(2.*s2),j1).*((1./(2.*1i.*sqrt(etta))).^(m2+j2)).*((1./(po.^2)).^j1).*((b./2.*wo).^(m1-2.*s1-j1));
            end
            T3=T3+T4.*(((factorial(m1).*(-1).^s1)./(factorial(s1).*factorial(m1-2.*s1))).*(((2.*1i)./(delta.^(0.5))).^(m1-2.*s1)));
        end
        T2=T2+T3*nchoosek(M,m2)*((-1i)^(M-m2));
    end
    T1=T1+T2*nchoosek(M,m1)*((1i)^(M-m1));
end

Error using .*
Array sizes must match.

I0 = -real(X.*T1);

% convert to numerical double
I_numerical = double(vpa(I0))

Answer 1

[KSSV]

Check this:

nchoosek(m1-(2.*s2),j1)

this is giving you an expty matrix.

A = rand(1,5) ;
B = zeros([],1) 

B =

  0×1 empty double column vector

A*B

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix
individually, use TIMES (.*) for elementwise multiplication.

Use nchoosek use wisely.