function sol_3D()
    
    %
    % Draw the solution of a 3D system
    %
    
    close all
    
     Parameters
    b = 8/3;
    sigma = 10;
    r = 0.5;
    
     Draw the steady states
    figure
    hold on
    plot3(0, 0, 0, 'ro')
    plot3(sqrt(b*(r-1)), sqrt(b*(r-1)), r-1, 'ro')
    plot3(-sqrt(b*(r-1)), -sqrt(b*(r-1)), r-1, 'ro')
    
     Plot the solution in the phase space    
    plot3(x, y, z, 'g', 'LineWidth', 2) 
    hold on
    plot3(x0, y0, z0, 'g*') % Initial condition
    
    xlabel('x')
    ylabel('y')
    zlabel('z')
    grid on
    
    %% Plot the solution over time
    figure
    plot(t, x, 'r')
    hold on
    plot(t, y, 'b')  
    plot(t, z, 'k') 
    
    legend('x(t)', 'y(t)', 'z(t)')
    xlabel('t')
    
end
 
% This function corresponds to the right hand side of the ODE
function du = f(t, u, sigma, r, b)
    x = u(1);
    y = u(2);
    z = u(3);
    
    dx = sigma*(y-x);
    dy = -x*z+r*x-y;
    dz = x*y-b*z;    
    
    du = [dx; dy; dz];
end