%
%Integrates two coupled systems ODEs with synchronization
%

format compact
set(0,'defaulttextinterpreter','none')

options= odeset('Reltol',1.e-5);
%options= odeset('Reltol',1.e-22,'OutputFcn','odephas3','Refine',222);
%options= odeset('Reltol',1.e-3,'OutputFcn','odephas3','Refine',2);
%options= odeset('Reltol',1.e-3);

 A = 0.1, tf=200;fac=4; DC = 0.03;
 %A = [2 -4 ; 6 -1], tf=25;fac=4; DC = 0.2;
 %A = [2 -4 1; 3 0 -1 ; 0 -1 1]/10, tf=120;fac=4; DC = 0.05;
 
 N=length(A);
 EigA = eig(A)'
 CR = max(real(EigA))/2
 fprintf('Press any key to continue\n');pause
 
 S = 5;
 rng('default');
 x0 = S*(2*(rand(1,N)-0.5)); 
 y0 = S*(2*(rand(1,N)-0.5)); 
 y0 = -x0+(S/3)*(2*(rand(1,N)-0.5));
 x0 = x0+S;
 if(N==1)
  y0 = y0-S;
 end
 v0 = [x0 y0]
 tspan=[0:tf/1000:tf];

 [T1,Y1]=ode45('linearmio',tspan,v0,options,A,CR-DC,N);
 [T2,Y2]=ode45('linearmio',tspan,v0,options,A,CR,N);
 [T3,Y3]=ode45('linearmio',tspan,v0,options,A,CR+DC,N);
 
 ifac=floor(length(tspan)/fac);
 Tfac = T1(1:ifac);
 
 figure(1);
 clf;
 if(N==1)
  subplot(3,2,1)
  plot(Tfac,Y1(1:ifac,1),Tfac,Y1(1:ifac,2))
  axis tight; ylabel(['x and y (C=',num2str(CR-DC),')']);
  subplot(3,2,3)
  plot(Tfac,Y2(1:ifac,1),Tfac,Y2(1:ifac,2))
  axis tight; ylabel(['x and y (C=',num2str(CR),')']);
  subplot(3,2,5);xlabel('t')
  plot(Tfac,Y3(1:ifac,1),Tfac,Y3(1:ifac,2))
  axis tight; ylabel(['x and y (C=',num2str(CR+DC),')']);
  subplot(3,2,2)
  semilogy(T1,abs(Y1(:,1)-Y1(:,2)))
  axis tight; ylabel(['||x-y|| (C=',num2str(CR-DC),')']);
  subplot(3,2,4)
  plot(T2,abs(Y2(:,1)-Y2(:,2)))
  axis tight; ylabel(['||x-y|| (C=',num2str(CR),')']);
  subplot(3,2,6)
  semilogy(T3,abs(Y3(:,1)-Y3(:,2)))
  axis tight; ylabel(['||x-y|| (C=',num2str(CR+DC),')']);
  xlabel('t')
 elseif(N==2)
  subplot(3,2,1)
  plot(Y1(1:ifac,1),Y1(1:ifac,2),Y1(1:ifac,N+1),Y1(1:ifac,N+2))
  hold on
  plot(Y1(1,1),Y1(1,2),'*',Y1(1,N+1),Y1(1,N+2),'*')
  plot(Y1(ifac,1),Y1(ifac,2),'^',Y1(ifac,N+1),Y1(ifac,N+2),'^')
  axis tight; ylabel(['X(y) and Y(y) (C=',num2str(CR-DC),')']);
  subplot(3,2,3)
  plot(Y2(1:ifac,1),Y2(1:ifac,2),Y2(1:ifac,N+1),Y2(1:ifac,N+2))
  hold on
  plot(Y2(1,1),Y2(1,2),'*',Y2(1,N+1),Y2(1,N+2),'*')
  plot(Y2(ifac,1),Y2(ifac,2),'^',Y2(ifac,N+1),Y2(ifac,N+2),'^')
  axis tight; ylabel(['X(y) and Y(y) (C=',num2str(CR),')']);
  subplot(3,2,5);
  plot(Y3(1:ifac,1),Y3(1:ifac,2),Y3(1:ifac,N+1),Y3(1:ifac,N+2))
  hold on
  plot(Y3(1,1),Y3(1,2),'*',Y3(1,N+1),Y3(1,N+2),'*')
  plot(Y3(ifac,1),Y3(ifac,2),'^',Y3(ifac,N+1),Y3(ifac,N+2),'^')
  axis tight; xlabel('X(x) and Y(x)');ylabel(['X(y) and Y(y) (C=',num2str(CR+DC),')']);
  subplot(3,2,2)
  semilogy(T1,sqrt(sum((Y1(:,1:N)-Y1(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');
  subplot(3,2,4)
  plot(T2,sqrt(sum((Y2(:,1:N)-Y2(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');
  subplot(3,2,6)
  semilogy(T3,sqrt(sum((Y3(:,1:N)-Y3(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');xlabel('t')
 elseif(N>2)
  subplot(3,2,1)
  plot3(Y1(1:ifac,1),Y1(1:ifac,2),Y1(1:ifac,3),Y1(1:ifac,N+1),Y1(1:ifac,N+2),Y1(1:ifac,N+3))
  hold on
  plot3(Y1(1,1),Y1(1,2),Y1(1,3),'*',Y1(1,N+1),Y1(1,N+2),Y1(1,N+3),'*')
  plot3(Y1(ifac,1),Y1(ifac,2),Y1(ifac,3),'^',Y1(ifac,N+1),Y1(ifac,N+2),Y1(ifac,N+3),'^')
  axis tight; 
  xlabel(['X_x and Y_x']);
  ylabel(['X_y and Y_y']);
  zlabel(['X_z and Y_z (C=',num2str(CR-DC),')']);
  grid
  subplot(3,2,3)
  plot3(Y2(1:ifac,1),Y2(1:ifac,2),Y2(1:ifac,3),Y2(1:ifac,N+1),Y2(1:ifac,N+2),Y2(1:ifac,N+3))
  hold on
  plot3(Y2(1,1),Y2(1,2),Y2(1,3),'*',Y2(1,N+1),Y2(1,N+2),Y2(1,N+3),'*')
  plot3(Y2(ifac,1),Y2(ifac,2),Y2(ifac,3),'^',Y2(ifac,N+1),Y2(ifac,N+2),Y2(ifac,N+3),'^')
  xlabel(['X_x and Y_x']);
  ylabel(['X_y and Y_y']);
  zlabel(['X_z and Y_z (C=',num2str(CR),')']);
  grid
  subplot(3,2,5)
  plot3(Y3(1:ifac,1),Y3(1:ifac,2),Y3(1:ifac,3),Y3(1:ifac,N+1),Y3(1:ifac,N+2),Y3(1:ifac,N+3))
  hold on
  plot3(Y3(1,1),Y3(1,2),Y3(1,3),'*',Y3(1,N+1),Y3(1,N+2),Y3(1,N+3),'*')
  plot3(Y3(ifac,1),Y3(ifac,2),Y3(ifac,3),'^',Y3(ifac,N+1),Y3(ifac,N+2),Y3(ifac,N+3),'^')
  xlabel(['X_x and Y_x']);
  ylabel(['X_y and Y_y']);
  zlabel(['X_z and Y_z (C=',num2str(CR+DC),')']);
  grid
  subplot(3,2,2)
  semilogy(T1,sqrt(sum((Y1(:,1:N)-Y1(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');
  subplot(3,2,4)
  plot(T2,sqrt(sum((Y2(:,1:N)-Y2(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');
  subplot(3,2,6)
  semilogy(T3,sqrt(sum((Y3(:,1:N)-Y3(:,N+1:2*N)).^2,2)))
  axis tight; ylabel('||x-y||');xlabel('t')
 end

     
