function RJ

%%
%simulation of RJ model: dR/dt = a R + b J
%                        dJ/dt = c R + d J
  
%Define matrix: replace below with the relevant values
a=2;
b=0;
c=0;
d=3;


%a=.2;b=2;c=-2;d=0;
%a=1;b=1;c=4;d=1;

A=[a b ;...
   c d]; 


%%
%phase diagram :

R=-2:.2:2; LR=length(R); 
J=-2:.2:2; LJ=length(J); 
[matR,matJ]=meshgrid(R,J); %matrices with R in columns and J in rows (size LJxLR)


for iR=1:LR
    for iJ=1:LJ
        RJdot=A*[R(iR) ;J(iJ)];
        Rdot(iJ,iR)=RJdot(1); 
        Jdot(iJ,iR)=RJdot(2);
    end
end

phase_diagram=figure('position',[400 400 800 600]);
    quiver(matR,matJ,Rdot,Jdot); %all LJxLR
    grid on; line([0 0],ylim,'color','k'); line(xlim,[0 0],'color','k'); 
    xlabel('R')
    ylabel('J')
    title(['Phase diagram for A=[',num2str(A(1,:)),' ; ',num2str(A(2,:)),']']);


%%
%Compute a trajectory with initial condition X0 :

X0=[1;1];
Tmax=100;
dt=.01; 
t=0:dt:Tmax; 
%initial conditions: 
X=X0;
R(1)=X0(1);
J(1)=X0(2); 
%time evolution:
for it=2:length(t)
  X = X + dt*A*X;
  R(it)=X(1);
  J(it)=X(2);
end   
  
figure('position',[200 200 800 600]);
    subplot(211); %time evolution of R(t) and J(t) 
        plot(t,R,'b'); hold on; plot(t,J,'m--'); 
        plot(t(1),X0(1),'bx');plot(t(1),X0(2),'mo');
        legend('R(t)','J(t)',['R(0)=',num2str(X0(1))],['J(0)=',num2str(X0(2))]);
        title([{'Time evolution of R(t) and J(t)'}, ...
               {['A=[',num2str(A(1,:)),' ; ',num2str(A(2,:)),']']}]);
        xlabel('time');
        ylabel('R(t) and J(t)');
        xlim([0 Tmax]); 
        hold off;
    subplot(212); plot(R,J,'r'); %trajectory in the phase diagram 
        hold on; plot(X0(1),X0(2),'rx','linewidth',2);
        title('Phase diagram: R versus J (red cross=initial condition)');
        xlabel('R');
        ylabel('J');
        hold off;


figure(phase_diagram); %Ajout de cette trajectoire dans le diagramme des phases
xx=xlim; yy=ylim; hold on;plot(R,J,'r'); plot(R(1),J(1),'rx'); xlim(xx); ylim(yy); 

end