diff --git a/examples/Data/DoubleLayerCapacity/Compressibility.npz b/examples/Data/DoubleLayerCapacity/Compressibility.npz
new file mode 100644
index 0000000000000000000000000000000000000000..c55f3df75b022d4ec4f5bbdbc61455cf642727b2
Binary files /dev/null and b/examples/Data/DoubleLayerCapacity/Compressibility.npz differ
diff --git a/examples/Data/DoubleLayerCapacity/Lambda.npz b/examples/Data/DoubleLayerCapacity/Lambda.npz
index 3d1c273d650aff43776c4de88ceb8e7b09d2aa3d..519216e909dc9a2de0775695b9bab97e598e29a5 100644
Binary files a/examples/Data/DoubleLayerCapacity/Lambda.npz and b/examples/Data/DoubleLayerCapacity/Lambda.npz differ
diff --git a/examples/Data/DoubleLayerCapacity/Molarity.npz b/examples/Data/DoubleLayerCapacity/Molarity.npz
index a44c6681d29e00040ba493f9c81ca6567bce63e5..6a5f7c009272eab67b3a30201b6de7b919076ce9 100644
Binary files a/examples/Data/DoubleLayerCapacity/Molarity.npz and b/examples/Data/DoubleLayerCapacity/Molarity.npz differ
diff --git a/examples/Data/DoubleLayerCapacity/Solvation.npz b/examples/Data/DoubleLayerCapacity/Solvation.npz
index c12d668ec75b43cfddd7e557daba0c3aa4e3a049..1b81d1f77a880bb2d4d60acac2f5a69b7d3b9174 100644
Binary files a/examples/Data/DoubleLayerCapacity/Solvation.npz and b/examples/Data/DoubleLayerCapacity/Solvation.npz differ
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/Compressibility.py b/examples/ReproducableCode/DoubleLayerCapacity/Compressibility.py
index 9f1eff086027c9d788b24c1420a02311acdc0925..bf3c2f9d1953f2d1ab9068767e0f643fa71df928 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/Compressibility.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/Compressibility.py
@@ -3,6 +3,8 @@ Jan Habscheid
 Jan.Habscheid@rwth-aachen.de
 
 This script is used to analyze the influence of the compressibility on the charge of the system and the double-layer capacity.
