function [CT, xsol, CR] = compute_correlation_coefficient(Gamma, Lambda)
    % computeJacobian: Computes the Jacobian of a competitive ecological system
    % Input:
    %   Gamma: 2x3 matrix of interaction coefficients for species x1 and x2
    %   Lambda: 3x2 matrix of interaction coefficients for resources R1, R2, R3
    % Output:
    %   J: 5x5 Jacobian matrix evaluated at equilibrium
    
    % Validate input dimensions
    if ~isequal(size(Gamma), [2,3])
        error('Gamma must be a 2x3 matrix');
    end
    if ~isequal(size(Lambda), [3,2])
        error('Lambda must be a 3x2 matrix');
    end
    
    % Define symbolic variables
    syms x1 x2 R1 R2 R3 real
    
    % Define system equations
    Sx1 = -1 + Gamma(1,1)*R1 + Gamma(1,2)*R2 + Gamma(1,3)*R3;
    Sx2 = -1 + Gamma(2,1)*R1 + Gamma(2,2)*R2 + Gamma(2,3)*R3;
    
    Sy1 = 1 - R1 - Lambda(1,1)*x1 - Lambda(1,2)*x2;
    Sy2 = 1 - R2 - Lambda(2,1)*x1 - Lambda(2,2)*x2;
    Sy3 = 1 - R3 - Lambda(3,1)*x1 - Lambda(3,2)*x2;
    
    % Equilibrium conditions (set the growth terms to zero)
    eqns = [Sx1 == 0, Sx2 == 0, Sy1 == 0, Sy2 == 0, Sy3 == 0];
    vars = [R1, R2, R3, x1, x2];
    
    % Solve for equilibrium points
    sol = vpasolve(eqns, vars);
    
    if isempty(sol)
        error('No equilibrium solution found.');
    end
    
    % Extract equilibrium values into a numeric vector
    xsol = [
        double(sol.x1), ...
        double(sol.x2), ...
        double(sol.R1), ...
        double(sol.R2), ...
        double(sol.R3)
    ];
    
     % Display equilibrium values
%      disp('Equilibrium Values:');
%      fprintf('x1 = %.6f, x2 = %.6f, R1 = %.6f, R2 = %.6f, R3 = %.6f\n', ...
%          xsol(1), xsol(2), xsol(3), xsol(4), xsol(5));
%     
    % Define the dynamic equations
    x1d = x1 * (-1 + Gamma(1,1)*R1 + Gamma(1,2)*R2 + Gamma(1,3)*R3);
    x2d = x2 * (-1 + Gamma(2,1)*R1 + Gamma(2,2)*R2 + Gamma(2,3)*R3);
    R1d = R1 * (1 - R1 - Lambda(1,1)*x1 - Lambda(1,2)*x2);
    R2d = R2 * (1 - R2 - Lambda(2,1)*x1 - Lambda(2,2)*x2);
    R3d = R3 * (1 - R3 - Lambda(3,1)*x1 - Lambda(3,2)*x2);
    
    % Combine into a vector field
    F = [x1d; x2d; R1d; R2d; R3d];
    
    % Compute the Jacobian matrix
    J_sym = jacobian(F, [x1, x2, R1, R2, R3]);
    J = double(subs(J_sym, [x1, x2, R1, R2, R3], xsol));
    
%     % Display Jacobian
%     disp('Jacobian Matrix (J):');
%     disp(J);

    % Compute Eigenvalues and Eigenvectors
    [P, D] = eig(J);
    ev = diag(D); % Extract eigenvalues
    PP = inv(P); % Inverse eigenvector matrix
    
%     % Display eigenvalues and matrices
%       disp('Eigenvalues (ev):');
%       disp(ev);
%       disp('Right Eigenvector Matrix (P):');
%       disp(P);
%     disp('Diagonal Matrix (D):');
%     disp(D);
%      disp('Inverse Eigenvector Matrix (PP):');
%      disp(PP);
%  



    % Compute the Correlation Matrix
    Corrmat = zeros(5,5);
    for l = 1:5
        for m = 1:5
            sumVal = 0;
            for i = 1:5
                for j = 1:5
                    for k = 3:5 % Indices 3:5 correspond to R1, R2, R3 in xsol
                        sumVal = sumVal + (P(l,i) * P(m,j) * PP(i,k) * PP(j,k) * xsol(k)^2) / (ev(i) + ev(j));
                    end
                end
            end
            Corrmat(l,m) = sumVal;
        end
    end

    % Ensure numerical stability
    Corrmat = real(Corrmat);
    
%     % Display Correlation Matrix
%     disp('Correlation Matrix (Corrmat):');
%     disp(Corrmat);

   CT =  -Corrmat(1,2)/sqrt(Corrmat(1,1)*Corrmat(2,2));
   
   
   % from now on, compute corr coeff for growth rates
   
   T1 = 0;
   T2 = 0; 
   T12 = 0;
   for i=1:3
       for j = 1:3
   T1 = T1 + Gamma(1,i)*Gamma(1,j)*Corrmat(2+i,2+j);
   T2 = T2 + Gamma(2,i)*Gamma(2,j)*Corrmat(2+i,2+j);
   T12 = T12 + Gamma(1,i)*Gamma(2,j)*Corrmat(2+i,2+j);
       end
   end
   
   CR = -T12/sqrt(T1*T2);
end