%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Math 3520: Differential Equations and Dynamical Systems
%  Compute and plot the Rossler system strange attractor. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%% Runge-Kutta 4th order solver in 3D:
function X = rk4b(X0,F,steps,h)
 % Input
 %  X0 = [x0,y0,z0]   Row vector, initial position
 %  F  = function handle, computes Xdot(x,y,z)
 %  steps = number of time steps to compute (e.g., 100)
 %  h = time increment per step (e.g., 1/steps )
 % Output
 %  X = (1+steps)x(3) matrix of (x,y,z) values along
 %                    the trajectory starting at X0
  X = zeros(1+steps,3);  % Solution to be computed
  X(1,:) = X0;  % Initial point
  for n=1:steps
    f1 = F(X(n,:));
    f2 = F(X(n,:)+f1*h/2);
    f3 = F(X(n,:)+f2*h/2);
    f4 = F(X(n,:)+f3*h);
    X(n+1,:) = X(n,:) +h*(f1+2*f2+2*f3+f4)/6;
  end
end

%%%% From Strogatz, 3rd edition, p.472
X0=[1,1,1]; % same as Lorenz's initial value
a=0.2; % from Strogatz
b=0.2; % from Strogatz, same b=a
c=5.7; % from Strogatz, to get chaos
Xdot = @(V) [-V(2)-V(3), V(1)+a*V(2), b+V(3).*(V(1)-c)];

%%%% Trajectories
h=0.01;   % Use same time step for same accuracy
X = rk4b(X0, Xdot, 20000,h);  % Complete trajectory to t=200.
figure; hold on; 
init=1:100;   plot3(X(init,1),X(init,2),X(init,3));  % 0<t<1;
pause(2);
next=101:1000; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
next=1001:2000; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
next=2001:3000; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
next=3001:5000; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
next=5001:10000; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
next=10001:20001; plot3(X(next,1),X(next,2),X(next,3)); pause(2);
xlabel("X"); ylabel("Y"); zlabel("Z");
title("Trajectory for the Rossler System")