% Dinamica de sistemas
% Dx = A x  com  x(0) = xo

% exemplos:

A=[-0.5 1;0 -2];
%A = [-0.01 -5; 4 -5]        % (-a +/- b i)
%A = [-0.01 -5; 4 5]         % (a +/- b i)
%A = [-0.01 -5; 4 -0.02]     % (~0 +/- b i)
%A = [5 -1; -2 -4]           % (-a, b)
%A = [-5 -1; -2 -4]          % (-a, -b)
%A = [5 -1; -2 4]            % (a, b)
%A = [3 4; -5 -3]            % (+/- b i)
%A = [-0.4 0.6; -0.2 0.3]    % (-a, 0)
%A = [0.4 0.6; 0.2 0.3]      % (a, 0)

[V,D] = eig(A)              % A = V * D * inv(V)
iV = inv(V);

Dm = min(abs(diag(D)));
if Dm < 1e-8, Dm = 1 / (2 * pi), end;
tf = 9 / Dm;
dt = tf / 50;

% condicao inicial

if min(diag(D)) > 0, f = 0.1; else f = 1; end

S = 'ymcrgbw';
xx = -1:0.2:1; yy = [1;-1];
XXo = [];

%for j = 1:2,
%for i = 1:11,
for i = 1:30,
%   XXo = [XXo; xx(i) yy(j)];
   XXo = [XXo; (2 .* rand(1,2) - 1) * f];
end
%end

n = size(XXo,1);

% dominio do tempo

figure(1); clf;
subplot(211); plot(0,0); AXIS([0 tf -1 1]); hold;
subplot(212); plot(0,0); AXIS([0 tf -1 1]); hold;

for i = 1:min(10,n),
   Xo = XXo(i,:)';
   X = [Xo'];
   t = [0];
   for ti = dt : dt : tf,
      X = [X; real((V * expm (D * ti) * iV * Xo)')];
      t = [t; ti];
   end

   subplot(211); plot(t,X(:,1),S(rem(i,7)+1));
   subplot(212); plot(t,X(:,2),S(rem(i,7)+1));
end

% plano de fases

figure(2); clf;
plot(0,0,'ro'); AXIS([-1 1 -1 1]); hold;

for i = 1:n,
   Xo = XXo(i,:)';
   X = [Xo'];
   t = [0];
   for ti = dt : dt : tf,
      Xt = real((V * expm (D * ti) * iV * Xo)');
      if max(abs(Xt)) > 1, break; end;
      X = [X; Xt];
      t = [t; ti];
   end

   cometa(X(:,1),X(:,2));
%   plot(X(:,1),X(:,2),S(rem(i,7)+1));
end
