%
% Resolver um sistema diferencial da forma
%          y'=Ay
% ou "A historia de Romeu e Julieta"
%
clear;
clc;
%
% Matriz A de "amor"
%
%Julieta apaixonada (instavel)
tfinal=30;
A=[-0.1 -0.2
   0.1 0.3];
%Romeu apaixonado (instavel)
tfinal=90;
A=[-0.1 0.8
   -0.1 0.2];
%Romeu apaixonado (estavel)
tfinal=90;
A=[-0.1 0.8
   -0.1 0.0];
%Julieta apaixonada (estavel)
tfinal=400;
A=[0.0 -0.1
   0.8 0.0];
% condicoes iniciais
r0=0.5; j0=0.5;
n=1000;
deltat=tfinal/n;
t=0;
x=[r0 j0];
x0=x;
r=r0; j=j0;
% metodo de Runge-Kutta de quarta ordem
tempo=[0:deltat:tfinal];
for i=1:n
   k1=romeu(x0,A);
   k2=romeu(x0+deltat/2*k1,A);
   k3=romeu(x0+deltat/2*k2,A);
   k4=romeu(x0+deltat*k1,A);
   x0=x0+(deltat/6)*(k1+2*k2+2*k3+k4);
   r=[r x0(1)];
   j=[j x0(2)];
   x=[x x0];
end;
% graficos
set(gcf,'name','Romeu & Julieta');
subplot(211)
plot(tempo,r,'r',tempo,j,'g'), xlabel('Dias'), ylabel('Temperatura')
drawnow
title('A paixao em casa de Romeu e Julieta'), legend('Romeu', 'Julieta')
subplot(212)
for i=1:n
   plot(r(i),j(i),'EraseMode','none','MarkerSize',4)
   xlabel('Romeu'), ylabel('Julieta'), hold on
   drawnow
   axis([min(r)-0.1 max(r)+0.1 min(j)-0.1 max(j)+0.1])
end;
title('Psicanalise')
hold off