clear all; close all; clc;
% Diagrama de fase

A = [1 2;-2 1];
ej4=@(x,t) A*x;
x0 = [1;1];
xs = fsolve(@(x) ej4(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);

[R,V] = eig(A);disp('vectores propios =  ');disp(R);disp('valores propios =  ');disp([V(1,1);V(2,2)])
hold on
if imag(V(1,1)) == 0 && imag(V(2,2)) == 0
  if R(1,1) ~=0 && R(1,2) ~= 0 then
    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))
    plot(xf',[y1' y2'],'k');
   end
end

plot(xs(1),xs(2),'ko')

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