function [t, v] = Lorenz9_sim(x0,r,t_end)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% Numerical integration and visualization of the Lorenz system in 9D %
% Input variables:                                                   %
% x0           -              initial condition                      %
% r            -              Rayleigh parameter                     %
% t_transient  -              Length of Integration                  %
% By Franz Bozsak                                                    %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% declaration                                                        %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% Model parameters
sigma  = 10;
b      = 8/3;
a      = 4/b-1;          % this is a^2
c1     = (4+a)/(1+a);
c2     = (9+a)/(1+a);

i = 1;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% Integration using a Runge-Kutta Method                             %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

[t, v] = ode45('lorenz9ode', t_end, [x0 x0], [], sigma, r, b,c1,c2);

% Split out the trajectory matrix into x,y,z,x1,y1,z1,x2,y2,z2 vector variables
x = v(:, 1);
y = v(:, 2);
z = v(:, 3);
x1 = v(:, 4);
y1 = v(:, 5);
z1 = v(:, 6);
x2 = v(:, 7);
y2 = v(:, 8);
z2 = v(:, 9);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% Visualization                                                      %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %

% % Plot trajectory in 3x3D state-space 
% figure;
% plot3(x, y, z, 'b-');
% title('Lorenz State Space (X,Y,Z)');
% xlabel('x(t)');
% ylabel('y(t)');
% zlabel('z(t)');
% grid on;
% rotate3d on;
% view([37,26]);
% 
% % figure;
% % plot3(x1, y1, z1, 'b-');
% % title('Secondary State Space (X1,Y1,Z1)');
% % xlabel('x1(t)');
% % ylabel('y1(t)');
% % zlabel('z1(t)');
% % grid on;
% % rotate3d on;
% % view([37,26]);
% 
% figure;
% plot3(x2, y2, z2, 'b-');
% title('Secondary State Space (X2,Y2,Z2)');
% xlabel('x2(t)');
% ylabel('y2(t)');
% zlabel('z2(t)');
% grid on;
% rotate3d on;
% view([37,26]);
% 
% % Display time series plots
% figure;
% subplot(3,3,1);
% semilogy(t, x, 'b-');
% title('Time series');
% xlabel('t');
% ylabel('x(t)');
% % 
% subplot(3,3,2);
% semilogy(t, y, 'b-');
% xlabel('t');
% ylabel('y(t)');
% % 
% subplot(3,3,3);
% semilogy(t, z, 'b-');
% xlabel('x(t)');
% ylabel('z(t)');
% % 
% subplot(3,3,4);
% semilogy(t, x1, 'b-');
% xlabel('t');
% ylabel('x1(t)');
% % 
% subplot(3,3,5);
% semilogy(t, y1, 'b-');
% xlabel('t');
% ylabel('y1(t)');
% % 
% subplot(3,3,6);
% semilogy(t, z1, 'b-');
% xlabel('t');
% ylabel('z1(t)');
% % 
% subplot(3,3,7);
% semilogy(t, x2, 'b-');
% xlabel('t');
% ylabel('x2(t)');
% % 
% subplot(3,3,8);
% semilogy(t, y2, 'b-');
% xlabel('t');
% ylabel('y2(t)');
% % 
% subplot(3,3,9);
% semilogy(t, z2, 'b-');
% xlabel('t');
% ylabel('z2(t)');
% 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% % Poincare Maps                                                      %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% 
% % Now find next maxima to create next-maximum map of z(t)
% max_points = find((z(2:end-1) > z(3:end)) & (z(2:end-1) > z(1:end-2))) + 1;
% z_max      = z(max_points);
% t_max      = t(max_points);
% 
% % Plot Poincare map
% figure;
% plot(z_max(1:end-1), z_max(2:end), 'r.');
% title('Lorenz (z) poincare map');
% xlabel('z_{max}(n)');
% ylabel('z_{max}(n+1)');
% 
% % Now find next maxima to create next-maximum map of z2(t)
% max2_points = find((z2(2:end-1) > z2(3:end)) & (z2(2:end-1) > z2(1:end-2))) + 1;
% z2_max      = z2(max_points);
% t2_max      = t(max_points);
% 
% % Plot Poincare map
% figure;
% plot(z2_max(1:end-1), z2_max(2:end), 'r.');
% title('z2 Poincare map');
% xlabel('z2_{max}(n)');
% ylabel('z2_{max}(n+1)');

