matlab does not give any answer.

Here is my code:

tic

clc; all clear;

a=18897.2598;

w1=2*a; % in a0

w2=3.78*a; % in a0

lambda=780*0.001*a; %in a0

z_R1 =(pi/lambda).*(w1.^2); % in micrometer

z_R2 =(pi/lambda).*(w2.^2);

R=1:1000:30001; % in a0

angle_CM=0:1:1;

[R_CM,Cos_theta_R] = meshgrid(R,angle_CM);

l_S=0; l_P=1; s=1/100;

syms r

syms theta

for i=1

for j=1:2

Part_1_1(i,j)= vpa(1+((1./z_R1).*(R_CM(i,j).*Cos_theta_R(i,j)+r.*cos(theta))).^2);

Part_2_1(i,j)= vpa(exp(((-1/(w1.^2)).*(R_CM(i,j).^2).*(1-(Cos_theta_R(i,j).^2)))./Part_1_1(i,j).^2));

Part_3_1(i,j)= vpa(exp(((-1/(w1.^2)).*(r.^2).*((sin(theta)).^2))./Part_1_1(i,j).^2));

Part_4_1(i,j) =vpa(exp(((-1/(w1.^2)).*2.*R_CM(i,j).*r.*sin(theta).*sqrt(1-Cos_theta_R(i,j).^2))./Part_1_1(i,j).^2));

Total_1(i,j)=vpa(((Part_1_1(i,j).^(-0.5)).*(Part_2_1(i,j).*Part_3_1(i,j).*Part_4_1(i,j))));

Part_1_2(i,j)= vpa(1+((1./z_R2).*(R_CM(i,j).*Cos_theta_R(i,j)+r.*cos(theta))).^2);

Part_2_2(i,j)= vpa(exp(((-1/(w2.^2)).*(R_CM(i,j).^2).*(1-(Cos_theta_R(i,j).^2)))./Part_1_2(i,j).^2));

Part_3_2(i,j)= vpa(exp(((-1/(w2.^2)).*(r.^2).*((sin(theta)).^2))./Part_1_2(i,j).^2));

Part_4_2(i,j) =vpa(exp(((-1/(w2.^2)).*2.*R_CM(i,j).*r.*sin(theta).*sqrt(1-Cos_theta_R(i,j).^2))./Part_1_2(i,j).^2));

Total_2(i,j)=vpa(((Part_1_2(i,j).^(-0.5)).*(Part_2_2(i,j).*Part_3_2(i,j).*Part_4_2(i,j))));

Total(i,j) =vpa((abs((Total_1(i,j)-Total_2(i,j)))).^2);

del_0_S=3.1311806; del_2_S=0.1786; del_4_S=0; del_6_S=0; del_8_S=0; % For n S_1/2 state of Cs, Quantum defect parameter

n_S=10;

n_star_S=double(n_S-del_0_S-(del_2_S./((n_S-del_0_S).^2))-(del_4_S./((n_S-del_0_S).^4))...

-(del_6_S./((n_S-del_0_S).^6))-(del_8_S./((n_S-del_0_S).^8))); % Effective n

W_i= (whittakerW(n_star_S, l_S+0.5,(2.*(r))./(n_star_S)));

Gii= W_i./((sqrt(n_star_S.^2.*gamma(n_star_S+l_S+1).*gamma(n_star_S-l_S))));

%include angular part of S1/2, i.e.,

%Yj,MJ(theta,phi)=CG*Yl,ml(theta,phi)=((l+mj+0.5)/(2l+1))^0.5.*Yl,ml(theta,phi)

%For S1/2 state anular part=1*Y0,0(theta,phi)=Y0,0(theta,phi)=0.5*sqrt(1/pi)

fun_S_S(i,j)= 2*pi*(r.^2).*Gii.*Gii.*sin(theta).^2.*Total(i,j).*((0.5*sqrt(1/pi)).^2);

r_min=s*n_S*n_S/(n_S+n_S);

Matrix_element(i,j)=double(int(int(fun_S_S(i,j), 0, pi),r_min,1000));

end

end

toc