function Y = run_kepler_ode2(x,y,z,x2,y2,z2)


%clear
%clf
global rho sigma beta
rho = 28;
sigma = 10;
beta = 8/3;

% initial conditions
% x = 1;
% y = 1;
% z = 1;


Y = [x;y;z]; % by convention, should be a column vector
X = [x;y;z];
Z = [x;y;z];
Y2 = [x2;y2;z2]; % by convention, should be a column vector
X2 = [x2;y2;z2];
Z2 = [x2;y2;z2];

t    =  0.0;
tmax = 1.0;
dt   =  0.01;
n    = 10000; % number of steps

tt = zeros(n+1,1);
Xt = zeros(n+1,size(X,1));
Yt = zeros(n+1,size(Y,1));
Zt = zeros(n+1,size(Z,1));
tt(1)   = t;
Yt(1,:) = Y';
Xt(1,:) = X';
Zt(1,:) = Z';
X2t = zeros(n+1,size(X2,1));
Y2t = zeros(n+1,size(Y2,1));
Z2t = zeros(n+1,size(Z2,1));
Y2t(1,:) = Y2';
X2t(1,:) = X2';
Z2t(1,:) = Z2';
fff = 0;
for i=1:n

    f1 = kepler_ode2X(t     , X                );
    f2 = kepler_ode2X(t+dt/2, X+f1 .*  (dt/2)  );
    f3 = kepler_ode2X(t+dt/2, X+f2 .*  (dt/2)  );
    f4 = kepler_ode2X(t+dt  , X+f3 .*   dt     );
   
    X = X + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Xt(i+1,:) = X'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    f1 = kepler_ode2Y(t     , Y                );
    f2 = kepler_ode2Y(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = kepler_ode2Y(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = kepler_ode2Y(t+dt  , Y+f3 .*   dt     );
   
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Yt(i+1,:) = Y';
    
    
    f1 = kepler_ode2Z(t     , Z                );
    f2 = kepler_ode2Z(t+dt/2, Z+f1 .*  (dt/2)  );
    f3 = kepler_ode2Z(t+dt/2, Z+f2 .*  (dt/2)  );
    f4 = kepler_ode2Z(t+dt  , Z+f3 .*   dt     );
   
    Z = Z + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Zt(i+1,:) = Z'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    X(2) = Y(2);
    X(3) = Z(3);
    Y(1) = X(1);
    Y(3) = Z(3);
    Z(1) = X(1);
    Z(2) = Y(2);
    
    f1 = kepler_ode2X(t     , X2                );
    f2 = kepler_ode2X(t+dt/2, X2+f1 .*  (dt/2)  );
    f3 = kepler_ode2X(t+dt/2, X2+f2 .*  (dt/2)  );
    f4 = kepler_ode2X(t+dt  , X2+f3 .*   dt     );
   
    X2 = X2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    X2t(i+1,:) = X2'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    f1 = kepler_ode2Y(t     , Y2                );
    f2 = kepler_ode2Y(t+dt/2, Y2+f1 .*  (dt/2)  );
    f3 = kepler_ode2Y(t+dt/2, Y2+f2 .*  (dt/2)  );
    f4 = kepler_ode2Y(t+dt  , Y2+f3 .*   dt     );
   
    Y2 = Y2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Y2t(i+1,:) = Y2';
    
    
    f1 = kepler_ode2Z(t     , Z2                );
    f2 = kepler_ode2Z(t+dt/2, Z2+f1 .*  (dt/2)  );
    f3 = kepler_ode2Z(t+dt/2, Z2+f2 .*  (dt/2)  );
    f4 = kepler_ode2Z(t+dt  , Z2+f3 .*   dt     );
   
    Z2 = Z2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Z2t(i+1,:) = Z2'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    X2(2) = Y2(2);
    X2(3) = Z2(3);
    Y2(1) = X2(1);
    Y2(3) = Z2(3);
    Z2(1) = X2(1);
    Z2(2) = Y2(2);

    
    t = t+dt;
    
    tt(i+1)   = t;  % this records the time at each step
    % this records the vector Y at each step --> Yt(:,1) has all x coordinates

if rem(i,100)==0
    subplot(1,3,1)
plot(Xt(:,1),Yt(:,2),X2t(:,1),Y2t(:,2),'m')
title('Y(t) v X(t)')
ylabel('Y(t)')
xlabel('X(t)')
    
   
subplot(1,3,2)
plot3(Xt(:,1),Yt(:,2),Zt(:,3),X2t(:,1),Y2t(:,2),Z2t(:,3),'m')
title('3D plot')
zlabel('Z(t)')
ylabel('Y(t)')
xlabel('X(t)')
%     if sign(Zt(i,3))~=sign(Zt(i+1,3))
% %     hold on
%         subplot(1,3,3)
%          plot(Zt(i,3),Zt(i+1,),'ro')
%          pause(0.000001)
%      end
fff = fff+1;
saveas(gcf,['lorenz',sprintf('%03d',fff),'.png'])
end

end

hold on

for i = 1:100
    rho = i;
    for j = 1:50
    f1 = kepler_ode2X(t     , X                );
    f2 = kepler_ode2X(t+dt/2, X+f1 .*  (dt/2)  );
    f3 = kepler_ode2X(t+dt/2, X+f2 .*  (dt/2)  );
    f4 = kepler_ode2X(t+dt  , X+f3 .*   dt     );
   
    X = X + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Xt(i+1,:) = X'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    f1 = kepler_ode2Y(t     , Y                );
    f2 = kepler_ode2Y(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = kepler_ode2Y(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = kepler_ode2Y(t+dt  , Y+f3 .*   dt     );
   
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Yt(i+1,:) = Y';
    
    
    f1 = kepler_ode2Z(t     , Z                );
    f2 = kepler_ode2Z(t+dt/2, Z+f1 .*  (dt/2)  );
    f3 = kepler_ode2Z(t+dt/2, Z+f2 .*  (dt/2)  );
    f4 = kepler_ode2Z(t+dt  , Z+f3 .*   dt     );
   
    Z = Z + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Zt(i+1,:) = Z'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    X(2) = Y(2);
    X(3) = Z(3);
    Y(1) = X(1);
    Y(3) = Z(3);
    Z(1) = X(1);
    Z(2) = Y(2);
    
        f1 = kepler_ode2X(t     , X2                );
    f2 = kepler_ode2X(t+dt/2, X2+f1 .*  (dt/2)  );
    f3 = kepler_ode2X(t+dt/2, X2+f2 .*  (dt/2)  );
    f4 = kepler_ode2X(t+dt  , X2+f3 .*   dt     );
   
    X2 = X2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    X2t(i+1,:) = X2'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    f1 = kepler_ode2Y(t     , Y2                );
    f2 = kepler_ode2Y(t+dt/2, Y2+f1 .*  (dt/2)  );
    f3 = kepler_ode2Y(t+dt/2, Y2+f2 .*  (dt/2)  );
    f4 = kepler_ode2Y(t+dt  , Y2+f3 .*   dt     );
   
    Y2 = Y2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Y2t(i+1,:) = Y2';
    
    
    f1 = kepler_ode2Z(t     , Z2                );
    f2 = kepler_ode2Z(t+dt/2, Z2+f1 .*  (dt/2)  );
    f3 = kepler_ode2Z(t+dt/2, Z2+f2 .*  (dt/2)  );
    f4 = kepler_ode2Z(t+dt  , Z2+f3 .*   dt     );
   
    Z2 = Z2 + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    Z2t(i+1,:) = Z2'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    X2(2) = Y2(2);
    X2(3) = Z2(3);
    Y2(1) = X2(1);
    Y2(3) = Z2(3);
    Z2(1) = X2(1);
    Z2(2) = Y2(2);


    
    t = t+dt;
    
    tt(i+1)   = t;
    
 
    hold on
        if sign(Yt(i,2))~=sign(Yt(i+1,2))

        subplot(1,3,3)
         plot(rho,Y(2),'b*','MarkerSize',20);
        title('Bifurcation Plot')
         xlabel('Rho')
         ylabel('Y(t)')
        end
        if sign(Y2t(i,2))~=sign(Y2t(i+1,2))

        subplot(1,3,3)
         plot(rho,Y(2),'m*','MarkerSize',15);
        
        end
        
        
    end
    fff = fff +1;
 saveas(gcf,['lorenz',sprintf('%03d',fff),'.png'])
end


end
