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

%deff('[xprima] = ej1(t,x)','xprima = [-1 0;0 -4]*x') ;
ej1=@(x,t) [-1 0;0 -4]*x;
t = 0:0.1:100;
xf = linspace(-4,4,10);yf = linspace(-4,4,6);
for i = 1:length(xf)
    for j = 1:length(yf)
      sol = lsode(ej1,[xf(i);yf(j)],t);
      plot(sol(:,1),sol(:,2),'k');
      hold on
    end
end

[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);