%
%   Systems of Discrete-Time Equations
%
%     R[n+1] = aR[n] + pJ[n]
%     J[n+1] = bJ[n] + qR[n]
%
%   MATH 499, Feb. 2007, Jae-Hun Jung

R0 = input(' initial values of R0 of J0, [R0; J0] ');
a = input(' coefficients of amplication matrix [a1 a2; a3 a4] ');
n = 1; 
nmax = 300;

subplot(1,2,1)
plot(n,R0(1),'ro-',n,R0(2),'bx-'); axis([0 nmax -nmax nmax])
xlabel('n'), ylabel('R_n')
hold on
subplot(1,2,2)
plot(R0(1),R0(2),'.-r'); axis([ -nmax nmax -nmax nmax])
xlabel('R_n(1)'), ylabel('R_n(2)')
hold on
pause

for n=2:nmax
R1 = a*R0; x = [n;n-1]; y1 = [R1(1); R0(1)]; y2 = [R1(2); R0(2)];
subplot(1,2,1)
plot(x,y1,'ro-',x,y2,'bx-'); axis([0 nmax -nmax nmax])
subplot(1,2,2)
plot(y1,y2,'.-r'); axis([ -nmax nmax -nmax nmax])
pause
R0 = R1;
end