%  -----------------------------------------------------------
%                       Quasiperiodic Motion
%  -----------------------------------------------------------
function example_3

close all;
clear all;
clc;

%  ------- Parameters  -------
tmin        = 0.0;   
tmax        = 80.0;
h           = 0.1;
trans       = 10.0;
tspan_trans = [tmin:h:trans];
tspan       = [tmin:h:tmax];

%  ---------------  Integrate ODE  ----------------
V0 = [0.2,0.6,-0.5,0.2]';
[t_trans,v_trans] = ode45(@rhs_ode,tspan_trans,V0);

V0    = v_trans(end,:)'; 
[t,v] = ode45(@rhs_ode,tspan,V0);


%  ---------------  Plot Results  --------------
xmin = -5.0;  
xmax = 5.0;

%  ------  Time Series  -------
figure(1);
subplot(4,1,1);
plot(t,v(:,1),'r','LineWidth',2);
ylabel('{x_1}');
axis([tmin tmax xmin xmax]);

subplot(4,1,2);
plot(t,v(:,2),'g','LineWidth',2);
ylabel('{x_2}');
axis([tmin tmax xmin xmax]);

subplot(4,1,3);
plot(t,v(:,3),'b','LineWidth',2);
ylabel('{x_3}');
axis([tmin tmax xmin xmax]);

subplot(4,1,4);
plot(t,v(:,4),'b','LineWidth',2);
ylabel('{x_4}');
axis([tmin tmax xmin xmax]);


%  ------  Phase-Space  -------
figure(2);
plot3(v(:,1),v(:,2),v(:,3));
xlabel('x1'); ylabel('x2'); zlabel('x3');
grid on;

figure(3);
plot3(v(:,1),v(:,2),v(:,4));
xlabel('x1'); ylabel('x2'); zlabel('x4');
grid on;



function f = rhs_ode(t,v)

  % --------  Matrix A  ---------
  A = [10.67221819995419 -16.04166063326462   5.40277576664376 -12.25971759661399;
      5.34722021108821  -8.13749698329898   2.75694343887744  -3.65972081442841;
      2.15740660295380  -3.51944323776403   1.19536996814357  -0.31435165073845;
      5.32870169256969  -7.99305253885453   2.66435084628484  -3.73009118479878];
  
  % B = [0 -0.1 0 0; 0.1 0 0 0; 0 0 0 4.7958; 0 0 -4.7958 0];
  x = v;
  
  f = A*x;