clc
clear all
% **** System Inputs ****
par = v4_new_4d_parameters();
wG_in0 = 0.36;
wL_in0 = 8.64;
% **** Automatic ********
zs = [0.01:0.01:0.05 0.055 0.06:0.02:0.12 0.2:0.1:1];

ns0 = size(zs);
ns = ns0(2);
y1min = zeros(ns,1);
y1max = zeros(ns,1);
y1init = zeros(ns,1);
y1end = zeros(ns,1);

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

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

dx = [-0.1;0;0;0];
M=eye(4);
options=odeset('Mass',M,'BDF','on','AbsTol',1e-12,'RelTol',1e-12,'MaxStep',1);

for s = 1:ns
% **** System Inputs ****
z0 = zs(s);
disp(['Z=' num2str(z0)])
u0 = [z0;wG_in0;wL_in0];
% ******** Initializing *************
[x0,y0,par] =  v4_new_4d_initialize(u0,par);
% **********************

tf = 4;  
u=[zs(s);wG_in0;wL_in0];

%** Simulation using Simple method *****
[t,x]=ode15s(@v4_new_4d_model,[0;3600*tf],x0 + dx,options,u,'derivatives',par);
% *** Calculating measurements ****
n = length(t);
y = zeros(3,n);
for k = 1:n
    yt = v4_new_4d_model(t(k),x(k,:),u,'measurements',par);
    y(:,k) = yt(1:3);
end
n1 = floor(n/2);
y1min(s) = min(y(1,n1:n));
y1max(s) = max(y(1,n1:n));
y1init(s) = y(1,1);
y1end(s) = y(1,end);

y2min(s) = min(y(2,n1:n));
y2max(s) = max(y(2,n1:n));
y2init(s) = y(2,1);
y2end(s) = y(2,end);

y3min(s) = min(y(3,n1:n));
y3max(s) = max(y(3,n1:n));
y3init(s) = y(3,1);
y3end(s) = y(3,end);
clear y x t;
end

save('v01_bif_data.mat','y1min','y1max','y1init','y1end','y2min','y2max','y2init','y2end','y3min','y3max','y3init','y3end');
Steady_State_auto;