Building PDFs of an eigenvalue of Wishart distribution in MATLAB

I have been trying to build MATLAB code from PDFs in a paper, "On the marginal distribution of the eigenvalues of Wishart matrices." Specifically, I am working on

1. the PDF of the largest eigenvalue of the uncorrelated central Wishart matrix
2. the PDF of a generic eigenvalue of the uncorrelated central Wishart matrix

which are Equation (38) in the paper.

This is the PDF of the largest eigenvalue of the uncorrelated central Wishart matrix

This is the PDF of a generic (unordered) eigenvalue of the uncorrelated central Wishart matrix

And I have constructed a MATLAB script for these two PDFs like this.

input_end=30; % Define the length of time domain
x=0:0.01:input_end; % Define time domain
y_lmda1=zeros(1,length(x));   % Define f_lmda1(x)
y_generic=zeros(1,length(x));   % Define f_generic(x)

M=5;
N=5;
nmax=max(M,N);
nmin=min(M,N);

%%% Define K %%%

product_of_K_when_max=1;
for i=1:nmin
    fact_of_prdct_nmax=factorial(nmax-i);
    product_of_K_when_max=product_of_K_when_max*fact_of_prdct_nmax;
end

product_of_K_when_min=1;
for i=1:nmin
    fact_of_prdct_nmin=factorial(nmin-i);
    product_of_K_when_min=product_of_K_when_min*fact_of_prdct_nmin;
end


K=(product_of_K_when_min.*product_of_K_when_max)^(-1);

%%% Define the PDF of the largest eigenvalue %%%


PDF_sum=0;  %%% Define the initial state of the matrix omega_lambda1 %%%
PDF_generic_doublesummation=0;

omega_lambda1=zeros(nmin,nmin);    %%% Allocating initial values to the matrix Omega%%

%%% Define the matrix omega_lambda1 %%%

omega_lambda1=zeros(nmin,nmin); %%% Allocating initial values to the matrix Omega %%%


for k=1:length(x); %% output will be calculated with corresponding x(k)%% 
    for n=1:nmin; 
        for m=1:nmin; 
            for i=1:nmin;
                for j=1:nmin;
                    if i<n && j<m
                        omega_lambda1(i,j)=gamma(i+j-2+nmax-nmin+1)*gammainc(x(k),i+j-2+nmax-nmin+1);
                        %%% gamma() is multiplied because the incomplete Gamma function in MATLAB and the incomplete Gamma function in the paper are different.%%%
                    elseif i>=n && j >= m
                        omega_lambda1(i,j)=gamma(i+j+nmax-nmin+1)*gammainc(x(k),i+j+nmax-nmin+1);
                    else
                        omega_lambda1(i,j)=gamma(i+j-1+nmax-nmin+1)*gammainc(x(k),i+j-1+nmax-nmin+1);
                    end  
                end
            end
            PDF_sum=PDF_sum+((-1).^(n+m)).*(x(k).^(n+m-2+nmax-nmin)).*(exp(-x(k))).*(det(omega_lambda1)); 
            omega_lambda1=zeros(nmin,nmin);
        end
    end
      
    y_lmda1(k)=K.*PDF_sum; %%% After double summations, the normalizing factor is multiplied and it's input to the output vector.%%%
    PDF_sum=0; %%% Initializing the summation loop with new x value. %%%
end
 
%%% Define the matrix omega_generic %%%

sigma_nmin_1=nmin.^(-1/(nmin-1));

omega_generic=zeros(nmin,nmin);
for k=1:length(x);
    for n=1:nmin;
        for m=1:nmin;
            for i=1:nmin;
                for j=1:nmin;
                    if i<n && j<m
                        omega_generic(i,j)=factorial(i+j-2+nmax-nmin)*sigma_nmin_1;
                    elseif i>=n && j >= m
                        omega_generic(i,j)=factorial(i+j+nmax-nmin)*sigma_nmin_1;
                    else
                        omega_generic(i,j)=factorial(i+j-1+nmax-nmin)*sigma_nmin_1;
                    end
                end
            end
            PDF_generic_doublesummation=PDF_generic_doublesummation+...
                ((-1).^(n+m)).*(x(k).^(n+m-2+nmax-nmin)).*(exp(-x(k)))*(det(omega_generic));
            omega_generic=zeros(nmin,nmin);
        end
    end
    
    y_generic(k)=K.*PDF_generic_doublesummation; % After double summations, the normalizing factor is multiplied and it's input to the output vector.
    PDF_generic_doublesummation=0; %%% Initializing the summation loop with new x value.%%
end


figure;

plot(x,y_lmda1,'k-.');
legend("lambda1", "generic",'FontSize', 14);
ylabel('f(x)');
xlabel('x');

%%% PDF verification %%%
sum(y_lmda1)
sum(y_generic)

I normalized the PDFs with K. However, my result is here , which is not valid because both PDF values go way over 1. I have been trying to figure out what is the problem but it seems it is beyond my capability. I want to know what I am missing now from others.