function rj

ar = 1;
pr = .2;
aj = 1;
pj = -.2;

r = .5;
j = .6;

for i=1:200
  
  r(i+1) = ar*r(i) + pr*j(i);
  j(i+1) = aj*j(i) + pj*r(i);
  
end

r(50)
j(50)
r(50)/j(50)

subplot(2,2,1)
plot(r,'.-')
xlabel('time'); ylabel('Romeo');

subplot(2,2,2)
plot(j,'.-')
xlabel('time'); ylabel('Juliet');

subplot(2,2,3)
plot(r,j,'.-')
xlabel('Romeo'); ylabel('Juliet');
title('Phase Space')