diff --git a/.gitignore b/.gitignore
index 7aa5fc0a2f338b13b8bd7f9765e9a051d98a96a0..1d2761d31d0db25e70455691eb1a8e97e306c075 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
+BayesValidRox/.spyderproject/
 tests/BeamTest/Polynomials/
 tests/BeamTest/surrogateBeam/*.zip
 tests/BraninTest/Polynomials/*.npy
@@ -56,6 +57,7 @@ tests/MECBM/SurrogateModel/SeqDesign_25_100_DPP/
 tests/MECBM/SurrogateModel/SeqDesign_25_150_DKL/
 tests/MECBM/SurrogateModel/__pycache__/read_MECBM.cpython-37.pyc
 tests/MECBM/model_exe_orig.sh
+tests/Fetzer/Outputs_*
 BayesValidRox/BayesInference/__pycache__/Discrepancy.cpython-36.pyc
 BayesValidRox/surrogate_models/__pycache__/Exploration.cpython-36.pyc
 BayesValidRox/tests/AnalyticalFunction/AnalyticFunc_Results.pkl
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index ff1c83604bba1e847b62138eff8787d1815fdb6a..af4a7443bf8291039d3f4756f86441a850bec60f 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -590,6 +590,7 @@ class PostProcessing:
         
         """
         PCEModel = self.PCEModel
+        totalNSamples = PCEModel.ExpDesign.X.shape[0]
         
         plotList = ['Modified LOO error', 'Validation error', 'KLD', 'BME', 'RMSEMean', 'RMSEStd']
         seqList = [PCEModel.SeqModifiedLOO,PCEModel.seqValidError, PCEModel.SeqKLD, PCEModel.SeqBME, 
@@ -647,7 +648,7 @@ class PostProcessing:
                         patch.set(facecolor=fill_color[idx])
 
                     
-                step = 5
+                step = PCEModel.ExpDesign.NrofNewSample if PCEModel.ExpDesign.NrofNewSample!=1 else 1
                 edge_color = ['red', 'blue', 'green']
                 fill_color = ['tan', 'cyan', 'lightgreen']
                 fig, ax = plt.subplots()
@@ -672,7 +673,7 @@ class PostProcessing:
                             plt.axhline(y=0.0, xmin=0, xmax=1, c='green', ls='--', lw=2)
                     
                     # Plot each UtilFuncs
-                    labels = np.arange(NrofInitSamples, allErrors.shape[1]+NrofInitSamples, step) 
+                    labels = np.arange(NrofInitSamples, totalNSamples, step)
                     draw_plot(allErrors[:,::step], labels, edge_color, fill_color, idx)
                     
                  
@@ -708,8 +709,8 @@ class PostProcessing:
                 
                 for idx, name in enumerate(nameUtil):
                     SeqValues = SeqDict[name]
-                    x_idx = np.arange(NrofInitSamples, SeqValues.shape[0]+NrofInitSamples)
-                    
+                    step = PCEModel.ExpDesign.NrofNewSample if PCEModel.ExpDesign.NrofNewSample!=1 else 1
+                    x_idx = np.arange(NrofInitSamples, totalNSamples+1, step)
                     if plot == 'KLD' or plot == 'BME':
                         # BME convergence if refBME is provided
                         if refBME_KLD is not None:
@@ -736,7 +737,7 @@ class PostProcessing:
                         plt.semilogy(x_idx, SeqValues.mean(axis=1), marker=markers[idx],
                                      ls='--', lw=2, color=colors[idx], label=name.split("_rep",1)[0])
                 
-                plt.xticks(x_idx[::5])
+                plt.xticks(x_idx)
                 plt.xlabel('Number of runs')
                 plt.ylabel(plotLabel)
                 plt.legend()
diff --git a/BayesValidRox/surrogate_models/ExpDesigns.py b/BayesValidRox/surrogate_models/ExpDesigns.py
index 72e17c1bfdb36a76dea56da633796ed2d2aa2bb8..5f8904cd1d130c7d69f72d77b108aa2f13f8e6a2 100644
--- a/BayesValidRox/surrogate_models/ExpDesigns.py
+++ b/BayesValidRox/surrogate_models/ExpDesigns.py
@@ -29,6 +29,7 @@ class ExpDesigns:
         self.NCandidate = 1
         self.PostSnapshot = False
         self.stepSnapshot = 2
+        self.UtilityFunction = 'None'
         self.nReprications = 1
         self.parNames = []
         self.MAP = []
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 3d5ff6c4abf8a6b9f24ea614b0bd5d140a6c7919..c0cad88de871f3e89255f2dd62ac58bc44232928 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -995,19 +995,21 @@ class aPCE:
                 temp = totalScore.copy()
                 temp[np.isnan(totalScore)] = -np.inf
                 sorted_idxtotalScore = np.argsort(temp)[::-1]
-                bestIdx=sorted_idxtotalScore[:NrofNewSample][0]
+                bestIdx = sorted_idxtotalScore[:NrofNewSample]
                 
                 # select the requested number of samples
                 if ExploreMethod == 'Voronoi':
-                    X_can = exploration.closestPoints[bestIdx]
-                    
-                    #plotter(self.ExpDesign.X, X_can, ExploreMethod, scoreExploration=None)
-                    
-                    # Calculate the maxmin score for the region of interest
-                    newSamples, maxminScore = exploration.getMCSamples(allCandidates=X_can)
-                    
-                    # select the requested number of samples
-                    Xnew = np.array([newSamples[np.argmax(maxminScore)]])
+                    Xnew = np.zeros((NrofNewSample,ndim))
+                    for i, idx in enumerate(bestIdx):
+                        X_can = exploration.closestPoints[idx]
+                        
+                        #plotter(self.ExpDesign.X, X_can, ExploreMethod, scoreExploration=None)
+                        
+                        # Calculate the maxmin score for the region of interest
+                        newSamples, maxminScore = exploration.getMCSamples(allCandidates=X_can)
+                        
+                        # select the requested number of samples
+                        Xnew[i] = newSamples[np.argmax(maxminScore)]
                 else:
                     Xnew = allCandidates[sorted_idxtotalScore[:NrofNewSample]]
                 
@@ -1071,19 +1073,23 @@ class aPCE:
                 
                 #print("ExploitScore:\n", ExploitScore)
                 
-                # Find the index for max in ExploreScore
-                max_idx = np.where(ExploitScore == np.nanmax(ExploitScore))[0][0]
-                
-                # select the samples in the Vornoi cell with the highest ExploreScore
-                selectedCandidates = exploration.closestPoints[max_idx]
-                #print("selectedCandidates:\n", selectedCandidates)
-                
-                # Calculate the maxmin score for the region of interest
-                newSamples, maxminScore = exploration.getMCSamples(allCandidates=selectedCandidates)
-
-                # select the requested number of samples
-                Xnew = np.array([newSamples[np.argmax(maxminScore)]])
+                # find an optimal point subset to add to the initial design
+                # by maximization of the utility score and taking care of NaN values
+                temp = ExploitScore.copy()
+                bestIdx = np.argsort(temp)[::-1][:NrofNewSample]
+                #bestIdx = sorted_idxtotalScore[:NrofNewSample]
                 
+                Xnew = np.zeros((NrofNewSample,ndim))
+                for i, idx in enumerate(bestIdx):
+                    X_can = exploration.closestPoints[idx]
+                    
+                    #plotter(self.ExpDesign.X, X_can, ExploreMethod, scoreExploration=None)
+                    
+                    # Calculate the maxmin score for the region of interest
+                    newSamples, maxminScore = exploration.getMCSamples(allCandidates=X_can)
+                    
+                    # select the requested number of samples
+                    Xnew[i] = newSamples[np.argmax(maxminScore)]
             
             elif ExploitMethod == 'alphabetic':
                 # ------- EXPLOITATION: ALPHABETIC -------
@@ -1778,7 +1784,7 @@ class aPCE:
                         for OutName in OutputName:
                             Scores_all.append(list(validError[OutName].values()))
                         ValidError = [item for sublist in Scores_all for item in sublist]
-                        print("\nValidError:", ValidError)
+                        print("\nUpdated ValidError:", ValidError)
                         print("\n")
                         
                     # Extract Modified LOO from Output
@@ -1798,7 +1804,7 @@ class aPCE:
                     # if it increases consistently stop the iterations
                     n_checks = 3
                     itrNr = Xfull.shape[0] - Xinit.shape[0]
-                    if itrNr > n_checks:
+                    if itrNr > n_checks * self.ExpDesign.NrofNewSample:
                         ss = np.sign(SeqModifiedLOO - ModifiedLOO) #ss<0 error increasing
                         errorIncreases = np.sum(np.mean(ss[-2:], axis=1)) <= -1*n_checks
                     
@@ -1812,7 +1818,7 @@ class aPCE:
                         break
                     
                     else:
-                         prevPCEModel = PCEModel #PCEModel
+                         prevPCEModel = PCEModel
 
                     
                     # Store updated ModifiedLOO 
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
index 30f156b57d355f521fdc9c865046690b195cd024..b288801ec56976e83abead7731754500fd4531e3 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
@@ -107,12 +107,13 @@ if __name__ == "__main__":
     # hypercube sampling of the input model or user-defined values of X and/or Y:
     MetaModelOpts.addExpDesign()
     
-    NrofInitSamples = 50
+    NrofInitSamples = 10
     MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
     MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)MC 2)LHS 3)PCM 4)LSCM 5)user
     MetaModelOpts.ExpDesign.Method = 'sequential'  # 1) normal  2) sequential
     
     # Sequential experimental design (needed only for sequential ExpDesign)
+    MetaModelOpts.ExpDesign.NrofNewSample = 5
     MetaModelOpts.ExpDesign.MaxNSamples = 75
     MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-6
     
diff --git a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
index f0b1a250f35bd1515e196242343fe8dfcc7a240c..e30bf4ec24703d68ece3f5f21696bafb1b62f556 100644
--- a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
+++ b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
@@ -115,6 +115,7 @@ if __name__ == "__main__":
     #MetaModelOpts.ExpDesign.Y = ExpDesign_Y
     
     # Sequential experimental design (needed only for sequential ExpDesign)
+    MetaModelOpts.ExpDesign.NrofNewSample = 10
     MetaModelOpts.ExpDesign.MaxNSamples = 500
     MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-4