+
+# ! Efforts were made to simplify this similar to the incompressible case. However, this seems not to be robust, yet. Further investigation is needed. This is why it is solved numerically for now.
 '''
 
 # import the src file needed to solve the system of equations
@@ -10,7 +12,7 @@ import sys
 import os
 
 # Add the src directory to the sys.path
-src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../..', 'src')
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../..', 'src')
 sys.path.insert(0, src_path)
 
 from Eq02 import solve_System_2eq
@@ -27,7 +29,7 @@ import numpy as np
 # Define Parameter
 e0 = 1.602e-19 # [As]
 k = 1.381e-23 # [J/K]
-T = 293.15 # [K]
+T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
@@ -38,8 +40,6 @@ LR = 20e-8
 chi = 80 # [-]
 
 # Parameter and bcs for the electrolyte
-a2 = (pR)/(nR_m * k * T)
-K = 'incompressible'
 kappa = 0
 Molarity = 0.01
 y_R = Molarity / nR_mol
@@ -53,13 +53,14 @@ a2 = (pR)/(nR_m * k * T)
 number_cells = 1024
 p_right = 0
 phi_right = 0
+rtol = 1e-4 # ! Change back to 1e-8
 
 # phi^L domain
-Vol_start = 0
+Volt_start = 0
 Volt_end = 0.75
-n_Volts = 10
+n_Volts = 7 # ! Change back to 101
 
-phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
+phi_left_dimless = np.linspace(Volt_start, Volt_end, n_Volts) * e0/(k*T)
 
 y_A, y_C, y_S, phi, p, x, n_sol = [], [], [], [], [], [], []
 for index, K_ in enumerate(K_vec):
@@ -67,18 +68,17 @@ for index, K_ in enumerate(K_vec):
     print('K: ', str(K_))
     print('-------------------------------------------------------------------')
     phi_inner, y_A_inner, y_C_inner, y_S_inner, p_inner,x_inner, n_inner = [], [], [], [], [], [], []
-    for i, phi_bcs in enumerate(phi_left):
+    for i, phi_bcs in enumerate(phi_left_dimless):
         print('-------------------------------------------------------------------')
         print('K: ', str(K_))
         print('phi_bcs: ', phi_bcs)
         print('-------------------------------------------------------------------')
-        # ToDo: relax_param
-        relax_param = 0.3
+        relax_param = 0.05
         if phi_bcs > 20:
-            relax_param = 0.05
+            relax_param = 0.02
         if phi_bcs > 25:
-            relax_param = 0.03
-        y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K_, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style='hard_hard_log', rtol=1e-4, max_iter=10_000, return_type='Vector')
+            relax_param = 0.01
+        y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K_, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style='hard_log', rtol=rtol, max_iter=10_000, return_type='Vector')
         y_A_inner.append(y_A_)
         y_C_inner.append(y_C_)
         y_S_inner.append(1 - y_A_ - y_C_)
@@ -93,7 +93,7 @@ for index, K_ in enumerate(K_vec):
     p.append(p_inner)
     x.append(x_inner)
     n_sol.append(n_inner)
-phi_left *= (k*T)/e0
+phi_left_dim = phi_left_dimless * (k*T)/e0
     
 
 
@@ -101,24 +101,24 @@ phi_left *= (k*T)/e0
 Q = []
 for i in range(len(K_vec)):
     Q_inner = []
-    for j in range(len(phi_left)):
+    for j in range(len(phi_left_dimless)):
         # ToDo: nR_m -> n
         Q_inner.append(Q_num_(y_A[i][j], y_C[i][j], n_sol[i][j], x[i][j]))
     Q.append(Q_inner)
 Q = np.array(Q)
 
 # Double layer capacity
-dx_ = phi_left[1] - phi_left[0] # [1/V], uniformly distributed
+dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], uniformly distributed
 C_DL = (Q[:,1:] - Q[:,:-1])/dx_ # [µAs/cm³] # ? Right unit?
 C_DL = np.array(C_DL)
-phi_left_center = (phi_left[1:] + phi_left[:-1])/2 # Center points for plotting C_DL
+phi_left_center = (phi_left_dimless[1:] + phi_left_dimless[:-1])/2 # Center points for plotting C_DL
 
 
 
 # Visualizations
 # Charge
 plt.figure()
-[plt.plot(phi_left, Q[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
+[plt.plot(phi_left_dimless, Q[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
 plt.grid()
 plt.legend()
 plt.xlabel('$\delta \\varphi$ [-]')
@@ -143,9 +143,9 @@ colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
 plt.figure()
 for i in range(len(K_vec)):
     clr = colors[i]
-    plt.plot(phi_left, y_A_np[i,:,0], '--', color=clr)
-    plt.plot(phi_left, y_C_np[i,:,0], '-', color=clr)
-    plt.plot(phi_left, y_S_np[i,:,0], ':', color=clr)
+    plt.plot(phi_left_dimless, y_A_np[i,:,0], '--', color=clr)
+    plt.plot(phi_left_dimless, y_C_np[i,:,0], '-', color=clr)
+    plt.plot(phi_left_dimless, y_S_np[i,:,0], ':', color=clr)
 dummy, = plt.plot(0, 0, color='grey', linestyle='--', label='$y_A$')
 dummy, = plt.plot(0, 0, color='grey', linestyle='-', label='$y_C$')
 dummy, = plt.plot(0, 0, color='grey', linestyle=':', label='$y_S$')
@@ -164,9 +164,9 @@ plt.figure()
 clr = colors[0]
 for i in range(len(K_vec)):
     clr = colors[i]
-    plt.plot(phi_left, n_A_np[i,:,0], '--', color=clr)
-    plt.plot(phi_left, n_C_np[i,:,0], '-', color=clr)
-    plt.plot(phi_left, n_S_np[i,:,0], ':', color=clr)
+    plt.plot(phi_left_dimless, n_A_np[i,:,0], '--', color=clr)
+    plt.plot(phi_left_dimless, n_C_np[i,:,0], '-', color=clr)
+    plt.plot(phi_left_dimless, n_S_np[i,:,0], ':', color=clr)
 dummy, = plt.plot(0, 0, color='grey', linestyle='--', label='$n_A$')
 dummy, = plt.plot(0, 0, color='grey', linestyle='-', label='$n_C$')
 dummy, = plt.plot(0, 0, color='grey', linestyle=':', label='$n_S$')
@@ -178,3 +178,8 @@ plt.xlabel('$\delta \\varphi$ [-]')
 plt.ylabel('$n_\\alpha^L [-]$')
 plt.tight_layout()
 plt.show()
+
+
+
+# Store data
+np.savez('../../Data/DoubleLayerCapacity/Compressibility.npz', Lambda2=Lambda2, a2=a2, K_vec=K_vec, kappa=kappa, Molarity=Molarity, y_R=y_R, z_A=z_A, z_C=z_C, number_cells=number_cells, p_right=p_right, phi_right=phi_right, rtol=rtol, Volt_start=Volt_start, Volt_end=Volt_end, n_Volts=n_Volts, phi_left_dim=phi_left_dim, phi_left_dimless=phi_left_dimless, y_A=y_A, y_C=y_C, y_S=y_S, phi=phi, p=p, x=x, n_sol=n_sol, Q=Q, C_DL=C_DL, phi_left_center=phi_left_center, y_A_np=y_A_np, y_C_np=y_C_np, y_S_np=y_S_np, n_sol_np=n_sol_np, n_A_np=n_A_np, n_C_np=n_C_np, n_S_np=n_S_np)
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
new file mode 100644
index 0000000000000000000000000000000000000000..cf5abb92597d4273da9257e3de8a446285e43617
--- /dev/null
+++ b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
@@ -0,0 +1,159 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script is used to validate the analytical method for the calculation of the double layer capacity.
+
+# ! Problems, if solvation != 0. Further investigation is needed.
+'''
+
+# import the src file needed to solve the system of equations
+import sys
+import os
+
+# Add the src directory to the sys.path
+src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../../', 'src')
+sys.path.insert(0, src_path)
+
+from Eq02 import solve_System_2eq
+from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
+
+
+# Remove the src directory from sys.path after import
+del sys.path[0]
+
+# Import plotting library
+import matplotlib.pyplot as plt
+import numpy as np
+
+# # Define Parameter
+# e0 = 1.602e-19 # [As]
+# k = 1.381e-23 # [J/K]
+# T = 293.75 # [K]
+
+# # Load numerical data
+# data_numerical = np.load('../../Data/DoubleLayerCapacity/Compressibility.npz')
+# Lambda2 = data_numerical['Lambda2']
+# a2 = data_numerical['a2']
+# K_vec = data_numerical['K_vec']
+# temp = K_vec[1:].astype(np.float64)
+# K_vec_converted = [K_vec[0]]
+# [K_vec_converted.append(temp_) for temp_ in temp]
+# K_vec = K_vec_converted
+# kappa = data_numerical['kappa']
+# Molarity = data_numerical['Molarity']
+# y_R = data_numerical['y_R']
+# z_A = data_numerical['z_A']
+# z_C = data_numerical['z_C']
+# p_right = data_numerical['p_right']
+# Volt_start = data_numerical['Volt_start']
+# Volt_end = data_numerical['Volt_end']
+# n_Volts = data_numerical['n_Volts']
+# phi_left_dim = data_numerical['phi_left_dim']
+# phi_left_dimless = data_numerical['phi_left_dimless']
+# phi_right = data_numerical['phi_right']
+# y_A_num = data_numerical['y_A']
+# y_C_num = data_numerical['y_C']
+# y_S_num = data_numerical['y_S']
+# phi_num = data_numerical['phi']
+# p_num = data_numerical['p']
+# x_num = data_numerical['x']
+# Q_num = data_numerical['Q']
+# C_DL_num = data_numerical['C_DL']
+# phi_left_center = data_numerical['phi_left_center']
+
+# Define Parameter
+e0 = 1.602e-19 # [As]
+k = 1.381e-23 # [J/K]
+T = 293.75 # [K]
+epsilon0 = 8.85e-12 #[F/m]
+F = 9.65e+4 # [As/mol]
+NA = 6.022e+23 # [1/mol] - Avogadro constant
+nR_mol = 55
+nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
+pR = 1.01325 * 1e+5 # [Pa]
+LR = 20e-8
+chi = 80 # [-]
+
+# Parameter and bcs for the electrolyte
+Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
+a2 = (pR)/(nR_m * k * T)
+K = 20_000#'incompressible'
+kappa = 0
+Molarity = 0.01
+y_R = Molarity / nR_mol
+z_A, z_C = -1.0, 1.0
+phi_right = 0.0
+p_right = 0
+
+# Solver settings
+number_cells = 1024
+relax_param = 0.03
+p_right = 0
+
+
+# phi^L domain
+Vol_start = 0
+Volt_end = 0.75
+n_Volts = 10#0
+
+phi_left_dimless = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
+
+
+# Numerical calculations
+# Solution vectors
+y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
+for i, phi_bcs in enumerate(phi_left_dimless):
+    y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
+    y_S_ = 1 - y_A_ - y_C_
+    y_A_num.append(y_A_)
+    y_C_num.append(y_C_)
+    y_S_num.append(y_S_)
+    phi_num.append(phi_)
+    p_num.append(p_)
+    x_num.append(x_)
+    
+Q_num = []
+for j in range(len(phi_left_dimless)):
+    Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
+Q_num = np.array(Q_num)
+
+dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], Assumption: phi^L is uniformly distributed
+C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+C_DL_num = np.array(C_DL_num)
+C_dl_num = C_dl(Q_num, phi_left_dimless)
+
+# Analytical calculations
+# K = K_vec[0]
+# print(K)
+# print(phi_left_dimless)
+Q_ana = Q_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left_dimless, phi_right, p_right, K, Lambda2, a2, kappa)
+C_DL_ana = C_dl(Q_ana, phi_left_dimless)
+
+
+# Plotting
+plt.figure()
+# plt.plot(phi_left, Q_num[1] - Q_ana, label='Difference')
+# plt.plot(phi_left_dimless, Q_num - Q_ana, label='Difference')
+# plt.plot(phi_left_dimless, Q_num[1], label='Numerical')
+plt.plot(phi_left_dimless, Q_num, label='Numerical')
+plt.plot(phi_left_dimless, Q_ana, label='Analytical')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$Q_{num} - Q_{ana} [-]$')
+plt.tight_layout()
+plt.show()   
+
+plt.figure()
+# plt.plot(Phi_pot_center(phi_left), C_DL_num[1] - C_DL_ana, label='Difference')
+# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num - C_DL_ana, label='Difference')
+# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num[1], label='Numerical')
+plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left_dimless), C_DL_ana, label='Analytical')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$C_{dl,num} - C_{dl,ana} [-]$')
+plt.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/ValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
similarity index 73%
rename from examples/ReproducableCode/DoubleLayerCapacity/ValidationAnalyticalMethod.py
rename to examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
index 75f11a76f6ee1664a4949e5333b5e8a6d78647de..f2cfffa6ce66efb6811e721d5cce97ea18a13d14 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/ValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
@@ -3,6 +3,8 @@ Jan Habscheid
 Jan.Habscheid@rwth-aachen.de
 
 This script is used to validate the analytical method for the calculation of the double layer capacity.
