% Parameters
gamma_values = linspace(0.01, 0.99, 50);
t_values = linspace(0, 2, 100);
difference_values_gamma_t = zeros(length(t_values), length(gamma_values));

% Loop through gamma and t values
for i = 1:length(gamma_values)
    gamma_val = gamma_values(i);
    a_val = (1 - sqrt(gamma_val))^2;
    b_val = (1 + sqrt(gamma_val))^2;
    
    for j = 1:length(t_values)
        t_val = t_values(j);
        
        % Define the density function for Marchenko-Pastur distribution
        density_function = @(s) (1 / (2 * pi * gamma_val * s)) * sqrt(max((b_val - s) * (s - a_val), 0));
        
        % Define the integrands for LHS and RHS
        lhs_integrand = @(s) (1 - gamma_val * (1 - exp(-t_val * s))).^2 .* exp(-t_val * s) .* density_function(s);
        rhs_integrand = @(s) s .* exp(-2 * t_val * s) .* density_function(s);
        
        % Perform numerical integration for LHS and RHS
        lhs_value = integral(lhs_integrand, a_val, b_val, 'ArrayValued', true);
        rhs_value = integral(rhs_integrand, a_val, b_val, 'ArrayValued', true);
        
        % Store the difference between LHS and RHS
        difference_values_gamma_t(j, i) = lhs_value - rhs_value;
    end
end

% Create meshgrid for contour plot
[gamma_values_mesh, t_values_mesh] = meshgrid(gamma_values, t_values);

% Plotting the contour plot for the difference as a function of gamma and t
figure;
contourf(t_values_mesh, gamma_values_mesh, difference_values_gamma_t, 15, 'LineColor', 'none');
% title('Difference between LHS and RHS as a Function of t and gamma');
xlabel('t');
ylabel('gamma');
colorbar;
colormap('winter');
