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

ej2=@(x,t) [-1 0;0 4]*x;
t = 0:0.1:10;
xf = linspace(-4,4,10);yf = linspace(-4,4,6);

for i = 1:length(xf)
  sol = lsode(ej2,[xf(i);0.1],t);
  plot(sol(:,1),sol(:,2),'k');
  hold on
  sol = lsode(ej2,[xf(i);-0.1],t);
  plot(sol(:,1),sol(:,2),'k');
end

sol = lsode(ej2,[xf(1);0],t);
plot(sol(:,1),sol(:,2),'k');

sol = lsode(ej2,[0.01;0.01],t);
plot(sol(:,1),sol(:,2),'k');

sol = lsode(ej2,[-0.01;-0.01],t);
plot(sol(:,1),sol(:,2),'k');

sol = lsode(ej2,[xf(end);0],t);
plot(sol(:,1),sol(:,2),'k');

    
[x,y]= meshgrid(-4:0.5:4);
xprima1=@(x) [-1 0]*x;
xprima2=@(x) [0 4]*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,0.5);
xlabel('x1');ylabel('x2')
axis([-4,4,-4,4]);
set (h, "maxheadsize", 0.3);