% Run RELAXATION algorithm.
%
% Based on:
%   Trimborn, Koch, Steger (2007) Multidimensional Transitional Dynamics: A
%   Simple Numerical Procedure, forthcoming in Macroeconomic Dynamics
%   Copyright by Trimborn, Koch, Steger, 2008
%   visit http://www.relaxation.uni-siegen.de

%clear all
%clear global
disp('Initialize Relaxation algorithm');

tic

%globalpar                           % Initializes the global parameter
% done in script instead
%parini                              % Loads the Parameter Values
% done in script instead

relaxsetting                        % Loads the settings for the Relaxation algorithm, i.e. dimensions, boundary conditions etc. 

% Converts settings to a form suitable for relax.m
[guess, start, errorcode]=initrelax(@funcODE, @funcSTAT, n, n1, n3, nu, y, M, statev);       


if errorcode==0            % Executes the relaxation algorithm if no error occured during initilization
    
    [t, x]=relax(@funcODE, @funcSTAT, @funcINI, @funcfinal, n, n1, n3, nu, guess, M, start, Endcond, maxit, tol, damp, dampfac);    
    
    %Normalization of specified variables
    for i=1:M
        x(normal,i)=x(normal,i)./x(normal,end);
    end
    
    varex                               % Extracts the variables and stores them in the memory
 
%    If you want to calculate and display the eigenvalues at the steady
%    state, remove the comment of the two subsequent lines.
%    [EVa EVe Jac]=eigDAS(@funcODE, @funcSTAT, x(:,end));
%    disp(['Eigenvalues: ',num2str(EVa')]);disp([' ']);
    
end


%Calculation time
time=toc;
timesec=mod(time,60);
timemin=floor(time/60);
disp(['Calculation time: ',num2str(time),' seconds (',num2str(timemin),' min ',num2str(timesec),' sec)']);

%To get a first impression of the results, remove the comment of the
%subsequent line
%plotrelax(t, x, n1, 100)

Tfinal=100;

%Extract equilibrium paths
k_path=x(1,:);
c_path=x(2,:);
m_path=x(3,:);