+
+# ! Problems, if solvation != 0. Further investigation is needed.
 '''
 
 # import the src file needed to solve the system of equations
@@ -27,7 +29,7 @@ import numpy as np
 # Define Parameter
 e0 = 1.602e-19 # [As]
 k = 1.381e-23 # [J/K]
-T = 293.15 # [K]
+T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
@@ -41,7 +43,7 @@ chi = 80 # [-]
 Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
 a2 = (pR)/(nR_m * k * T)
 K = 'incompressible'
-kappa = 0
+kappa = 5
 Molarity = 0.01
 y_R = Molarity / nR_mol
 z_A, z_C = -1.0, 1.0
@@ -57,7 +59,7 @@ p_right = 0
 # phi^L domain
 Vol_start = 0
 Volt_end = 0.75
-n_Volts = 20#0
+n_Volts = 5#0
 
 phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 
@@ -66,7 +68,7 @@ phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 # Solution vectors
 y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
 for i, phi_bcs in enumerate(phi_left):
-    y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
+    y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-10, max_iter=2500, return_type='Vector', relax_param=relax_param)
     y_S_ = 1 - y_A_ - y_C_
     y_A_num.append(y_A_)
     y_C_num.append(y_C_)
@@ -80,21 +82,23 @@ for j in range(len(phi_left)):
     Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
 Q_num = np.array(Q_num)
 
-dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
-C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
-C_DL_num = np.array(C_DL_num)
+# dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
+# C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+# C_DL_num = np.array(C_DL_num)
 C_dl_num = C_dl(Q_num, phi_left)
 
 
 
 # Analytical calculations
-Q_ana = Q_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, Lambda2, a2, kappa)
+Q_ana = Q_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, K, Lambda2, a2, kappa)
 C_DL_ana = C_dl(Q_ana, phi_left)
 
 
 # Plotting
 plt.figure()
-plt.plot(phi_left, Q_num - Q_ana, label='Difference')
+# plt.plot(phi_left, Q_num - Q_ana, label='Difference')
+plt.plot(phi_left, Q_num, label='Numerical')
+plt.plot(phi_left, Q_ana, label='Analytical')
 plt.grid()
 plt.legend()
 plt.xlabel('$\delta \\varphi$ [-]')
@@ -103,7 +107,9 @@ plt.tight_layout()
 plt.show()   
 
 plt.figure()
-plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
+# plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
+plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left), C_DL_ana, label='Analytical')
 plt.grid()
 plt.legend()
 plt.xlabel('$\delta \\varphi$ [-]')
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/Lambda.py b/examples/ReproducableCode/DoubleLayerCapacity/Lambda.py
index cbb20e75c248ddc3b59c59f04541479202209c75..851295a814c41c954107c7d6b66928bf4ac85599 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/Lambda.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/Lambda.py
@@ -26,7 +26,7 @@ import numpy as np
 # Define Parameter
 e0 = 1.602e-19 # [As]
 k = 1.381e-23 # [J/K]
-T = 293.15 # [K]
+T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
@@ -66,8 +66,8 @@ Q_DL_dimless_ = []
 for Lambda2 in Lambda2_vec:
     y_A_R, y_C_R = Molarity / nR_mol, Molarity / nR_mol
     y_N_R = 1 - y_A_R - y_C_R
-    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
-    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, nR_m, e0, LR, kappa))
     
 C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
 C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/Molarity.py b/examples/ReproducableCode/DoubleLayerCapacity/Molarity.py
index 1ed04415d1b556414466eb844aec6828e2df6849..101b6c889728d38cec9b0e29c89c7f87579d3f39 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/Molarity.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/Molarity.py
@@ -25,7 +25,7 @@ import numpy as np
 # Define Parameter
 e0 = 1.602e-19 # [As]
 k = 1.381e-23 # [J/K]
-T = 293.15 # [K]
+T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
@@ -63,8 +63,8 @@ Q_DL_dimless_ = []
 for mol in Molarity:
     y_A_R, y_C_R = mol / nR_mol, mol / nR_mol
     y_N_R = 1 - y_A_R - y_C_R
-    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
-    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, nR_m, e0, LR, kappa))
     
 C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
 C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/Solvation.py b/examples/ReproducableCode/DoubleLayerCapacity/Solvation.py
index 93b535ff072e2ed959416f485ee2832f328598a4..c2c30ddabcd1d8e0a59bbc00d5b46af141668dc9 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/Solvation.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/Solvation.py
@@ -28,7 +28,7 @@ import numpy as np
 # Define Parameter
 e0 = 1.602e-19 # [As]
 k = 1.381e-23 # [J/K]
-T = 293.15 # [K]
+T = 293.75 # [K]
 epsilon0 = 8.85e-12 #[F/m]
 F = 9.65e+4 # [As/mol]
 NA = 6.022e+23 # [1/mol] - Avogadro constant
@@ -72,8 +72,8 @@ Q_DL_dimless_ = []
 for kappa in kappa_vec:
     y_A_R, y_C_R = Molarity / nR_mol, Molarity / nR_mol
     y_N_R = 1 - y_A_R - y_C_R
-    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, kappa))
-    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, Lambda2, a2, nR_m, e0, LR, kappa))
+    Q_DL_dimless_.append(Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, kappa))
+    Q_DL_dim_.append(Q_DL_dim_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, Phi_Pot_Diff_dimless, phi_R, p_R, K, Lambda2, a2, nR_m, e0, LR, kappa))
     
 C_DL_dim = [C_dl(q_dl, Phi_Pot_Diff_dim) for q_dl in Q_DL_dim_]
 C_DL_dimless = [C_dl(q_dl, Phi_Pot_Diff_dimless) for q_dl in Q_DL_dimless_]
diff --git a/examples/Visualizations/DoubleLayerCapacity.py/VisCompressibility.py b/examples/Visualizations/DoubleLayerCapacity.py/VisCompressibility.py
new file mode 100644
index 0000000000000000000000000000000000000000..7a99e021a7ffa8fdb89de239fc3588f51b66ca38
--- /dev/null
+++ b/examples/Visualizations/DoubleLayerCapacity.py/VisCompressibility.py
@@ -0,0 +1,72 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script visualizes the incompressible solutions for the charge of the systema and the double layer capacity.
+'''
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+# Molarity
+
+# Load data
+data = np.load('../../Data/DoubleLayerCapacity/Compressibility.npz')
+
+Q = data['Q']
+C_DL = data['C_DL']
+y_A = data['y_A']
+y_C = data['y_C']
+y_S = data['y_S']
+phi_left = data['phi_left']
+phi_left_center = data['phi_left_center']
+n_sol = data['n_sol']
+K_vec = data['K_vec']
+
+
+# Charge
+plt.figure()
+[plt.plot(phi_left, Q[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('Q $[-]$')
+plt.tight_layout()
+plt.show()   
+
+# Double layer capacity
+plt.figure()
+[plt.plot(phi_left_center, C_DL[i], label=f'K = {K_vec[i]}') for i in range(len(K_vec))]
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$C_{dl} [-]$')
+plt.tight_layout()
+plt.show()
+
+# y_alpha_L
+y_A_np, y_C_np, y_S_np = np.array(y_A), np.array(y_C), np.array(y_S)
+markers = ['--', '-', ':']
+colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple']
+
+n_sol_np = np.array(n_sol)
+n_A_np, n_C_np, n_S_np = y_A_np * n_sol_np, y_C_np * n_sol_np, y_S_np * n_sol_np
+plt.figure()
+clr = colors[0]
+for i in range(len(K_vec)):
+    clr = colors[i]
+    plt.plot(phi_left, n_A_np[i,:,0], '--', color=clr)
+    plt.plot(phi_left, n_C_np[i,:,0], '-', color=clr)
+    plt.plot(phi_left, n_S_np[i,:,0], ':', color=clr)
+dummy, = plt.plot(0, 0, color='grey', linestyle='--', label='$n_A$')
+dummy, = plt.plot(0, 0, color='grey', linestyle='-', label='$n_C$')
+dummy, = plt.plot(0, 0, color='grey', linestyle=':', label='$n_S$')
+for index, K_ in enumerate(K_vec):
+    dummy, = plt.plot(0, 0, color=colors[index], linestyle='-', label=f'K = {K_}')
+plt.grid()
+plt.legend()
+plt.xlabel('$\delta \\varphi$ [-]')
+plt.ylabel('$n_\\alpha^L [-]$')
+plt.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/src/Helper_DoubleLayerCapacity.py b/src/Helper_DoubleLayerCapacity.py
index 50df4ba64fd2b9919a679efa94794033cfd901bd..d8bd6970498e7bb1f73f9a367cfd10f903346359 100644
--- a/src/Helper_DoubleLayerCapacity.py
+++ b/src/Helper_DoubleLayerCapacity.py
@@ -110,7 +110,7 @@ def Q_num_(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float
     nF_int = -np.trapz(nF_dimensionless, x)
     return nF_int
 
-def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, Lambda2:float, a2:float, kappa:float) -> float:
+def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
     '''
     Calculates charge of the system using the analytical method in dimensionless units
 
@@ -132,16 +132,18 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
     z_C : float
         Charge number of cations
     phi_L : float
-        Electric potential at the left boundary
+        Electric potential at the left boundary [V]
     phi_R : float
         Electric potential at the right boundary
     p_R : float
         Pressure at the right boundary
+    K : str | float
+        Bulk modulus, use 'incompressible' for an incompressible mixture
     Lambda2 : float
         Dimensionless parameter
     a2 : float
         Dimensionless parameter
-    kappa : float
+    solvation : float
         Solvation number
 
     Returns
@@ -150,17 +152,34 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
         Charge of the system in dimensionless units
     '''
     z_N = 0
-    D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
-    D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
-    D_N = y_N_R / (np.exp(-a2*p_R-z_N*phi_R))
-    E_p = p_R
-    CLambda_L = np.log(D_A * np.exp(-z_A * phi_L + E_p * a2) + D_C * np.exp(-z_C * phi_L + E_p * a2) + D_N * np.exp(-z_N * phi_L + E_p * a2))
-    CLambda_R = np.log(D_A * np.exp(-z_A * phi_R + E_p * a2) + D_C * np.exp(-z_C * phi_R + E_p * a2) + D_N * np.exp(-z_N * phi_R + E_p * a2))
-    Lambda = np.sqrt(Lambda2)
-    Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(kappa+1)) * (np.sqrt(CLambda_L) - np.sqrt(CLambda_R))
+    if K == 'incompressible':
+        D_A = y_A_R / (np.exp(-(solvation+1)*a2*p_R-z_A*phi_R))
+        D_C = y_C_R / (np.exp(-(solvation+1)*a2*p_R-z_C*phi_R))
+        D_N = y_N_R / (np.exp(-(solvation+1)*a2*p_R-z_N*phi_R))
+        E_p = p_R
+        CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - (solvation+1) * E_p) + D_C * np.exp(-z_C * phi_L - (solvation+1) * E_p) + D_N * np.exp(-z_N * phi_L - (solvation+1) * E_p))
+        CLambda_R = np.log(D_A * np.exp(-z_A * phi_R + E_p * a2) + D_C * np.exp(-z_C * phi_R + E_p * a2) + D_N * np.exp(-z_N * phi_R + E_p * a2))
+        Lambda = np.sqrt(Lambda2)
+        Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L) - np.sqrt(CLambda_R))
+    else:
+        C_A = y_A_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_A*phi_R))
+        C_C = y_C_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_C*phi_R))
+        C_N = y_N_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_N*phi_R))
+        E_p = p_R
+
+        Lambda_tilda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
+        Lambda_tilda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+
+        Lambda_hat_L = np.exp(np.log(Lambda_tilda_L) / (a2 * K))
+        Lambda_hat_R = np.exp(np.log(Lambda_tilda_R) / (a2 * K))
+
+        dx_Phi_L = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_L))
+        dx_Phi_R = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_R))
+
+        Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
     return Q_DL
 
-def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, kappa:float) -> float:
+def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float:
     '''
     Calculates charge of the system using the analytical method in dimensionless units
 
@@ -187,6 +206,8 @@ def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, ph
         Value of electric potential at the right
     p_R : float
         Value of pressure at the right boundary
+    K : str | float
+        Bulk modulus, use 'incompressible' for an incompress
     Lambda2 : float
         Dimensionless parameter
     a2 : float
@@ -197,7 +218,7 @@ def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, ph
         Dielectric constant
     LR : float
         Reference length in m
-    kappa : float
+    solvation : float
         Solvation number
 
     Returns
@@ -205,7 +226,7 @@ def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, ph
     float
         Charge of the system in µAs/cm³
     '''
-    Q_DL = Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, Lambda2, a2, kappa)
+    Q_DL = Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
     Q_DL *= nR_m * e0 * LR
     Q_DL *= 1e+6 
     Q_DL *= 1/(1e+4)
diff --git a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc
index 278a45c3840a132f67cf7d53ce418765cdeb8a48..0520aeca0b59a17e6054ae07e36af289020f6477 100644
Binary files a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc and b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc differ