clear all; close all; clc;
% Diagrama de fase
A = [2 1;2 -1];
ej3=@(x,t) A*x;

x0 = [1;1];
xs = fsolve(@(x) ej3(x,0) , x0);

[x,y]= meshgrid(-8:8);
xprima1=@(x) A(1,:)*x;
xprima2=@(x) A(2,:)*x;
for i=1:size(x,1)
  for j=1:size(x,2)
    tangen1(i,j)= xprima1([x(i,j);y(i,j)]);
    tangen2(i,j)= xprima2([x(i,j);y(i,j)]);
  end
end

h=quiver(x,y,tangen1,tangen2);
axis([-8,8,-8,8]);
xlabel('x1');ylabel('x2')
set (h, "maxheadsize", 0.3);


xf = linspace(-8,8,8);yf = linspace(-8,8,6);
[R,diagevals] = eig(A);
pend1 = R(2,1)/R(1,1);y1 = xs(2) + pend1*(xf-xs(1));
pend2 = R(2,2)/R(1,2);y2 = xs(2) + pend2*(xf-xs(1));
hold on
plot(xf',[y1' y2'],'k');
plot(xs(1),xs(2),'ko');

t = 0:0.1:10;
for i = 1:16
  [x10,x20] = ginput(1);
  sol = lsode(ej3,[x10;x20],t);
  plot(sol(:,1),sol(:,2),'k');
end
axis([-8,8,-8,8]);