clc
clear all
close all
% **** System Inputs ****
par = v1_new_6d_parameters();
% **** Automatic ********
zs = [0.03:0.01:0.1 0.12 0.15 0.2 0.3:0.1:1];
ns0 = size(zs);
ns = ns0(2);
y1min = zeros(ns,1);
y1max = zeros(ns,1);

y2min = zeros(ns,1);
y2max = zeros(ns,1);

y3min = zeros(ns,1);
y3max = zeros(ns,1);

y4min = zeros(ns,1);
y4max = zeros(ns,1);

y5min = zeros(ns,1);
y5max = zeros(ns,1);

y6min = zeros(ns,1);
y6max = zeros(ns,1);

y7min = zeros(ns,1);
y7max = zeros(ns,1);

y8min = zeros(ns,1);
y8max = zeros(ns,1);

dx = [0;0;-0.1;0;0;0];
z1 = 1;
M=eye(6);
options=odeset('Mass',M,'BDF','on','AbsTol',1e-12,'RelTol',1e-12,'MaxStep',1);
tf = 6;
for s = 1:ns
    % **** System Inputs ****
    z2 = zs(s);
    u0 = [z1;z2;0;0];
    % ******** Initializing *************
    [x0,y0,par] =  v1_new_6d_initialize(u0,par);
    %[A,B,C,D]=v2_new_6d_linmod(x0,y0,u0,par);
    [A,B,C,D] = v1_new_6d_linmodl_num(x0,u0,par); %numerical
    % **********************
    [eigval]=real(eig(A));
    unstable=0;
    for i=1:length(eigval)
        if eigval(i) >=0
            unstable=1;
        end
    end
    if(unstable) %If unstable system
        disp(['Z2=' num2str(z2)])
        disp('Unstable')
        %** Simulation using Simple method *****
        [t,x]=ode15s(@v1_new_6d_model,[0;3600*tf],x0 + dx,options,u0,'derivatives',par);
        % *** Calculating measurements ****
        n = length(t);
        y = zeros(8,n);
        for k = 1:n
            yt = v1_new_6d_model(t(k),x(k,:),u0,'measurements',par);
            y(:,k) = yt([1 3 4 7 9 8 2 5]);
        end
        n1 = floor(n/2);
        y1min(s) = min(y(1,n1:n));
        y1max(s) = max(y(1,n1:n));
        
        y2min(s) = min(y(2,n1:n));
        y2max(s) = max(y(2,n1:n));
        
        y3min(s) = min(y(3,n1:n));
        y3max(s) = max(y(3,n1:n));
        
        y4min(s) = min(y(4,n1:n));
        y4max(s) = max(y(4,n1:n));
        
        y5min(s) = min(y(5,n1:n));
        y5max(s) = max(y(5,n1:n));
        
        y6min(s) = min(y(6,n1:n));
        y6max(s) = max(y(6,n1:n));
        
        y7min(s) = min(y(7,n1:n));
        y7max(s) = max(y(7,n1:n));
        
        y8min(s) = min(y(8,n1:n));
        y8max(s) = max(y(8,n1:n));
    else %If stable system
        disp(['Z2=' num2str(z2)])
        disp('Stable')
        y1min(s) = y0(1);
        y1max(s) = y0(1);
        
        y2min(s) = y0(3);
        y2max(s) = y0(3);
        
        y3min(s) = y0(4);
        y3max(s) = y0(4);
        
        y4min(s) = y0(7);
        y4max(s) = y0(7);
        
        y5min(s) = y0(9);
        y5max(s) = y0(9);
        
        y6min(s) = y0(8);
        y6max(s) = y0(8);
        
        y7min(s) = y0(2);
        y7max(s) = y0(2);
        
        y8min(s) = y0(5);
        y8max(s) = y0(5);
    end
    clear y x t;
end

save('v1_6d_bif_data.mat','zs','y1min','y1max','y2min','y2max','y3min','y3max','y4min','y4max','y5min','y5max','y6min','y6max','y7min','y7max','y8min','y8max');
Steady_State_auto;