diff --git a/ERD_boundary1.m b/Boundary1.m
similarity index 100%
rename from ERD_boundary1.m
rename to Boundary1.m
diff --git a/ERD_boundary2.m b/Boundary2.m
similarity index 100%
rename from ERD_boundary2.m
rename to Boundary2.m
diff --git a/ERD_boundary3.m b/Boundary3.m
similarity index 100%
rename from ERD_boundary3.m
rename to Boundary3.m
diff --git a/ERD_boundary4.m b/Boundary4.m
similarity index 100%
rename from ERD_boundary4.m
rename to Boundary4.m
diff --git a/DATA_case2.mat b/DATA_case2.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f14208f56b7a771b6bfd8753081021c22ffb5451
Binary files /dev/null and b/DATA_case2.mat differ
diff --git a/ERD_model.m b/ODEsystem.m
similarity index 100%
rename from ERD_model.m
rename to ODEsystem.m
diff --git a/Sim_0a_data.m b/Sim_0a_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..bc0f0792704d8a7b0260a1943906df7f3da4a579
--- /dev/null
+++ b/Sim_0a_data.m
@@ -0,0 +1,125 @@
+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0a_data(input1, input2)
+%%  Sim_0a_data(input)
+%
+%   Data for Simulation in chapter 3
+%
+%   SWRO-PRO hybrid system I (realistic ERD)
+%
+%   option_data = 0
+%
+%   Input:
+%       input1        -   x
+%       input2        -   x
+
+%% model versions
+version=zeros(1,10);
+% version(1)=0 if co-current, 1 otherwise
+version(2)=1;  % 0 = SWRO beta fixed
+version(3)=0;  % 0 = PRO beta fixed
+version(4)=1;  % 0 = ideal SWRO
+version(5)=1;  % 0 = ideal PRO
+% ERD
+version(6)=3;
+% version(6) = 0 --> only SWRO (no ERD)
+% version(6) = 1 --> only SWRO (with ERD)
+% version(6) = 2 --> only PRO
+% version(6) = 3 --> SWRO-PRO hybrid system (with one ERD)
+% version(6) = 4 --> SWRO-PRO hybrid system (with two ERDs) 
+version(7)=1; % 1 = ICP and ECP for SWRO
+version(8)=1; % 1 = ICP and ECP for PRO
+
+%% Membrane unit properties
+H = 1e-3;           % height of the membrane [m]
+ro_water = 1000;    % mass density of water  [kg/m^3] 
+ro_salt = 2165;     % mass density of salt [kg/m^3] 
+Mw = 18;            % molecular weight of water [kg/kmol]
+Ms = 58.44;         % molecular weight of salt  [kg/kmol]
+Rw = 462;           % gas constant of water  [J/(kg K)] 
+T0 = 297;           % temperature [K] 
+eta = 1.3e-3;       % seawater viscosity [kg/(m s)]
+p_r = 1e5;          % pressure [Pa]=[kg/ms^2]
+rho_r = 1e3;        % density [kg/m^3]
+C_r = 1;            % salt concentration [%]
+
+%% SWRO
+swro_Z=1;               % width of SWRO membrane [m]
+swro_L=4;               % length of SWRO membrane [m]
+if version(6)== 2; swro_L=1; end
+swro_alpha = 5.0815e-9; % SWRO water permeablity co-efficient [s/m]
+swro_R =0.96;           % SWRO salt rejection rate 
+swro_KK = 1e-2;         % SWRO ICP mass transfer coefficient
+swro_KD = 1/swro_KK;    % SWRO ECP draw side mass transfer coefficient
+swro_KF = 1/swro_KK;    % SWRO ECP fresh side mass transfer coefficient
+swro_x_r=swro_L^2;      % x_r=swro_L^2 since x = linspace(0,1,n) (if x = linspace(0,swro_L,n) then x_r=swro_L;)     
+swro_b1 = H/swro_x_r;                           % H/swro_L ratio
+swro_b2 = swro_Z/swro_x_r;                      % Z/swro_L ratio
+J_r = sqrt(H^3/swro_x_r*p_r*rho_r);             % flux [kg/s^2]
+swro_gamma= swro_x_r * p_r *  swro_alpha /J_r;  % SWRO scaling factor - mass balance
+swro_gamma2= J_r^2./(swro_x_r^2 * p_r * rho_r); % SWRO scaling factor - momentum balance
+swro_W_r=J_r*p_r/rho_r;                         % net work [W/m^2]
+sigma=0.999 ;                                   % Rejection coefficient
+swro_beta_fix=4.43e-4/J_r*swro_x_r;             % value for fixed SWRO beta [kg/sm^2]
+ 
+%% PRO
+Z=1;                    % width of the PRO membrane [m]
+L=1.5;                  % length of the PRO membrane [m]      
+alpha = 5.47e-9;        % water permeablity co-efficient [s/m]
+R =0.94;                % salt rejection rate [1]
+KK=7.13e2;              % mass transfer coefficient [sm^2/kg]
+KD=1/KK;                % PRO ECP draw side mass transfer coefficient [sm^2/kg]
+KF=1/KK;                % PRO ECP fresh side mass transfer coefficient [sm^2/kg]
+x_r=L^2;                % x_r=L^2 since x = linspace(0,1,n) (if x = linspace(0,L,n) then x_r=L;)              
+b1 = H/x_r;             % H/L ratio
+b2 = Z/x_r;                           % Z/L ratio
+Q_r = sqrt(H^3/x_r*p_r*rho_r);        % flux [kg/s^2]
+gamma= x_r * p_r * alpha /Q_r ;       % PRO scaling factor - mass balance
+gamma2= Q_r^2./(x_r^2 * p_r * rho_r); % PRO scaling factor - momentum balance
+W_r=Q_r*p_r/rho_r;                    % net work [W/m^2]
+beta_fix =1.71e-4/Q_r*x_r;            % value for fixed PRO beta [kg/sm^2]
+
+%% Sea Water
+cE= 35/983/C_r;                            % salt concentration in seawater
+pE= 1e5/p_r;                               % external pressure
+rho_E=(cE + 1)./(ro_water*cE/ro_salt + 1); % density of incomming seawater
+
+%% SWRO operating conditions
+J_sf_0 = 0;              % salt flux in fresh side at 0 
+J_wf_0 = 0;              % water flux in fresh side at L               
+Pf_L = pE;               % pressure fresh side at L
+
+Pd_0 =  66e5/p_r;        % pressure draw side at 0
+Pd_L =  60e5/p_r;        % pressure draw side at L (not need in the hybrid system)
+
+%% PRO operation conditions 
+Q_sf_0 = 0;              % salt flux in fresh side at 0
+pf_L = pE;               % pressure of fresh side at L
+
+pd_0 = 12e5/p_r;         % pressure draw side at 0
+pd_L = 14.6e5/p_r;       % pressure of fresh side at 0
+pf_0 = 1.3e5/p_r;        % pressure draw side at L 
+
+%% ERD/Turbine/Pump parameters
+T_eff  = .95;               % turbine efficiency
+HP_eff = .9;            	% high pressure pump efficiency
+LP_eff = .95;               % low pressure pump efficiency
+V_m=0.052;                  % Volumetric mixing
+ERD_eff=.96;                % ERD unit pressure efficiency
+ERD_fric=5e-04;             % ERD friction coefficient
+A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
+eta_ERD=0.01;               % leak of high pressure brine
+mix_density=997/rho_r;  	% density of mixture in ERD
+pw=2;                       % water price [$/m^3]
+pe=0.5;                     % electricity price [$/kWh]
+mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
+
+%% display figures
+fig=[1,0,0,0,0,0,0,0,1,0,0]; % f(i)=1 --> figure i will be displayed
+%fig=[0,0,0,0,0,0,0,0,0,0,0]; % f(i)=1 --> figure i will be displayed
+
+%% model specific changes:
+% co-current
+if pd_0 > pd_L; version(1)=0; else; version(1)=1; end            
+% no SWRO needed (--> use trivial data)
+if version(6) == 2; version(4)=0; Pd_0=40; Pd_L=35; end
+
+end
\ No newline at end of file
diff --git a/Sim_0_data.m b/Sim_0b_data.m
similarity index 97%
rename from Sim_0_data.m
rename to Sim_0b_data.m
index 4e8caebe9b26957995d19dadf0c99af049d0d6f9..2cd106e8cf641985104f9132dcc9b0cf6755c6bb 100644
--- a/Sim_0_data.m
+++ b/Sim_0b_data.m
@@ -1,11 +1,11 @@
-function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0_data(input1, input2)
-%%  Sim_32_data(input)
+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0b_data(input1, input2)
+%%  Sim_0b_data(input)
 %
 %   Data for Simulation in chapter 3
 %
 %   SWRO-PRO hybrid system I (realistic ERD)
 %
-%   option_data = 0
+%   option_data = 0.1
 %
 %   Input:
 %       input1        -   x
@@ -104,7 +104,7 @@ HP_eff = .9;            	% high pressure pump efficiency
 LP_eff = .95;               % low pressure pump efficiency
 V_m=0.052;                  % Volumetric mixing
 ERD_eff=.96;                % ERD unit pressure efficiency
-ERD_fric=5e-06;             % ERD friction coefficient
+ERD_fric=5e-04;             % ERD friction coefficient
 A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
 eta_ERD=0.01;               % leak of high pressure brine
 mix_density=997/rho_r;  	% density of mixture in ERD
