function flow2D_ode_linear()
    
    %
    % Draw the flow of a 2D linear system
    %
    
    close all
    
     Resolution
    nPoints = 10;
    
     Evaluate the flow for each point
    dX = a*X + b*Y;
    dY = c*X + d*Y;
    
     Draw the steady states
    hold on
    plot(0, 0, 'ro')
    
     Evaluate the solution
    x0 = 0.6; % Initial condition for x
    y0 = 0.7; % Initial condition for y
    T = 8; % Time horizon
    [t, u] = ode45(@f, [0, T], [x0, y0], [], a, b, c, d); % Solve!
    x = u(:, 1);
    y = u(:, 2);
    
    plot(x, y, 'c', 'LineWidth', 2) % Draw the solution in the phase space
    plot(x0, y0, 'go') % Draw the initial condition
    
    xlim([xmin xmax])
    ylim([ymin ymax])
    xlabel('x')
    ylabel('y')
    legend('Flow', 'Steady state', 'Nullcline 1', 'Nullcline 2', 'Trajectory', 'Initial condition');
    
    %% Draw the solution over time
    figure
    plot(t, x)
    hold on
    plot(t, y)
    xlabel('t')
    legend('x(t)', 'y(t)')
    
end
 
function du = f(t, u, a, b, c, d)
    x = u(1);
    y = u(2);
    
    dx = a*x + b*y;
    dy = c*x + d*y;
    
    du = [dx; dy];
end