function [t,Xh] = lorenz(Pr,b,r,t1,t2,XYZinit,plotting)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function [t,Xh] = lorenz(Pr,b,r,t1,t2,XYZinit,plotting)
%
%   INPUT
%   Pr          Prandtl number
%   b           wave number
%   r           pseudo-Rayleigh parameter
%   t1          start time for plotting and history recording
%   t2          end time for plotting and history recording
%   XYZinit     1 x 3 vector with initial conditions for (X,Y,Z)
%   plotting    flag for on-line visualization: 1 .. plot ; 0 .. no plot
%
%   OUTPUT
%   t           vector of timestamps at which history is recorded
%   Xh          matrix of with values of (X, Y, Z) at time t1 < t < t2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% initial conditions
X = XYZinit';

%% parameters
% Pr = 10;
% b = 8/3.;
% r = 0;

%% simulation control
tau = 0;
tauend = t2;
dtau = 1e-3;

%% graphical output
tout = t1;
dtout = 1e-2;

%% integration in time

% Runge-Kutta 3rd order
rk1 = [0 -5/9 -153/128];
rk2 = [1/3 15/16 8/15];
rk3 = [0 1/3 3/4];

Xh = [];
t = [];

while (tau <= tauend)

    g = zeros(3,1);
    for k = 1:3

        % right-hand sides
        g(1) = Pr*(X(2)-X(1));
        g(2) = -X(1)*X(3) + r*X(1) - X(2);
        g(3) = X(1)*X(2) - b*X(3);

        X = X + rk2(k)*dtau*g;

    end


    if (tau >=tout)
        tout = tout + dtout;

        Xh(:,end+1) = X;
        t(end+1) = tau;

        if (plotting)
            plot3(Xh(1,:),Xh(2,:),Xh(3,:))
            hold on
            plot3(X(1),X(2),X(3),'r.')
            hold off
            view(51,26)
            drawnow
        end
    end


    tau = tau + dtau;

end

end