diff --git a/Skript.asv b/Skript.asv
new file mode 100644
index 0000000000000000000000000000000000000000..e92418fdd16e09a87cd0f886271eab02f3ba8c20
--- /dev/null
+++ b/Skript.asv
@@ -0,0 +1,166 @@
+%% Matlab skript for the Paper 2024:
+%
+% - reproduce the simulations
+% - reproduce optimisation results
+% - reproduce all figures.png (converted to .eps using LibreOfficeDraw)
+% - test average computation time
+%
+
+%% Numerical Simulations
+%
+% Create figures form chapter: Numerical simulations
+% 
+close all;
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4),
+%
+%figure(1); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_6.png');
+%figure(2); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_7.png');
+%figure(3); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_8_ERD.png');
+%
+
+%% Average simulation time:
+%
+% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data"
+%
+T=zeros(1,100); 
+for i=1:length(T)
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4);
+T(i)=time;
+end
+mean(T) %= 0.9600
+%
+
+%% Optimisation - case2 (with epsilon constraints)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case2.mat")
+% 
+% load("DATA_case2.mat")
+% save DATA_case2.mat X_2
+% 
+% X2  -  set of operating conditions for epsilon constraint method (2x40)
+% Y2  -  system output using X2 (4x40)
+%
+% X2_P  -  set of operating conditions from paretosearch (2x200)
+% Y2_P  -  system output using X2_P (4x200)
+%
+% time2    - time for epsilon constraint method (40points)
+% time2_P  - time for paretosearch solver (200 points)
+%
+
+X2=zeros(2,40);Y2=zeros(4,40);X2_P=zeros(2,200);
+
+
+
+load("DATA_T1.mat")
+c1=datetime("now");
+x0=[69.9 51.5];
+FWmin=linspace(0.3770, 0.0015, 40);
+for i=1:40
+clc,
+i, 
+close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
+option_mesh = 1e4; option_BVP = 1e-4; option_data = 1;
+[x_opt, f_opt]=fmincon(@(x)fun_1(x,1,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',1,option_mesh,option_BVP, FWmin(i)),foptions);
+end
+c2=datetime("now"); %von 2 bis 39
+save DATA_T1.mat c1 c2 % 00:02:30
+
+%% Optimisation - case2 (with paretosearch)
+%
+% Create Paretofront for SWRO+ERD 
+%
+
+
+
+%% Optimisation-2 (red)
+close all;clear all;clc
+c3=datetime("now");
+x0=[69.996767760342830 49.753887785404920];
+FWmin=linspace(0.3862, 0.0066, 40);
+for i=2:39
+clc,
+i,
+close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
+option_mesh = 1e4; option_BVP = 1e-4; option_data = 2;
+[x_opt, f_opt]=fmincon(@(x)fun_1(x,2,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',2,option_mesh,option_BVP, FWmin(i)),foptions);
+end
+c4=datetime("now"); %von 2 bis 39
+save DATA_T2.mat c3 c4
+(c4-c3)/38, %00:01:54
+
+%% hybrid paretosearch - old
+close all; clear all; clc
+rng default % For reproducibility
+A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
+% initial points:
+load('DATA_Approx.mat');
+X0 = X_approx'; %size(X0)= (93,5)
+%remove infeasible poitns
+z=length(X_approx(1,:));
+for i=1:z
+    k=0;
+    for j=1:4
+    if X_approx(j,z+1-i) > ub(j); k=1;end
+    if X_approx(j,z+1-i) < lb(j); k=1;end
+    end
+    if k==1;X0(z+1-i,:)=[]; end
+end
+%After removing points:
+%size(X0)= (61,5)
+%
+options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
+% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
+[X_paretosearch_3, Y_paretosearch_3] = paretosearch(@(x)fun_1(x,3,'Pareto',1e4,1e-4),5,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
+
+%%
+load DATA_Pareto.mat
+scatter(SEC_Pareto_1,FW_Pareto_1);hold on
+scatter(SEC_Pareto_2,FW_Pareto_2)
+
+%% red - Paretosearch
+%
+close all; clear all; clc
+rng default % For reproducibility
+%
+A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+%A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
+%
+% initial points:
+load DATA_Pareto.mat
+X0 = [X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2'];
+%
+options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
+c5=datetime("now");
+% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
+[X, Y] = paretosearch(@(x)fun_1(x,2,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
+c6=datetime("now");
+scatter(-Y(:,1),-Y(:,2));hold on
+scatter(SEC_Pareto_2,FW_Pareto_2)
+save DATA_T3 c5 c6
+%c6-c5 = 00:56:57
+
+%% blue - Paretosearch
+%
+close all; clear all; clc
+rng default % For reproducibility
+%
+A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
+%A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
+%
+% initial points:
+load DATA_Pareto.mat
+X0 = [X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1'];
+%
+options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
+c7=datetime("now");
+% x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
+[X, Y] = paretosearch(@(x)fun_1(x,1,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
+c8=datetime("now");
+scatter(-Y(:,1),-Y(:,2));hold on
+scatter(SEC_Pareto_1,FW_Pareto_1)
+save DATA_T4 c7 c8
+%c8-c7 = 01:05:44
diff --git a/Skript.m b/Skript.m
index dd695017ee29c3e105d3ce19b3a7d75f43551d72..2a11bef1161e88f8b43d32db22e8f65b14d1b996 100644
--- a/Skript.m
+++ b/Skript.m
@@ -20,21 +20,40 @@ close all;
 
 %% Average simulation time:
 %
-% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data"
+% calculate average computation time of previous simulation
 %
 T=zeros(1,100); 
 for i=1:length(T)
-[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4);
+[Output1, Output2, time]=fun_1(1,0.1,'sol',1e4,1e-4);
 T(i)=time;
 end
 mean(T) %= 0.9600
 %
 
-%% Optimisation-1(blue)
+%% Optimisation - case2 (with epsilon constraints)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case2.mat")
+% 
+% load("DATA_case2.mat")
+% save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P
+% 
+% X2  -  set of operating conditions for epsilon constraint method (2x40)
+% Y2  -  system output using X2 (4x40)
+%
+% X2_P  -  set of operating conditions from paretosearch (2x200)
+% Y2_P  -  system output using X2_P (4x200)
+%
+% time2    - time for epsilon constraint method (40points)
+% time2_P  - time for paretosearch solver (200 points)
+%
+
+load("DATA_T1.mat")
 c1=datetime("now");
 x0=[69.9 51.5];
 FWmin=linspace(0.3770, 0.0015, 40);
-for i=2:39
+for i=1:40
 clc,
 i, 
 close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
@@ -45,6 +64,13 @@ end
 c2=datetime("now"); %von 2 bis 39
 save DATA_T1.mat c1 c2 % 00:02:30
 
+%% Optimisation - case2 (with paretosearch)
+%
+% Create Paretofront for SWRO+ERD 
+%
+
+
+
 %% Optimisation-2 (red)
 close all;clear all;clc
 c3=datetime("now");
diff --git a/fun_1.asv b/fun_1.asv
new file mode 100644
index 0000000000000000000000000000000000000000..15e4f1a08dbbad837c03eb4dc451cdf7521090fa
--- /dev/null
+++ b/fun_1.asv
@@ -0,0 +1,836 @@
+function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,option_BVP)
+%%  fun   Solves BVP for SWRO-PRO hybrid system using bvp5c
+%
+%       fun(input,option_data,obj,option_mesh,option_BVP) will solve the
+%       nonlinear BVP of the selected SWRO-PRO system. In "option_data" 
+%       all Informations of the system are given with exception of the 
+%       additional "input". 
+%
+%   Input:
+%       input         -   Input for the Data function
+%       option_data   -   Selects Data function used: data(input)
+%       obj           -   'SEC', 'FW', 'Rev', 'PD_net' , 'SE_net', 'approx'  or 'sol' ;
+%       option_mesh   -   NMax of bvp5c
+%       option_BVP    -   RelTol of bvp5c
+%
+%   Output:
+%       output1 = [SEC_net, FW, Rev, SWRO_Rec, PRO_Rec, PD_net, SE_net];
+%       output2 = [W_net, -W_p1, -W_p2,  -W_p3, -W_p4, W_t];
+%       figures may be displayed if obj='sol'
+%
+%   Example:
+%       x0=1;       
+%       D=.11;      
+%       obj='sol';  
+%       mesh=2e4;
+%       tol=1e-6;  
+%       [Sim_01_output1, Sim_01_output2]=fun(x0,D,obj,mesh,tol),
+%
+tic;
+%% Read data
+if option_data == 0; DATA = @(x)Sim_0a_data(input); end
+if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end
+
+
+
+% Simulations
+if option_data == 0.11; DATA = @(x)Sim_11_data(input); end
+if option_data == 0.12; DATA = @(x)Sim_12_data(input); end
+if option_data == 0.13; DATA = @(x)Sim_13_data(input); end
+if option_data == 0.21; DATA = @(x)Sim_21_data(input); end
+if option_data == 0.22; DATA = @(x)Sim_22_data(input); end
+if option_data == 0.23; DATA = @(x)Sim_23_data(input); end
+if option_data == 0.24; DATA = @(x)Sim_24_data(input); end
+if option_data == 0.31; DATA = @(x)Sim_31_data(input); end
+if option_data == 0.32; DATA = @(x)Sim_32_data(input); end
+% PRO Membrane Comparison
+if option_data == 0.51; DATA = @(x)Sim_PRO_01_data(input); end
+if option_data == 0.52; DATA = @(x)Sim_PRO_02_data(input); end
+if option_data == 0.53; DATA = @(x)Sim_PRO_03_data(input); end
+% Rough vs. Detailled
+if option_data == 0.91; DATA = @(x)Sim_rough_approx_data(input); end
+% Optimizations
+if option_data==1; DATA = @(x)Opt_1_data(input); end
+if option_data==2; DATA = @(x)Opt_2_data(input); end
+if option_data==3; DATA = @(x)Opt_3_data(input); end
+if option_data==4; DATA = @(x)Opt_4_data(input); end
+if option_data==5; DATA = @(x)Opt_L_data(input); end
+% for 3d plot
+if option_data==11; DATA = @(x)plot3d_data(input); end
+
+[H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD, KF, KD] ...
+    = DATA(1);
+
+%% Set options
+ode_options = bvpset('Stats','off','NMax', option_mesh, 'RelTol', option_BVP);      
+
+try
+n=25;
+x = linspace(0,1,n);    
+  
+%% Create initial guess for the BVP
+J_wd_0 = (983/1018)/J_r;          % guess for J_wd(0)
+J_wf_0 = 0.01353/J_r;             % should be close to 0
+Pf_0   = 1.1*Pf_L;                % guess for P_f(0)
+Q_wd_0 = 0.01353*(983/1018)/Q_r;  % guess for Q_wd(0) 
+Q_wf_0 = 0.01353/Q_r;             % guess for Q_wf(0)
+
+if version(1)==0; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE;  Q_wd_0; Q_sf_0; Q_wf_0; pd_0; pf_0]; end
+if version(1)==1; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE; -Q_wd_0; Q_sf_0; Q_wf_0; pd_0; pf_0]; end
+solinit = bvpinit(x,y_init);
+
+%% Solve the BVP
+if version(6) == 0; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 1; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary1(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 2; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 3; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary3(ya,yb, DATA),solinit,ode_options); end       
+if version(6) == 4; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary4(ya,yb, DATA),solinit,ode_options); end   
+            
+%% Evaluate the solution
+n=100;
+x=linspace(0,1,n);
+y = deval(sol,x); Y = y'; Y = real(Y);
+
+%% SWRO evaluation
+C_d = Y(:,1);
+C_f = Y(:,3)./Y(:,4);
+J_sd = Y(:,1).*Y(:,2);
+J_wd = Y(:,2); 
+J_d = Y(:,1).*Y(:,2)+Y(:,2);
+J_sf = Y(:,3);
+J_wf = Y(:,4);
+J_f = Y(:,3)+Y(:,4);
+P_d = Y(:,5); 
+P_f = Y(:,6); 
+% same quantities as in "ODEsystem.m" 
+swro_p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1)/Ms)/p_r;
+swro_p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3)/Ms./Y(:,4))/p_r;
+swro_del_c = Y(:,1)./(Y(:,1) + 1) - Y(:,3)./(Y(:,3) + Y(:,4));
+swro_local_ro_d = (Y(:,1) + 1)./(ro_water*Y(:,1)/ro_salt + 1);
+swro_local_ro_f = (Y(:,3) + Y(:,4))./(ro_water*Y(:,3)/ro_salt + Y(:,4));
+swro_beta = (1 - swro_R)*((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))/swro_R;
+if version(2) == 0 ;swro_beta= swro_beta_fix.*ones(1,n); end
+if version(4) == 0 ;swro_beta=zeros(1,n); end
+% J_w,in and J_s,in
+J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f));
+if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f));  end 
+J_sin=swro_beta.*swro_del_c;
+if version(7)==1
+J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) ));
+J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD +  Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD));
+end
+
+%% PRO evaluation
+c_d = Y(:,1+6);                    % concentration of draw side
+c_f = Y(:,3+6)./Y(:,4+6);          % concentration of fresh side
+Q_sd = Y(:,1+6).*Y(:,2+6);          % salt flow rate in draw side
+Q_wd = Y(:,2+6);                   % water flow rate in draw side
+Q_d = Y(:,1+6).*Y(:,2+6)+Y(:,2+6); % total flow rate draw side
+Q_sf = Y(:,3+6);                   % salt flow rate in fresh side
+Q_wf = Y(:,4+6);                   % water flow rate in fresh side
+Q_f = Y(:,3)+Y(:,4+6);            % total flow rate fresh side
+p_d = Y(:,5+6);                     % Pressure draw side
+p_f = Y(:,6+6);                     % Pressure fresh side
+% same quantities as in "ODEsystem.m" 
+p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1+6)/Ms)/p_r;
+p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3+6)/Ms./Y(:,4+6))/p_r;
+del_c = Y(:,1+6)./(Y(:,1+6) + 1) - Y(:,3+6)./(Y(:,3+6) + Y(:,4+6));
+local_ro_d = (Y(:,1+6) + 1)./(ro_water*Y(:,1+6)/ro_salt + 1);            	
+local_ro_f = (Y(:,3+6) + Y(:,4+6))./(ro_water*Y(:,3+6)/ro_salt + Y(:,4+6)); 
+beta = (1 - R)*((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)))./R;
+if version(3) == 0 ; beta =  beta_fix.*ones(1,n); end
+if version(5) == 0 ; beta = zeros(1,n); end
+% Q_w,in and Q_s,in
+Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))');
+if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end
+Q_sin = beta'.*del_c;
+if version(8)==1 
+Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))';
+Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD +  Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta');
+end
+    
+%% Freswater production and recovery rates
+Freshwater = J_wf(end)*J_r; %in [kg/s] 
+FW = (Freshwater*swro_Z/(swro_local_ro_f(end)*rho_r))*3600; %in [m^3/h] 
+SWRO_Recovery =(J_f(end)./J_d(1))*100;     % SWRO freshwater recovery [%]
+PRO_Recovery = (1- Q_f(end)./Q_f(1))*100;  % PRO recovery rate [%]
+% warning flags:
+    if FW < 0 
+        fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n');
+    end
+    if SWRO_Recovery > 100 
+        fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n');
+    end
+    if SWRO_Recovery < 0 
+        fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
+    end
+    if PRO_Recovery > 100 
+        fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
+    end
+    if PRO_Recovery < 0 
+        fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+    end
+
+%% Output of System
+J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1);
+W_p1 = 1/HP_eff * (P_d(1)-pE)*(J_E *swro_Z)/rho_E; W_p1 = W_p1*swro_W_r;
+W_p2 = 0; W_p2 = W_p2*swro_W_r; 
+W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
+
+    %% Version(6)=0 -->  only SWRO (no ERD)
+    if version(6)==0; W_p3=W_p1; W_p4=0; W_t=0; end
+
+    %% Version(6)=1 -->  only SWRO (with ERD)
+    if version(6)==1
+        rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E;
+        C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);
+        J_wERD = J_E/(C_ERD+1);
+        Q_exit = J_d(end)*(1-eta_ERD);
+        c_exit = (C_d(end)*J_wd(end)*(1 - eta_ERD)+cE*J_wE-C_ERD*J_wERD)/(J_wd(end)*(1 - eta_ERD)+J_wE-J_wERD);
+        rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE;
+        f_1= (Q_exit*swro_Z/rho_exit*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 );
+        f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); 
+        f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+        pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0;
+            if W_p3 < 0 
+                fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
+            end
+
+    end
+    %% Version(6)=2 --> only PRO
+    if version(6)==2
+        if version(1) ==0; W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); end
+        if version(1) ==1; W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1); end
+        W_p1=0; W_p3=0; W_t=W_t*W_r;
+        if version(1) ==0; W_p2= 1/LP_eff * (p_d(1)-pE)*(Q_d(1)*Z)./local_ro_d(1); end 
+        if version(1) ==1; W_p2= 1/LP_eff * (p_d(end)-pE)*(abs(Q_d(end))*Z)./local_ro_d(end); end 
+        W_p2 = W_p2*W_r;
+    end
+
+    %% Version(6)=3 --> SWRO-PRO hybrid system (with one ERD)
+    if version(6)==3
+        rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E;
+        C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);  
+        qq=Q_r/J_r;
+            if version(1) == 0 %co-current
+                f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2;  
+                f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+                W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
+                W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); 
+            end
+            if version (1) == 1 %counter-current   
+                f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2  ;
+                f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;     
+                W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
+                W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1);
+            end
+        W_p3=W_p3*swro_W_r; W_t=W_t*W_r;
+    end
+
+    %% Version(6)=4 --> SWRO-PRO hybrid system (with two ERDs)
+    if version(6)==4 
+        if version(1) == 0; rho_ERD2= V_m*(local_ro_d(end) - rho_E)+rho_E; end
+        if version(1) == 1; rho_ERD2= V_m*(local_ro_d(1) - rho_E)+rho_E; end
+        C_ERD2= -ro_salt*(rho_ERD2-1)/(rho_ERD2*ro_water-ro_salt);
+        rho_ERD= V_m*(swro_local_ro_d(end) - rho_ERD2)+rho_ERD2;
+        C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);
+        J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; %J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1);
+        qq=Q_r/J_r;
+        J_E_sea = 2*J_E;
+        J_wE_sea = J_E_sea/(cE+1);
+        if version(1) == 0
+            Q_exit = mixer_ERD*Q_d(end)*(1-eta_ERD);
+            c_exit = (c_d(end)*mixer_ERD*Q_wd(end)*(1 - eta_ERD)+cE*J_wE_sea-C_ERD2*2*J_E)/(mixer_ERD*Q_wd(end)*(1 - eta_ERD)+J_wE_sea-2*J_E);
+        end
+        if version(1) == 1
+            Q_exit = mixer_ERD*abs(Q_d(1))*(1-eta_ERD);
+            c_exit = (c_d(1)*mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+cE*J_wE_sea-C_ERD2*2*J_E)/(mixer_ERD*abs(Q_wd(1))*(1 - eta_ERD)+J_wE_sea-2*J_E);
+        end
+
+
+%% Removed comments
+
+%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here??
+
+% -3.9191    0.3247    0.0131   32.8727    7.1201       NaN       NaN
+%
+% 1.0e+03 *
+% -1.2742   -0.9744         0   -0.4873   -0.0010    0.1885
+
+% In ideal cases take mean of some quantities. 
+% In ideal cases salt will fluctuate around 1e-16. We take mean here for 
+% 
+% This is just for figures. (), otherwise one could use ylim ) 
+%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end
+%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end
+
+        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
+        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
+        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
+
+        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
+
+        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;       
+        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
+        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end     
+
+        %J_wERD = J_E/(C_ERD+1);
+
+
+
+
+
+
+
+
+
+
+
+        rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE;
+        %calculation of Pressure P_ERD2
+            if version(1) == 0 %co-current
+                f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
+                f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2;  
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
+                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
+                if pERD2 > pMax; pERD2=pMax; end
+                %if pERD2 < 0; pERD2=pMax; end  
+                
+            end
+            if version (1) == 1 %counter-current 
+                f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
+                f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2;  
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
+                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
+                if pERD2 > pMax; pERD2=pMax; end
+                %if pERD2 < 0; pERD2=pMax; end  
+            end     
+            %calculation of Pressure P_ERD
+            if version(1) == 0 %co-current
+                f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2;  
+                f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+                f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+                if pERD > pMax; pERD=pMax; end
+                %if pERD < 0; pERD=pMax; end  
+                W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
+                W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); 
+            end
+            if version (1) == 1 %counter-current   
+                f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2  ;
+                f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
+                f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
+                if pERD > pMax; pERD=pMax; end 
+                %if pERD < 0; pERD=pMax; end  
+                W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
+                W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1);
+            end         
+    % Recalculation for Pump 1
+    W_p1 = 1/HP_eff * (P_d(1)-pERD2)*(J_E *swro_Z)/rho_ERD2; W_p1 = W_p1*swro_W_r;
+    end
+    
+%% final output:
+W_net= - W_p1 - W_p2 - W_p3 - W_p4 + W_t; % in [W]
+PD_net = W_net./(Z*L); %in [W/m^2]    
+%old SE_net
+%if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + Q_d(1)*Q_r*Z/(rho_r*local_ro_d(1)))/1000/3600; end %in [kWh/m^3]
+%if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)) + abs(Q_d(end))*Q_r*Z/(rho_r*local_ro_d(end)))/1000/3600; end %in [kWh/m^3]
+%new SE_tot
+if version(1)==0; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
+if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600; end %in [kWh/m^3]
+SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3]  
+Rev= pw*FW + pe*SEC_net*FW; %in [$/h]
+tt=toc;
+%% if BVP-solver failed
+catch
+    sol.stats.maxerr =1000;
+    SEC_net=1000;FW=0;Rev=-20;SWRO_Recovery=NaN;PRO_Recovery=NaN;PD_net=NaN;SE_net=NaN;
+    W_net=NaN;W_p1=NaN;W_p2=NaN;W_p3=NaN;W_p4=NaN;W_t=NaN;
+end
+%% Output of objective function
+switch (obj)
+        case 'SEC'
+            if sol.stats.maxerr > option_BVP
+            output1 = 20;
+            else; output1 = -SEC_net;
+            end 
+        case 'FW' 
+            if sol.stats.maxerr > option_BVP
+            output1 = 20;
+            else; output1 = -FW;
+            end 
+        case 'Rev' 
+            if sol.stats.maxerr > option_BVP
+            output1 = 0;
+            else; output1 = -Rev;
+            end
+        case 'Pareto' 
+            if sol.stats.maxerr > option_BVP    
+            output1 = [NaN, NaN];
+            else; output1 = [-SEC_net, -FW];
+                %if SWRO_Recovery < 30;output1 = [NaN, NaN];end
+                %if SWRO_Recovery > 60;output1 = [NaN, NaN];end
+                %if PRO_Recovery < 40;output1 = [NaN, NaN];end
+                %if PRO_Recovery > 95;output1 = [NaN, NaN];end
+            end 
+        case 'PD_net' 
+            if sol.stats.maxerr > option_BVP
+            output1 = 0;
+            else; output1 = -PD_net;
+            end
+        case 'SE_net' 
+            if sol.stats.maxerr > option_BVP
+            output1 = 0;
+            else; output1 = -SE_net;
+            end
+        case 'approx' 
+            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
+            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            %sol.stats.nmeshpoints
+            %sol.stats.maxerr
+            %sol.stats.nODEevals
+            %sol.stats.nBCevals
+            output2=[sol.stats.nmeshpoints sol.stats.maxerr];
+    case 'exe' 
+            if sol.stats.maxerr > option_BVP
+            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
+            output2=[NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN,NaN, NaN, NaN, NaN, NaN, NaN];
+            else
+            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
+            if version(6) == 3; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
+            if version(6) == 4; output1=[SEC_net, FW, SWRO_Recovery, PRO_Recovery]; end
+            output2=[x*swro_L; 100*C_d; J_r*J_d; p_r*P_d; p_r*P_f; x*L; 100*c_d; Q_r*Q_d; Q_r*Q_f; p_r*p_d; p_r*p_f]; % 11
+            output3=[x(2:end)*swro_L; 100*C_f(2:end); J_r*J_f(2:end); J_r/swro_x_r*J_cross(2:end); J_r/swro_x_r*J_sin(2:end); x(2:end)*L; 100*c_f(2:end); Q_r/x_r*Q_cross(2:end); Q_r/x_r*Q_sin(2:end)]; %9
+            %
+            output1=output1';output2=output2';output3=output3';
+            [a,b]=size(output1);
+            output1=[1:1:b; output1];
+            [a,b]=size(output2);
+            output2=[1:1:b; output2];
+            [a,b]=size(output3);
+            output3=[1:1:b; output3];
+            end
+    case 'sol' 
+            if sol.stats.maxerr > option_BVP
+            output1 = [NaN, NaN, NaN, NaN, NaN, NaN, NaN];
+            output2=[NaN, NaN, NaN, NaN, NaN, NaN];
+            output3=NaN;
+            else
+            if version(6) == 0; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 1; output1=[SEC_net, FW, Rev, SWRO_Recovery, NaN, NaN, NaN]; end
+            if version(6) == 2; output1=[NaN, NaN, NaN, NaN, PRO_Recovery, PD_net, SE_net]; end
+            if version(6) == 3; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            if version(6) == 4; output1=[SEC_net, FW, Rev, SWRO_Recovery, PRO_Recovery, NaN, NaN]; end
+            output2=[W_net, -W_p1, -W_p2,  -W_p3, -W_p4, W_t];
+            output3=tt;
+            end
+  
+            
+%% figure 1
+if sol.stats.maxerr ==1000
+fprintf(2,' \nERROR: bvp5c could not satisfy the relative error tolerance ---> no figures could be displayed\n');
+else
+lw=1.5; %Linewidth for all figures
+if fig(1) == 1
+f=figure(1); %
+f.Position = [276.3333 196.3333 1100 1000];
+x2=x(2:end);
+% total mass flows
+subplot(3,2,1);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% concentrations
+subplot(3,2,3);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
+yyaxis right; C_f=C_f(2:end);
+plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Pressures
+subplot(3,2,5);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% total mass flows
+subplot(3,2,2);lc='#0072bD';rc='#77AC30';
+plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% concentrations
+subplot(3,2,4);lc='#0072bD';rc='#77AC30';
+plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
+yyaxis right; c_f=c_f(2:end);
+plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','northwest');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
+% Pressures
+subplot(3,2,6);lc='#0072bD';rc='#77AC30';
+plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+f=figure(2); %
+f.Position = [199 232.3333 1.0987e+03 639.3333];
+% Permeate flows 
+subplot(2,2,1); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
+plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
+yyaxis right
+plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Osmotic/ Hydraulic pressure difference
+subplot(2,2,3);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
+plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
+yyaxis right; 
+plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Permeate flows 
+subplot(2,2,2); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
+plot(x2*L, -Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x [m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
+yyaxis right
+plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','southwest');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Osmotic/ Hydraulic pressure difference
+subplot(2,2,4);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
+plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
+yyaxis right; 
+plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','southeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
+
+
+
+end
+%% figure 2
+if fig(2) == 1
+figure(2); % pressure+concentration in both PRO channels
+subplot(2,2,1);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_d(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(2,2,2);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('p_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(2,2,3);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_d(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+x2=x(2:end);c_f2=100.*c_f(2:end);
+subplot(2,2,4);plot(x2*L, c_f2*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('c_f(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 3
+if fig(3) == 1
+figure(3); % Q_s,in and Q_w,in
+x2=x(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
+subplot(1,2,1);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{w,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(1,2,2);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('Q_{s,in}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 4
+lw=1.5; %Linewidth for all figures
+if fig(4) == 1
+figure(4); % Beta of PRO
+    if version(2)==0
+    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]);ylim([swro_beta(1)*J_r/swro_x_r*.9  swro_beta(1)*J_r/swro_x_r*1.1]); hold on
+    else
+    plot(x*swro_L, swro_beta*J_r/swro_x_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('B(x)','Location','northeast');xlim([0 swro_L]); hold on
+    end
+end
+%% figure 5
+if fig(5) == 1
+figure(5); % pressure+concentration in both SWRO channels
+subplot(2,2,1);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
+x2=x(2:end);C_f2=100.*C_f(2:end);
+subplot(2,2,2);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(2,2,4);plot(x2*swro_L, C_f2*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(2,2,3);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 6
+if fig(6) == 1
+figure(6); % J_s,in and J_w,in
+x2=x(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
+subplot(1,2,2);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(1,2,1);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 7
+if fig(7) == 1
+f=figure(7); % PRO quantities
+f.Position = [100 1 1100 1000];
+subplot(4,3,4);plot(x*L, Q_sd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,1);plot(x*L, Q_wd*Q_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,5);plot(x*L, Q_sf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,2);plot(x*L, Q_wf*Q_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0; 
+subplot(4,3,7);plot(x*L, p_d*p_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','Location','northeast'); xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(4,3,8);plot(x*L, p_f*p_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(4,3,10);plot(x*L, 100*c_d*C_r,'b','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+x2=x(2:end);c_f=100.*c_f(2:end);Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
+subplot(4,3,11);plot(x2*L, c_f*C_r,'g','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,[3 6]);plot(x2*L, Q_cross2*Q_r/x_r, 'b', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,[9 12]);plot(x2*L, Q_sin2*Q_r/x_r, 'g', 'LineWidth',lw);xlabel('x [m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 8
+if fig(8) == 1
+f=figure(8); % SWRO quantities
+f.Position = [100 1 1100 1000];
+x2=x(2:end);J_wf2=J_wf(2:end);J_sf2=J_sf(2:end);
+subplot(4,3,1);plot(x*swro_L, J_wd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,4);plot(x*swro_L, J_sd*J_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,d}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,2);plot(x2*swro_L, J_wf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{w,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,5);plot(x2*swro_L, J_sf2*J_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); legend('J_{s,f}^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;     
+subplot(4,3,7);plot(x*swro_L, P_d*p_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(4,3,8);plot(x*swro_L, P_f*p_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); legend('P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 5;
+subplot(4,3,10);plot(x*swro_L, 100*C_d*C_r,'Color', '#0072bD','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+C_f=100.*C_f(2:end);J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
+subplot(4,3,11);plot(x2*swro_L, C_f*C_r,'Color', '#77AC30','LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); legend('C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,[3 6]);plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', '#0072bD', 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+subplot(4,3,[9 12]);plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', '#77AC30', 'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay=gca; ay.YAxis.Exponent = 0;
+end
+%% figure 9
+if fig(9) == 1
+f=figure(3);  % ERD quantities
+f.Position = [100 1 900 500]; tiledlayout(1,4); 
+rot=0; rx=0; %rotation and shift of texts
+vals1 = [pE*p_r; pERD*p_r; 0; P_d(end)*p_r; p_d(end)*p_r];
+vals2 = 100*[cE; C_ERD; 0; C_d(end); c_d(1)];
+vals3 = [J_E*J_r; J_ERD*J_r; 0; J_d(end)*J_r; Q_d(1)*Q_r];
+%adjust for counter current / SWRO+ERD / 2ERDs
+if version(1)==1; vals2(5) = 100*c_d(end); vals3(5) =-Q_d(end)*Q_r;end
+if version(6)==1; vals1(5) = p_exit*p_r; vals2(5) = 100*c_exit; vals3(5) = Q_exit*J_r; end
+if version(6)==4; vals1(1)= pERD2*p_r; vals2(1)=100*C_ERD2; end
+nexttile
+b = bar(1,vals1, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[kg/ms^2]'}); ylim([0 1.1*max(vals1)]);ay=gca;ay.YAxis(1).Exponent=5;
+t=text(b(1).XEndPoints,b(1).YEndPoints,"P_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom');
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"P_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); 
+t=text(b(4).XEndPoints,b(4).YEndPoints,"P_{b}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); 
+t=text(b(5).XEndPoints,b(5).YEndPoints,"P_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); 
+nexttile
+b = bar(1,vals2, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[%]'}); ylim([0 1.12*max(vals2)]);ay=gca;ay.YAxis(1).Exponent=0;
+t= text(b(1).XEndPoints,b(1).YEndPoints,"C_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"C_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(4).XEndPoints,b(4).YEndPoints,"C_b^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(5).XEndPoints,b(5).YEndPoints,"C_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+nexttile
+b = bar(1,vals3, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[kg/ms]'}); ylim([0 1.1*max(vals3)]);ay=gca;ay.YAxis(1).Exponent=0;
+t=text(b(1).XEndPoints,b(1).YEndPoints,"J_{s}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"J_{s}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(4).XEndPoints ,b(4).YEndPoints,"J_{b}^{in}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(5).XEndPoints,b(5).YEndPoints,"J_{b}^{out}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot ;
+nexttile
+if version(6)>3
+lgd = legend([b(1) b(2) b(5) b(4)],'LP ERD 1 Seawater Inlet','HP ERD 1 Seawater Outlet','HP SWRO Brine Outlet','LP PRO Brine Inlet' , 'Location', 'EastOutside');
+else
+lgd = legend([b(1) b(2) b(5) b(4)],'Sea water inlet','Sea water outlet','Brine outlet','Brine inlet' , 'Location', 'EastOutside');
+end
+lgd.Layout.Tile = 4; axis off ; 
+end
+%% figure 10
+if fig(9) == 1 && version(6) ==4
+close all;f=figure(3);  % ERD quantities
+f.Position = [950 1 900 500]; tiledlayout(1,4); 
+rot=0; rx=0; %rotation and shift of texts
+vals1 = [pE*p_r; pERD2*p_r; 0; p_d(end)*p_r; p_exit*p_r];
+vals2 = 100*[cE; C_ERD2; 0; c_d(end); c_exit];
+vals3 = [J_E_sea*J_r; 2*J_E*J_r; 0; mixer_ERD*Q_d(end)*Q_r; Q_exit*Q_r];
+%adjust for counter current
+if version(1)==1; vals1(4) = p_d(1)*p_r; vals2(4)=c_d(1)*100;vals3(4)=-mixer_ERD*Q_d(1)*Q_r; end
+nexttile
+b = bar(1,vals1, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[kg/ms^2]'}); ylim([0 1.1*max(vals1)]);ay=gca;ay.YAxis(1).Exponent=5;
+t=text(b(1).XEndPoints,b(1).YEndPoints,"P_E",'HorizontalAlignment','center','VerticalAlignment','bottom');
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"P_{ERD 2}",'HorizontalAlignment','center','VerticalAlignment','bottom'); 
+if version(1) == 0; t=text(b(4).XEndPoints,b(4).YEndPoints,"p_d^1",'HorizontalAlignment','center','VerticalAlignment','bottom');  end
+if version(1) == 1; t=text(b(4).XEndPoints,b(4).YEndPoints,"p_d^0",'HorizontalAlignment','center','VerticalAlignment','bottom');  end
+t=text(b(5).XEndPoints,b(5).YEndPoints,"P_E",'HorizontalAlignment','center','VerticalAlignment','bottom'); 
+nexttile
+b = bar(1,vals2, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[%]'}); ylim([0 1.12*max(vals2)]);ay=gca;ay.YAxis(1).Exponent=0;
+t= text(b(1).XEndPoints,b(1).YEndPoints,"C_E",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints," C_{ERD 2}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+if version(1) == 0; t=text(b(4).XEndPoints,b(4).YEndPoints,"c_d^1",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; end
+if version(1) == 1; t=text(b(4).XEndPoints,b(4).YEndPoints,"c_d^0",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; end
+t=text(b(5).XEndPoints,b(5).YEndPoints,"C_{exit}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot; 
+nexttile
+b = bar(1,vals3, 'FaceColor', 'b'); b(4).FaceColor = "#77AC30"; b(5).FaceColor = [.2 .6 .5]; b(1).FaceColor = '#0072bD';
+xticklabels({'[kg/ms]'}); ylim([0 1.1*max(vals3)]);ay=gca;ay.YAxis(1).Exponent=0;
+t=text(b(1).XEndPoints,b(1).YEndPoints,"J_{sea}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+t=text(b(2).XEndPoints+rx ,b(2).YEndPoints,"   2J_{E}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+strg=[num2str(mixer_ERD), '*'];
+if version(1) == 0; t=text(b(4).XEndPoints ,b(4).YEndPoints,[strg,"$Q_d^1$"],'Interpreter','latex','HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;end
+if version(1) == 1; t=text(b(4).XEndPoints ,b(4).YEndPoints,[strg,"$-Q_d^0$"],'Interpreter','latex','HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;end
+t=text(b(5).XEndPoints,b(5).YEndPoints,"Q_{exit}",'HorizontalAlignment','center','VerticalAlignment','bottom'); t.Rotation = rot;
+nexttile
+lgd = legend([b(1) b(2) b(3) b(4)],'LP ERD 2 Seawater Inlet','HP ERD 2 Seawater Outlet','HP ERD 2 Brine Inlet','LP ERD 2 Brine Outlet' , 'Location', 'EastOutside');
+lgd.Layout.Tile = 4; axis off ; 
+end
+end
+%% figures for Paper
+if fig(11) == 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+f=figure(1); % SWRO quantities
+f.Position = [276.3333 196.3333 1100 1000];
+x2=x(2:end);
+% total mass flows
+subplot(3,2,1);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% concentrations
+subplot(3,2,2);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
+yyaxis right; C_f=C_f(2:end);
+plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Pressures
+subplot(3,2,3);lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Permeate flows 
+subplot(3,2,5); lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
+plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
+yyaxis right
+plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Osmotic/ Hydraulic pressure difference
+subplot(3,2,6);lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
+plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
+yyaxis right; 
+plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+f=figure(2); % PRO quantities
+f.Position = [936.3333 266.3333 1.0987e+03 632];
+x2=x(2:end);
+% total mass flows
+subplot(3,2,1);lc='#0072bD';rc='#77AC30';
+plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% concentrations
+subplot(3,2,2);lc='#0072bD';rc='#77AC30';
+plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
+yyaxis right; c_f=c_f(2:end);
+plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
+% Pressures
+subplot(3,2,3);lc='#0072bD';rc='#77AC30';
+plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Permeate flows 
+subplot(3,2,5); lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
+plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw);xlabel('x[m]','Fontsize',10);  ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
+yyaxis right
+plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+% Osmotic/ Hydraulic pressure difference
+subplot(3,2,6);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
+plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);xlabel('x[m]','Fontsize',10); ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
+yyaxis right; 
+plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
+end
+%% figures for Paper - all in individual plot
+if fig(11) == 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+f=figure(1);% SWRO quantities
+x2=x(2:end);
+% total mass flows
+lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, J_d*J_r,'Color', lc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0; xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*swro_L, J_f*J_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{RO}(x)','J_{f}^{RO}(x)','Location','northeast');xlim([0 swro_L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROa.png');
+% concentrations
+close all; f=figure(1);
+lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, 100*C_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;hold on
+yyaxis right; C_f=C_f(2:end);
+plot(x2*swro_L, 100*C_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{RO}(x)','C_f^{RO}(x)','Location','southeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROb.png');
+% Pressures
+close all; f=figure(6);
+lc='#0072bD';rc='#77AC30';
+plot(x*swro_L, P_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*swro_L, P_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{RO}(x)','P_f^{RO}(x)','Location','northeast');xlim([0 swro_L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROc.png');
+% Permeate flows 
+close all;f=figure(1);
+lc='k';rc='#b81414'; J_cross2=J_cross(2:end);J_sin2=J_sin(2:end);
+plot(x2*swro_L, J_cross2*J_r/swro_x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0;ylim([min( J_cross2*J_r/swro_x_r) 1.1*max( J_cross2*J_r/swro_x_r)]); hold on
+yyaxis right
+plot(x2*swro_L, J_sin2*J_r/swro_x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{RO}(x)','J_{s,in}^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([min( J_sin2*J_r/swro_x_r) 1.1*max( J_sin2*J_r/swro_x_r)]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROd.png');
+% Osmotic/ Hydraulic pressure difference
+close all;f=figure(1);
+lc='k';rc='#b81414'; osm_diff=swro_p_osm_d(2:end)-swro_p_osm_f(2:end); a=min([min(osm_diff*p_r) min((P_d-P_f)*p_r)]); b=max( [max(osm_diff*p_r), max((P_d-P_f)*p_r)]);
+plot(x*swro_L, (P_d-P_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ylim([a b]);ay=gca; ay.YAxis.Exponent = 6; hold on
+yyaxis right; 
+plot(x2*swro_L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{RO}(x)','\Delta \pi^{RO}(x)','Location','northeast');xlim([0 swro_L]);ylim([a b]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_SWROe.png');
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+f=figure(2); % PRO quantities
+x2=x(2:end);
+% total mass flows
+lc='#0072bD';rc='#77AC30';
+plot(x*L,  Q_d*Q_r,'Color', lc,'LineWidth',lw);ylabel('[kg/sm]','Fontsize',10); ay=gca; ay.YAxis.Exponent = 0;xlim([0 swro_L]); hold on
+yyaxis right
+plot(x*L, Q_f*Q_r,'Color', rc,'LineWidth',lw); ylabel('[kg/sm]','Fontsize',10); legend('J_{d}^{PRO}(x)','J_{f}^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROa.png');
+% concentrations
+close all; f=figure(2);lc='#0072bD';rc='#77AC30';
+plot(x*L, 100*c_d*C_r,'Color', lc,'LineWidth',lw); ylabel('[%]','Fontsize',10); ay=gca; ay.YAxis(1).Exponent = 0;hold on
+yyaxis right; c_f=c_f(2:end);
+plot(x2*L, 100*c_f*C_r,'Color', rc,'LineWidth',lw); ylabel('[%]','Fontsize',10); legend('C_d^{PRO}(x)','C_f^{PRO}(x)','Location','southeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;ay.YAxis(2).Exponent = 0;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROb.png');
+% Pressures
+close all; f=figure(2);lc='#0072bD';rc='#77AC30';
+plot(x*L, p_d*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 5; hold on
+yyaxis right
+plot(x*L, p_f*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('P_d^{PRO}(x)','P_f^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROc.png');
+% Permeate flows 
+close all; f=figure(2);lc='k';rc='#b81414'; Q_cross2=Q_cross(2:end);Q_sin2=Q_sin(2:end);
+plot(x2*L, Q_cross2*Q_r/x_r, 'Color', lc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10);ay=gca; ay.YAxis.Exponent = 0; hold on
+yyaxis right
+plot(x2*L, Q_sin2*Q_r/x_r, 'Color', rc, 'LineWidth',lw); ylabel('[kg/sm^2]','Fontsize',10); legend('J_{w,in}^{PRO}(x)','J_{s,in}^{PRO}(x)','Location','northeast');xlim([0 L]);ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc;
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROd.png');
+% Osmotic/ Hydraulic pressure difference
+close all; f=figure(2);lc='k';rc='#b81414';osm_diff=p_osm_d(2:end)-p_osm_f(2:end); %a=min([min(osm_diff*p_r) min((p_d-p_f)*p_r)]); b=max( [max(osm_diff*p_r), max((p_d-p_f)*p_r)]);
+plot(x*L, (p_d-p_f)*p_r,'Color', lc,'LineWidth',lw);ylabel('[Pa]','Fontsize',10); ay=gca;ay.YAxis.Exponent = 6; hold on; %ylim([a b]); 
+yyaxis right; 
+plot(x2*L, osm_diff*p_r,'Color', rc,'LineWidth',lw); ylabel('[Pa]','Fontsize',10); legend('\Delta P^{PRO}(x)','\Delta \pi^{PRO}(x)','Location','northeast');xlim([0 L]); ay.YAxis(1).Color = lc; ay.YAxis(2).Color = rc; %ylim([a b]);
+set(gcf,'color','w'); f = gcf; exportgraphics(f,'Simulation_PROe.png');
+end
+
+
+%ax = axes; yyaxis('left'); yyaxis('right');
+%ax.YAxis(1).Color = [0 0 0]; ax.YAxis(2).Color = [0 0 0];
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+end
+end
\ No newline at end of file
diff --git a/fun_1.m b/fun_1.m
index 406a9d28f58bf1a9b8ca095d5c53ec82da1e2141..1cfef86ce7f11b65b7582dca2d9b404a052a2ec3 100644
--- a/fun_1.m
+++ b/fun_1.m
@@ -28,8 +28,8 @@ function [output1, output2, output3] = fun_1(input,option_data,obj,option_mesh,o
 %
 tic;
 %% Read data
-if option_data == 0; DATA = @(x)Sim_0_data(input); end
-
+if option_data == 0; DATA = @(x)Sim_0a_data(input); end
+if option_data == 0.1; DATA = @(x)Sim_0b_data(input); end
 
 
 
@@ -80,108 +80,96 @@ if version(1)==1; y_init = [cE; J_wd_0; J_sf_0; J_wf_0; Pd_0; Pf_0; cE; -Q_wd_0;
 solinit = bvpinit(x,y_init);
 
 %% Solve the BVP
-if version(6) == 0; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary2(ya,yb, DATA),solinit,ode_options); end
-if version(6) == 1; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary1(ya,yb, DATA),solinit,ode_options); end
-if version(6) == 2; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary2(ya,yb, DATA),solinit,ode_options); end
-if version(6) == 3; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary3(ya,yb, DATA),solinit,ode_options); end       
-if version(6) == 4; sol = bvp5c(@(x,J_p)ERD_model(x, J_p, DATA), @(ya,yb)ERD_boundary4(ya,yb, DATA),solinit,ode_options); end   
+if version(6) == 0; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 1; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary1(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 2; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary2(ya,yb, DATA),solinit,ode_options); end
+if version(6) == 3; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary3(ya,yb, DATA),solinit,ode_options); end       
+if version(6) == 4; sol = bvp5c(@(x,J_p)ODEsystem(x, J_p, DATA), @(ya,yb)Boundary4(ya,yb, DATA),solinit,ode_options); end   
             
 %% Evaluate the solution
 n=100;
 x=linspace(0,1,n);
 y = deval(sol,x); Y = y'; Y = real(Y);
-C_d=zeros(1,n);C_f=zeros(1,n);J_sd=zeros(1,n);J_wd=zeros(1,n);J_d=zeros(1,n);J_sf=zeros(1,n);J_wf=zeros(1,n);J_f=zeros(1,n);P_d=zeros(1,n);P_f=zeros(1,n);swro_p_osm_d=zeros(1,n);swro_p_osm_f=zeros(1,n);swro_beta=zeros(1,n);Y(:,4+6)=abs(Y(:,4+6));swro_del_c=zeros(1,n);J_cross=zeros(1,n);J_sin=zeros(1,n);
-c_d=zeros(1,n);c_f=zeros(1,n);Q_sd=zeros(1,n);Q_wd=zeros(1,n);Q_d=zeros(1,n);Q_sf=zeros(1,n);Q_wf=zeros(1,n);Q_f=zeros(1,n);p_d=zeros(1,n);p_f=zeros(1,n);p_osm_d=zeros(1,n);p_osm_f=zeros(1,n);beta=zeros(1,n);del_c=zeros(1,n);Q_cross=zeros(1,n);Q_sin=zeros(1,n); local_ro_d=zeros(1,n);local_ro_f=zeros(1,n);swro_local_ro_d=zeros(1,n);swro_local_ro_f=zeros(1,n);
-
-    for i=1:n
-    %% SWRO evaluation
-    % Concentrations
-    C_d(i)= Y(i,1);                 % concentration of draw side
-    C_f(i)= Y(i,3)./Y(i,4);         % concentration of fresh side
-    % Flow rates in draw side
-    J_sd(i)= Y(i,1)*Y(i,2);         % salt flow rate in draw side
-    J_wd(i)= Y(i,2);                % water flow rate in draw side
-    J_d(i) = Y(i,1)*Y(i,2)+Y(i,2);  % total flow rate draw side
-    % Flow rates in fresh side
-    J_sf(i)= Y(i,3);                % salt flow rate in fresh side
-    J_wf(i)= Y(i,4);                % water flow rate in fresh side
-    J_f(i) = Y(i,3)+Y(i,4);         % total flow rate fresh side
-    % Pressures 
-    P_d(i)=Y(i,5);                  % Pressure draw side
-    P_f(i)=Y(i,6);                  % Pressure fresh side
-
-    % Osmotic pressure
-    swro_p_osm_d(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,1)/Ms)/p_r;
-    swro_p_osm_f(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,3)/Ms./Y(i,4))/p_r;
-    % Concentration difference
-    swro_del_c(i) = Y(i,1)./(Y(i,1) + 1) - Y(i,3)./(Y(i,3) + Y(i,4));
-    % Local density
-    swro_local_ro_d(i) = (Y(i,1) + 1)./(ro_water*Y(i,1)/ro_salt + 1);           % local density of draw side
-    swro_local_ro_f(i) = (Y(i,3) + Y(i,4))./(ro_water*Y(i,3)/ro_salt + Y(i,4)); % local density of fresh side
-    % B(x)
-    swro_beta(i) = (1 - swro_R)*((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)))/swro_R;
-    if version(2) == 0 ;swro_beta(i)= swro_beta_fix; end
-    if version(4) == 0 ;swro_beta(i)= 0; end
-    % J_w,in and J_s,in
-    
-    J_cross(i) = ((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)));
-    if version(4) == 0; J_cross(i) = ((Y(i,5) - Y(i,6)) - sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)));  end 
-    J_sin(i)=swro_beta(i)*swro_del_c(i);
- 
-    % J_w,in and J_s,in with ECP and ICP
-    if version(7)==1
-    J_cross(i) = ((Y(i,5)-Y(i,6))-sigma.*(swro_p_osm_d(i) - swro_p_osm_f(i)) + swro_beta(i)*(Y(i,5)-Y(i,6))*(swro_KK +1/swro_KD + 1/swro_KF)) ./   (1 + swro_alpha*swro_beta(i)*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d(i)/swro_KD -swro_p_osm_f(i)*(swro_KK + 1/swro_KF) ));
-    J_sin(i) = swro_beta(i)*( swro_del_c(i)-J_cross(i)*( Y(i,1)./(Y(i,1) + 1)./swro_KD +  Y(i,3)./(Y(i,3) + Y(i,4))*(swro_KK+1/swro_KF)))/(1+swro_alpha*swro_beta(i)*(swro_KK+1/swro_KF+1/swro_KD));
-    end
-    
+
+%% SWRO evaluation
+C_d = Y(:,1);
+C_f = Y(:,3)./Y(:,4);
+J_sd = Y(:,1).*Y(:,2);
+J_wd = Y(:,2); 
+J_d = Y(:,1).*Y(:,2)+Y(:,2);
+J_sf = Y(:,3);
+J_wf = Y(:,4);
+J_f = Y(:,3)+Y(:,4);
+P_d = Y(:,5); 
+P_f = Y(:,6); 
+% same quantities as in "ODEsystem.m" 
+swro_p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1)/Ms)/p_r;
+swro_p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3)/Ms./Y(:,4))/p_r;
+swro_del_c = Y(:,1)./(Y(:,1) + 1) - Y(:,3)./(Y(:,3) + Y(:,4));
+swro_local_ro_d = (Y(:,1) + 1)./(ro_water*Y(:,1)/ro_salt + 1);
+swro_local_ro_f = (Y(:,3) + Y(:,4))./(ro_water*Y(:,3)/ro_salt + Y(:,4));
+swro_beta = (1 - swro_R)*((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))/swro_R;
+if version(2) == 0 ;swro_beta= swro_beta_fix.*ones(1,n); end
+if version(4) == 0 ;swro_beta=zeros(1,n); end
+% J_w,in and J_s,in
+J_cross = ((Y(:,5) - Y(:,6)) - sigma.*(swro_p_osm_d - swro_p_osm_f))./(1 + p_r*swro_alpha*swro_KK*sigma.*(swro_p_osm_d - swro_p_osm_f));
+if version(4) == 0; J_cross = ((Y(:,5) - Y(:,6)) - sigma*(swro_p_osm_d - swro_p_osm_f));  end 
+J_sin=swro_beta.*swro_del_c;
+if version(7)==1
+J_cross = ((Y(:,5)-Y(:,6))-sigma.*(swro_p_osm_d - swro_p_osm_f) + swro_beta.*(Y(:,5)-Y(:,6)).*(swro_KK +1/swro_KD + 1/swro_KF)) ./ (1 + swro_alpha.*swro_beta.*(1/swro_KD + 1/swro_KF + swro_KK) + p_r*swro_alpha*sigma.*(-swro_p_osm_d./swro_KD -swro_p_osm_f.*(swro_KK + 1/swro_KF) ));
+J_sin = swro_beta.*( swro_del_c-J_cross.*( Y(:,1)./(Y(:,1) + 1)./swro_KD +  Y(:,3)./(Y(:,3) + Y(:,4)).*(swro_KK+1/swro_KF)))./(1+swro_alpha.*swro_beta.*(swro_KK+1/swro_KF+1/swro_KD));
+end
+
 %% PRO evaluation
-    % Concentrations
-    c_d(i)= Y(i,1+6);                    % concentration of draw side
-    c_f(i)= Y(i,3+6)./Y(i,4+6);          % concentration of fresh side
-    % Flow rates in draw side
-    Q_sd(i)= Y(i,1+6)*Y(i,2+6);          % salt flow rate in draw side
-    Q_wd(i)= Y(i,2+6);                   % water flow rate in draw side
-    Q_d(i) = Y(i,1+6)*Y(i,2+6)+Y(i,2+6); % total flow rate draw side
-    % Flow rates in fresh side 
-    Q_sf(i)= Y(i,3+6);                   % salt flow rate in fresh side
-    Q_wf(i)= Y(i,4+6);                   % water flow rate in fresh side
-    Q_f(i) = Y(i,3)+Y(i,4+6);            % total flow rate fresh side
-    % Pressures
-    p_d(i)=Y(i,5+6);                     % Pressure draw side
-    p_f(i)=Y(i,6+6);                     % Pressure fresh side
-
-    % Osmotic pressure
-    p_osm_d(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,1+6)/Ms)/p_r;
-    p_osm_f(i) = ro_water*Rw*T0*log(1 + 2*Mw*Y(i,3+6)/Ms./Y(i,4+6))/p_r;
-    % Concentration difference
-    del_c(i) = Y(i,1+6)./(Y(i,1+6) + 1) - Y(i,3+6)./(Y(i,3+6) + Y(i,4+6));
-    % local density
-    local_ro_d(i) = (Y(i,1+6) + 1)./(ro_water*Y(i,1+6)/ro_salt + 1);            	
-    local_ro_f(i) = (Y(i,3+6) + Y(i,4+6))./(ro_water*Y(i,3+6)/ro_salt + Y(i,4+6)); 
-    % B(x) 
-    beta(i) = (1 - R)*((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6)))/R;
-    if version(3) == 0 ; beta(i) =  beta_fix; end
-    if version(5) == 0 ; beta(i) = 0; end
-    % Q_w,in and Q_s,in
-    Q_cross(i) = (p_osm_d(i) - p_osm_f(i) - (Y(i,5+6) - Y(i,6+6)) - (Y(i,5+6) - Y(i,6+6))*beta(i)*KK*alpha*p_r)./(1 + KK*alpha*p_r*(beta(i) +p_osm_f(i)));
-    if version(5) == 0; Q_cross(i) = ((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6))) ; end
-    Q_sin(i)=beta(i)*del_c(i);
-    
-    if version(8)==1 %% version with ICP and ECP
-    Q_cross(i)=((p_osm_d(i) - p_osm_f(i)) - (Y(i,5+6) - Y(i,6+6)) - (Y(i,5+6) - Y(i,6+6))*beta(i)*(KK+1/KF+1/KD)*p_r*alpha)./(1 + alpha*beta(i)*(KK+1/KF+1/KD) + alpha*p_r*(p_osm_f(i)*(KK+1/KF) + p_osm_d(i)/KD));
-    Q_sin(i)=beta(i)*( del_c(i)-Q_cross(i)*( Y(i,1+6)./(Y(i,1+6) + 1)./KD +  Y(i,3+6)./(Y(i,3+6) + Y(i,4+6))*(KK+1/KF)))/(1+alpha*beta(i)*(KK+1/KF+1/KD));
-    end
-    
-    end
-    % In ideal cases take mean of some quantities (one could instead use ylim on y-axis) 
-    if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end
-    if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end
+c_d = Y(:,1+6);                    % concentration of draw side
+c_f = Y(:,3+6)./Y(:,4+6);          % concentration of fresh side
+Q_sd = Y(:,1+6).*Y(:,2+6);          % salt flow rate in draw side
+Q_wd = Y(:,2+6);                   % water flow rate in draw side
+Q_d = Y(:,1+6).*Y(:,2+6)+Y(:,2+6); % total flow rate draw side
+Q_sf = Y(:,3+6);                   % salt flow rate in fresh side
+Q_wf = Y(:,4+6);                   % water flow rate in fresh side
+Q_f = Y(:,3)+Y(:,4+6);            % total flow rate fresh side
+p_d = Y(:,5+6);                     % Pressure draw side
+p_f = Y(:,6+6);                     % Pressure fresh side
+% same quantities as in "ODEsystem.m" 
+p_osm_d = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,1+6)/Ms)/p_r;
+p_osm_f = ro_water*Rw*T0*log(1 + 2*Mw*Y(:,3+6)/Ms./Y(:,4+6))/p_r;
+del_c = Y(:,1+6)./(Y(:,1+6) + 1) - Y(:,3+6)./(Y(:,3+6) + Y(:,4+6));
+local_ro_d = (Y(:,1+6) + 1)./(ro_water*Y(:,1+6)/ro_salt + 1);            	
+local_ro_f = (Y(:,3+6) + Y(:,4+6))./(ro_water*Y(:,3+6)/ro_salt + Y(:,4+6)); 
+beta = (1 - R)*((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)))./R;
+if version(3) == 0 ; beta =  beta_fix.*ones(1,n); end
+if version(5) == 0 ; beta = zeros(1,n); end
+% Q_w,in and Q_s,in
+Q_cross = (p_osm_d - p_osm_f - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*KK.*alpha.*p_r)./((ones(1,n) + KK*alpha*p_r.*(beta +p_osm_f'))');
+if version(5) == 0; Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6))) ; end
+Q_sin = beta'.*del_c;
+if version(8)==1 
+Q_cross = ((p_osm_d - p_osm_f) - (Y(:,5+6) - Y(:,6+6)) - (Y(:,5+6) - Y(:,6+6)).*beta'.*(KK+1./KF+1./KD).*p_r.*alpha)./(ones(1,n) + alpha.*beta.*(KK+1./KF+1./KD) + alpha.*p_r.*(p_osm_f'.*(KK+1./KF) + p_osm_d'./KD))';
+Q_sin = beta'.*( del_c-Q_cross.*( Y(:,1+6)./(Y(:,1+6) + 1)./KD +  Y(:,3+6)./(Y(:,3+6) + Y(:,4+6)).*(KK+1/KF)))./(1+alpha.*beta');
+end
     
-%% Analysis of Recovery rates
+%% Freswater production and recovery rates
 Freshwater = J_wf(end)*J_r; %in [kg/s] 
 FW = (Freshwater*swro_Z/(swro_local_ro_f(end)*rho_r))*3600; %in [m^3/h] 
 SWRO_Recovery =(J_f(end)./J_d(1))*100;     % SWRO freshwater recovery [%]
 PRO_Recovery = (1- Q_f(end)./Q_f(1))*100;  % PRO recovery rate [%]
+% warning flags:
+    if FW < 0 
+        fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n');
+    end
+    if SWRO_Recovery > 100 
+        fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n');
+    end
+    if SWRO_Recovery < 0 
+        fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
+    end
+    if PRO_Recovery > 100 
+        fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
+    end
+    if PRO_Recovery < 0 
+        fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+    end
 
 %% Output of System
 J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; J_wE = J_E/(cE+1);
@@ -202,14 +190,13 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
         rho_exit= (ro_salt*(c_exit+1))/(c_exit*ro_water+ro_salt); p_exit=pE;
         f_1= (Q_exit*swro_Z/rho_exit*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2 );
         f_2= (J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2); 
-        f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+        f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
         pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-        if pERD > pMax; pERD=pMax; end
-        if pERD < 0; pERD=pMax; end  
-        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
         W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0;
+            if W_p3 < 0 
+                fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
+            end
+
     end
     %% Version(6)=2 --> only PRO
     if version(6)==2
@@ -220,33 +207,25 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
         if version(1) ==1; W_p2= 1/LP_eff * (p_d(end)-pE)*(abs(Q_d(end))*Z)./local_ro_d(end); end 
         W_p2 = W_p2*W_r;
     end
-    
+
     %% Version(6)=3 --> SWRO-PRO hybrid system (with one ERD)
     if version(6)==3
         rho_ERD = V_m*(swro_local_ro_d(end) - rho_E)+rho_E;
-        C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);
-        J_wERD = J_E/(C_ERD+1); qq=Q_r/J_r;
+        C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);  
+        qq=Q_r/J_r;
             if version(1) == 0 %co-current
                 f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2;  
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-                f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
-                if pERD > pMax; pERD=pMax; end
-                %if pERD < 0; pERD=pMax; end   
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*(Q_d(end)*Z)/local_ro_d(end); 
             end
             if version (1) == 1 %counter-current   
                 f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_E/A_ERD)^2  ;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;     
-                f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;       
-                %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
-                %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
-                if pERD > pMax; pERD=pMax; end
-                %if pERD < 0; pERD=pMax; end     
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs(Q_d(1))*Z)/local_ro_d(1);
             end
@@ -260,8 +239,7 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
         C_ERD2= -ro_salt*(rho_ERD2-1)/(rho_ERD2*ro_water-ro_salt);
         rho_ERD= V_m*(swro_local_ro_d(end) - rho_ERD2)+rho_ERD2;
         C_ERD = -ro_salt*(rho_ERD-1)/(rho_ERD*ro_water-ro_salt);
-        J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; 
-        J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1);
+        J_E = (C_d(1)*J_wd(1)+J_wd(1))/2; J_ERD = J_E; %J_wE = J_E/(C_ERD2+1); J_wERD= J_ERD/(C_ERD+1);
         qq=Q_r/J_r;
         J_E_sea = 2*J_E;
         J_wE_sea = J_E_sea/(cE+1);
@@ -278,45 +256,33 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
             if version(1) == 0 %co-current
                 f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end)/A_ERD)^2;  
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
-                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
-                if pERD2 > pMax; pERD2=pMax; end
-                %if pERD2 < 0; pERD2=pMax; end  
-                
             end
             if version (1) == 1 %counter-current 
                 f_1= qq*Q_exit*Z/rho_exit*mix_density*(J_E_sea*swro_Z/rho_E/A_ERD)^2;  
                 f_2= J_E*swro_Z/rho_ERD2*mix_density*(qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1)/A_ERD)^2;  
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD2 = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-                f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
-                if pERD2 > pMax; pERD2=pMax; end
-                %if pERD2 < 0; pERD2=pMax; end  
+
             end     
             %calculation of Pressure P_ERD
             if version(1) == 0 %co-current
                 f_1= qq*Q_d(1)*Z/local_ro_d(1)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2;  
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2;  
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
                 pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
-                f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
-                if pERD > pMax; pERD=pMax; end
-                %if pERD < 0; pERD=pMax; end  
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(end)-pE)*((1-mixer_ERD)*Q_d(end)*Z)/local_ro_d(end); 
             end
             if version (1) == 1 %counter-current   
                 f_1= abs(qq*Q_d(end))*Z/local_ro_d(end)*mix_density*(J_E*swro_Z/rho_ERD2/A_ERD)^2  ;
                 f_2= J_ERD*swro_Z/rho_ERD*mix_density*(J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end)/A_ERD)^2  ;
-                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1);
-                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
-                f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
-                if pERD > pMax; pERD=pMax; end 
-                %if pERD < 0; pERD=pMax; end  
+                f_ERD= ERD_fric*(ERD_eff*f_2 - f_1)/100;
+                pERD = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;           
                 W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD;
                 W_t = T_eff * (p_d(1)-pE)*(abs((1-mixer_ERD)*Q_d(1))*Z)/local_ro_d(1);
-            end         
+            end
     % Recalculation for Pump 1
     W_p1 = 1/HP_eff * (P_d(1)-pERD2)*(J_E *swro_Z)/rho_ERD2; W_p1 = W_p1*swro_W_r;
     end
@@ -333,6 +299,13 @@ if version(1)==1; SE_net = W_net./(Q_f(1)*Q_r*Z/(rho_r*local_ro_f(1)))/1000/3600
 SEC_net = W_net*(swro_local_ro_f(end)*rho_r)./(J_f(end)*J_r*swro_Z)/1000/3600; %in [kWh/m^3]  
 Rev= pw*FW + pe*SEC_net*FW; %in [$/h]
 tt=toc;
+
+
+
+
+
+
+
 %% if BVP-solver failed
 catch
     sol.stats.maxerr =1000;
@@ -808,5 +781,77 @@ end
 
 
 
+
+end
 end
-end
\ No newline at end of file
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%% Removed comments
+
+%Y(:,4+6)=abs(Y(:,4+6)) % ??why is this here??
+
+% -3.9191    0.3247    0.0131   32.8727    7.1201       NaN       NaN
+%
+% 1.0e+03 *
+% -1.2742   -0.9744         0   -0.4873   -0.0010    0.1885
+
+% In ideal cases take mean of some quantities. 
+% In ideal cases salt will fluctuate around 1e-16. We take mean here for 
+% 
+% This is just for figures. (), otherwise one could use ylim ) 
+%if version(4)==0 ; J_f=J_wf; J_sf=zeros(1,n); J_sd=mean(J_sd)*ones(1,n); C_f=zeros(1,n); end
+%if version(5)==0 ; Q_f=Q_wf; Q_sf=zeros(1,n); Q_sd=mean(Q_sd)*ones(1,n); c_f=zeros(1,n); end
+
+        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_exit*Q_exit*swro_Z/rho_exit ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
+        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
+        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
+
+        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
+
+        %f_ERD=0;pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pE*J_E*swro_Z/rho_E -f_ERD )/J_ERD/swro_Z;       
+        %display(['pERD:              ',num2str(pERD)]) % pERD use in the model
+        %display(['pERD_theoretical:  ',num2str(pMax)]) % theoretical pERD with no leak, friction perm and 100% efficiency
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end     
+
+        %J_wERD = J_E/(C_ERD+1);
+
+        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(end)*qq*mixer_ERD*Q_d(end)*Z*(1-eta_ERD)/local_ro_d(end) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;
+        %if pERD2 > pMax; pERD2=pMax; end
+        %if pERD2 < 0; pERD2=pMax; end  
+
+        %f_ERD=0;pMax = rho_ERD2*(ERD_eff*( p_d(1)*qq*mixer_ERD*abs(Q_d(1))*Z*(1-eta_ERD)/local_ro_d(1) - p_exit*qq*Q_exit*Z/rho_exit ) +  pE*J_E_sea*swro_Z/rho_E -f_ERD )/2/J_E/swro_Z;            
+        %if pERD2 > pMax; pERD2=pMax; end
+        %if pERD2 < 0; pERD2=pMax; end  
+
+        %f_ERD=0; pMax = rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(1)*qq*Q_d(1)*Z/local_ro_d(1) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;
+        %if pERD > pMax; pERD=pMax; end
+        %if pERD < 0; pERD=pMax; end  
+
+        %f_ERD=0; pMax=rho_ERD*(ERD_eff*( P_d(end)*J_d(end)*swro_Z*(1-eta_ERD)/swro_local_ro_d(end) - p_d(end)*abs(qq*Q_d(end))*Z/local_ro_d(end) ) +  pERD2*J_E*swro_Z/rho_ERD2 -f_ERD )/J_ERD/swro_Z;            
+        %if pERD > pMax; pERD=pMax; end 
+        %if pERD < 0; pERD=pMax; end  
\ No newline at end of file