function stability_3D()
 
    close all
 
     Define the range of the varying parameter
    nPoints = 100; % Resolution
    r = linspace(0, 30, nPoints); 
    
     For each value of the varying parameter, study the stability of a steady state
    for i=1:nPoints
        
         Extract the coordinates of the selected steady state
        x = E(1);
        y = E(2);
        z = E(3);
        
         Find the eigenvalues of J
        eigenvalues(:,i) = eig(J);
 
    end
    
    %% Plot
    figure
    
    subplot(1,2,1)
    for i=1:3
        plot(r, real(eigenvalues(i,:)), 'LineWidth', 2);
        hold on
    end
    title('Real part of eigenvalues');
    xlabel('r')
    legend('\lambda_1', '\lambda_2', '\lambda_3')
    grid on
    
    subplot(1,2,2)
    for i=1:3
        plot(r, imag(eigenvalues(i,:)), 'LineWidth', 2);
        hold on
    end
    title('Imaginary part of eigenvalues');
    xlabel('r')    
    legend('\lambda_1', '\lambda_2', '\lambda_3')
    grid on
end