int function is not working properly in matlab

Following is my script to integrate a function in matlab2023a. Here If I put n_S approximately 50 and i=50, the code does work. But for large value as n_S=100 and i=100 it does not work. Could you help me to what could be the reaso n behind that.

clc
clear
del_0_S=4.0493532;
n_S=100;
n_star_S=vpa(n_S-del_0_S); 
Energy_S=vpa(-1/(2*n_star_S^2));
l_s=0;  l_d=2;  s=1/100;          % l values of S, P, D, and F states
syms r
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)))));
del_0_D3=2.47545; 
% 
n_D3=5:115;
for i=100
    n_star_D3(i)=double(n_D3(i)-del_0_D3);% Effective n
    Energy_D3(i)=double((-1/(2*n_star_D3(i)^2)));
    W_f_D3(i)=(whittakerW(n_star_D3(i), l_d+0.5,(2.*(r))./(n_star_D3(i))));
    Gff_D3(i)=(W_f_D3(i)./((sqrt(n_star_D3(i).^2.*gamma(n_star_D3(i)+l_d+1).*gamma(n_star_D3(i)-l_d)))));
    fun_S_D3(i)=(((r).^2).*Gii.*Gff_D3(i));
    r_min_S_D3(i)=double((s*n_S*n_D3(i)/((n_D3(i)+n_S))));
    Matrix_element_S_D3(i)=double((int(fun_S_D3(i),r_min_S_D3(i),25000)));
    
end

Answer 1

[Torsten]

Evaluate

    gamma(n_star_D3(i)+l_d+1)
    gamma(n_star_D3(i)-l_d)    

and you will see the reason.

Consequently, Gff_D3(i) becomes 0.

Answer 2

[Walter Roberson]

Your Gii is pretty small.

Your Gff_D3 is pretty small.

Multiply the two together and you get a value small enough that it is smaller than realmin

Your integral comes out as 0.

Q = @(v) sym(v);

del_0_S = Q(40493532) / Q(10)^7;
n_S = Q(100);
n_star_S = n_S-del_0_S; 
Energy_S = -1/(2*n_star_S^2);
l_s = Q(0);  l_d = Q(2);  s = Q(1/100);          % l values of S, P, D, and F states
syms r
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)))));
del_0_D3 = Q(247545) / Q(10)^5;
% 
n_D3=5:115;
for i=100
    n_star_D3(i) = n_D3(i)-del_0_D3;% Effective n
    Energy_D3(i) = (-1/(2*n_star_D3(i)^2));
    W_f_D3(i) = (whittakerW(n_star_D3(i), l_d+0.5,(2.*(r))./(n_star_D3(i))));
    Gff_D3(i) = (W_f_D3(i)./((sqrt(n_star_D3(i).^2.*gamma(n_star_D3(i)+l_d+1).*gamma(n_star_D3(i)-l_d)))));
    fun_S_D3(i) = (((r).^2).*Gii.*Gff_D3(i));
    r_min_S_D3(i) = (s*n_S*n_D3(i)/((n_D3(i)+n_S)));
    %{
    temp = int(fun_S_D3(i),r_min_S_D3(i),25000);
    Matrix_element_S_D3(i) = double(temp); 
    %}
    vpa(Gii, 10)
    vpa(Gff_D3(i), 10)
    vpa(fun_S_D3(i), 10)
end

ans = 

ans = 

ans = 

%Matrix_element_S_D3(i)

[Anal]

This is the plot of Gii with r. Even though the denomenator of Gii is very small. It compensates with the numarator which is W_i.

Actually I am getting the following error when I choose n_S=100 an i=100. For smaller values of both quantities, the code works.

"Unable to convert expression containing remaining symbolic function calls into double array.

Argument must be expression that evaluates to number.

Error in sym/double (line 729)

Xstr = mupadmex('symobj::double', S.s, 0);"

Is there any other process to solve this problem?

[Anal]

I think that it is better to solve numerically by writting the whole code and not just by calling matlab function.