diff --git a/README.md b/README.md
index 05c83c73c0974ada4bd55a343470532594afcf07..0067dace29cb3ac2028e1c91a5ced5e4d531a620 100644
--- a/README.md
+++ b/README.md
@@ -27,7 +27,7 @@ All the calculations were performed on a Linux machine. According to the documen
 ```
 conda create --name fenicsx-env python=3.12.3 -y
 conda activate fenicsx-env
-conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 -y
+conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 pytest==8.3.3 -y
 ```
 
 ### Alternative installation
@@ -45,6 +45,16 @@ docker compose build
 docker compose run solver
 ```
 
+### Testing
+
+Use pytest with 
+```
+python -m pytest
+```
+
+to verify that everything was installed correctly.
+
+
 ## Usage
 
 Find the visualizations from the thesis and some extra calculations in the "examples" folder.
diff --git a/environment.yml b/environment.yml
index 4c54a3639c4b9d7a7ec0829c6ec66d6eb1ffced1..c5ec78808a6c5a0ade6b7d19e74f92688e4b3da8 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,320 +1,85 @@
-name: fenicsx
-channels:
-  - conda-forge
-  - defaults
-dependencies:
-  - _libgcc_mutex=0.1=conda_forge
-  - _openmp_mutex=4.5=2_gnu
-  - aiohttp=3.9.5=py312h98912ed_0
-  - aiosignal=1.3.1=pyhd8ed1ab_0
-  - alsa-lib=1.2.12=h4ab18f5_0
-  - aom=3.9.1=hac33072_0
-  - asttokens=2.4.1=pyhd8ed1ab_0
-  - attr=2.5.1=h166bdaf_1
-  - attrs=23.2.0=pyh71513ae_0
-  - binutils_impl_linux-64=2.40=ha1999f0_7
-  - binutils_linux-64=2.40=hb3c18ed_9
-  - blosc=1.21.5=hc2324a3_1
-  - brotli=1.1.0=hd590300_1
-  - brotli-bin=1.1.0=hd590300_1
-  - brotli-python=1.1.0=py312h30efb56_1
-  - bzip2=1.0.8=hd590300_5
-  - c-ares=1.28.1=hd590300_0
-  - c-blosc2=2.14.4=hb4ffafa_1
-  - ca-certificates=2024.6.2=hbcca054_0
-  - cairo=1.18.0=h3faef2a_0
-  - certifi=2024.6.2=pyhd8ed1ab_0
-  - cffi=1.16.0=py312hf06ca03_0
-  - charset-normalizer=3.3.2=pyhd8ed1ab_0
-  - comm=0.2.2=pyhd8ed1ab_0
-  - contourpy=1.2.1=py312h8572e83_0
-  - cycler=0.12.1=pyhd8ed1ab_0
-  - dav1d=1.2.1=hd590300_0
-  - dbus=1.13.6=h5008d03_3
-  - debugpy=1.8.1=py312h30efb56_0
-  - decorator=5.1.1=pyhd8ed1ab_0
-  - double-conversion=3.3.0=h59595ed_0
-  - eigen=3.4.0=h00ab1b0_0
-  - exceptiongroup=1.2.0=pyhd8ed1ab_2
-  - executing=2.0.1=pyhd8ed1ab_0
-  - expat=2.6.2=h59595ed_0
-  - fenics-basix=0.8.0=py312h2492b07_1
-  - fenics-dolfinx=0.8.0=py312hb4760eb_103
-  - fenics-ffcx=0.8.0=pyh4af843d_0
-  - fenics-libbasix=0.8.0=h9187eef_1
-  - fenics-libdolfinx=0.8.0=h714c792_103
-  - fenics-ufcx=0.8.0=h22f594c_0
-  - fenics-ufl=2024.1.0=pyhd8ed1ab_0
-  - ffmpeg=6.1.1=gpl_he44c6f3_112
-  - fftw=3.3.10=mpi_mpich_hbcf76dd_10
-  - font-ttf-dejavu-sans-mono=2.37=hab24e00_0
-  - font-ttf-inconsolata=3.000=h77eed37_0
-  - font-ttf-source-code-pro=2.038=h77eed37_0
-  - font-ttf-ubuntu=0.83=h77eed37_2
-  - fontconfig=2.14.2=h14ed4e7_0
-  - fonts-conda-ecosystem=1=0
-  - fonts-conda-forge=1=0
-  - fonttools=4.53.0=py312h9a8786e_0
-  - freetype=2.12.1=h267a509_2
-  - fribidi=1.0.10=h36c2ea0_0
-  - frozenlist=1.4.1=py312h98912ed_0
-  - gcc_impl_linux-64=12.3.0=h58ffeeb_11
-  - gcc_linux-64=12.3.0=h9528a6a_9
-  - gettext=0.22.5=h59595ed_2
-  - gettext-tools=0.22.5=h59595ed_2
-  - gl2ps=1.4.2=h0708190_0
-  - glew=2.1.0=h9c3ff4c_2
-  - glib=2.80.2=hf974151_0
-  - glib-tools=2.80.2=hb6ce0ca_0
-  - gmp=6.3.0=hac33072_2
-  - gnutls=3.7.9=hb077bed_0
-  - graphite2=1.3.13=h59595ed_1003
-  - gst-plugins-base=1.24.4=h9ad1361_0
-  - gstreamer=1.24.4=haf2f30d_0
-  - gxx_impl_linux-64=12.3.0=h2a574ab_11
-  - gxx_linux-64=12.3.0=ha28b414_9
-  - harfbuzz=8.5.0=hfac3d4d_0
-  - hdf4=4.2.15=h2a13503_7
-  - hdf5=1.14.3=mpi_mpich_h0f54ddc_5
-  - hypre=2.31.0=mpi_mpich_hd1da18f_1
-  - icu=73.2=h59595ed_0
-  - idna=3.7=pyhd8ed1ab_0
-  - importlib-metadata=7.2.0=pyha770c72_0
-  - importlib_metadata=7.2.0=hd8ed1ab_0
-  - ipykernel=6.29.4=pyh3099207_0
-  - ipython=8.25.0=pyh707e725_0
-  - jedi=0.19.1=pyhd8ed1ab_0
-  - jsoncpp=1.9.5=h4bd325d_1
-  - jupyter_client=8.6.2=pyhd8ed1ab_0
-  - jupyter_core=5.7.2=py312h7900ff3_0
-  - kahip=3.16=h6a42626_3
-  - kahip-python=3.16=py312had6fb94_3
-  - kernel-headers_linux-64=2.6.32=he073ed8_17
-  - keyutils=1.6.1=h166bdaf_0
-  - kiwisolver=1.4.5=py312h8572e83_1
-  - krb5=1.21.2=h659d440_0
-  - lame=3.100=h166bdaf_1003
-  - lcms2=2.16=hb7c19ff_0
-  - ld_impl_linux-64=2.40=hf3520f5_7
-  - lerc=4.0.0=h27087fc_0
-  - libabseil=20240116.2=cxx17_h59595ed_0
-  - libadios2=2.10.0=mpi_mpich_h3e60829_3
-  - libaec=1.1.3=h59595ed_0
-  - libasprintf=0.22.5=h661eb56_2
-  - libasprintf-devel=0.22.5=h661eb56_2
-  - libass=0.17.1=h8fe9dca_1
-  - libblas=3.9.0=22_linux64_openblas
-  - libboost=1.84.0=hba137d9_3
-  - libboost-devel=1.84.0=h00ab1b0_3
-  - libboost-headers=1.84.0=ha770c72_3
-  - libbrotlicommon=1.1.0=hd590300_1
-  - libbrotlidec=1.1.0=hd590300_1
-  - libbrotlienc=1.1.0=hd590300_1
-  - libcap=2.69=h0f662aa_0
-  - libcblas=3.9.0=22_linux64_openblas
-  - libclang-cpp18.1=18.1.7=default_h9bb3924_0
-  - libclang13=18.1.7=default_h087397f_0
-  - libcups=2.3.3=h4637d8d_4
-  - libcurl=8.8.0=hca28451_0
-  - libdeflate=1.20=hd590300_0
-  - libdrm=2.4.121=h4ab18f5_0
-  - libedit=3.1.20191231=he28a2e2_2
-  - libev=4.33=hd590300_2
-  - libexpat=2.6.2=h59595ed_0
-  - libffi=3.4.2=h7f98852_5
-  - libflac=1.4.3=h59595ed_0
-  - libgcc-devel_linux-64=12.3.0=h6b66f73_111
-  - libgcc-ng=13.2.0=h77fa898_11
-  - libgcrypt=1.10.3=hd590300_0
-  - libgettextpo=0.22.5=h59595ed_2
-  - libgettextpo-devel=0.22.5=h59595ed_2
-  - libgfortran-ng=13.2.0=h69a702a_11
-  - libgfortran5=13.2.0=h3d2ce59_11
-  - libglib=2.80.2=hf974151_0
-  - libglu=9.0.0=hac7e632_1003
-  - libgomp=13.2.0=h77fa898_11
-  - libgpg-error=1.49=h4f305b6_0
-  - libhwloc=2.10.0=default_h5622ce7_1001
-  - libiconv=1.17=hd590300_2
-  - libidn2=2.3.7=hd590300_0
-  - libjpeg-turbo=3.0.0=hd590300_1
-  - liblapack=3.9.0=22_linux64_openblas
-  - libllvm18=18.1.7=hb77312f_0
-  - libnetcdf=4.9.2=nompi_h135f659_114
-  - libnghttp2=1.58.0=h47da74e_1
-  - libnsl=2.0.1=hd590300_0
-  - libogg=1.3.4=h7f98852_1
-  - libopenblas=0.3.27=pthreads_h413a1c8_0
-  - libopenvino=2024.1.0=h2da1b83_7
-  - libopenvino-auto-batch-plugin=2024.1.0=hb045406_7
-  - libopenvino-auto-plugin=2024.1.0=hb045406_7
-  - libopenvino-hetero-plugin=2024.1.0=h5c03a75_7
-  - libopenvino-intel-cpu-plugin=2024.1.0=h2da1b83_7
-  - libopenvino-intel-gpu-plugin=2024.1.0=h2da1b83_7
-  - libopenvino-intel-npu-plugin=2024.1.0=he02047a_7
-  - libopenvino-ir-frontend=2024.1.0=h5c03a75_7
-  - libopenvino-onnx-frontend=2024.1.0=h07e8aee_7
-  - libopenvino-paddle-frontend=2024.1.0=h07e8aee_7
-  - libopenvino-pytorch-frontend=2024.1.0=he02047a_7
-  - libopenvino-tensorflow-frontend=2024.1.0=h39126c6_7
-  - libopenvino-tensorflow-lite-frontend=2024.1.0=he02047a_7
-  - libopus=1.3.1=h7f98852_1
-  - libpciaccess=0.18=hd590300_0
-  - libpng=1.6.43=h2797004_0
-  - libpq=16.3=ha72fbe1_0
-  - libprotobuf=4.25.3=h08a7969_0
-  - libptscotch=7.0.4=h2376d02_5
-  - libsanitizer=12.3.0=hb8811af_11
-  - libscotch=7.0.4=h3055ed5_5
-  - libsndfile=1.2.2=hc60ed4a_1
-  - libsodium=1.0.18=h36c2ea0_1
-  - libsqlite=3.46.0=hde9e2c9_0
-  - libssh2=1.11.0=h0841786_0
-  - libstdcxx-devel_linux-64=12.3.0=h6b66f73_111
-  - libstdcxx-ng=13.2.0=hc0a3c3a_11
-  - libsystemd0=255=h3516f8a_1
-  - libtasn1=4.19.0=h166bdaf_0
-  - libtheora=1.1.1=h7f98852_1005
-  - libtiff=4.6.0=h1dd3fc0_3
-  - libunistring=0.9.10=h7f98852_0
-  - libuuid=2.38.1=h0b41bf4_0
-  - libva=2.21.0=h4ab18f5_2
-  - libvorbis=1.3.7=h9c3ff4c_0
-  - libvpx=1.14.1=hac33072_0
-  - libwebp-base=1.4.0=hd590300_0
-  - libxcb=1.15=h0b41bf4_0
-  - libxcrypt=4.4.36=hd590300_1
-  - libxkbcommon=1.7.0=h662e7e4_0
-  - libxml2=2.12.7=hc051c1a_1
-  - libzip=1.10.1=h2629f0a_3
-  - libzlib=1.2.13=h4ab18f5_6
-  - loguru=0.7.2=py312h7900ff3_1
-  - lz4-c=1.9.4=hcb278e6_0
-  - matplotlib-base=3.8.4=py312h20ab3a6_2
-  - matplotlib-inline=0.1.7=pyhd8ed1ab_0
-  - metis=5.1.0=h59595ed_1007
-  - mpfr=4.2.1=h9458935_1
-  - mpg123=1.32.6=h59595ed_0
-  - mpi=1.0=mpich
-  - mpi4py=3.1.6=py312h60e9011_1
-  - mpich=4.2.1=h63d650b_101
-  - msgpack-python=1.0.8=py312h2492b07_0
-  - multidict=6.0.5=py312h98912ed_0
-  - mumps-include=5.7.2=ha770c72_0
-  - mumps-mpi=5.7.2=h09c71e5_0
-  - munkres=1.1.4=pyh9f0ad1d_0
-  - mysql-common=8.3.0=hf1915f5_4
-  - mysql-libs=8.3.0=hca2cd23_4
-  - ncurses=6.5=h59595ed_0
-  - nest-asyncio=1.6.0=pyhd8ed1ab_0
-  - nettle=3.9.1=h7ab15ed_0
-  - nlohmann_json=3.11.3=h59595ed_0
-  - numpy=1.26.4=py312heda63a1_0
-  - ocl-icd=2.3.2=hd590300_1
-  - openh264=2.4.1=h59595ed_0
-  - openjpeg=2.5.2=h488ebb8_0
-  - openssl=3.3.1=h4ab18f5_0
-  - p11-kit=0.24.1=hc5aa10d_0
-  - packaging=24.1=pyhd8ed1ab_0
-  - parmetis=4.0.3=h583469f_1006
-  - parso=0.8.4=pyhd8ed1ab_0
-  - pcre2=10.43=hcad00b1_0
-  - petsc=3.21.2=real_h7906ff3_102
-  - petsc4py=3.21.2=py312hf1b966b_1
-  - pexpect=4.9.0=pyhd8ed1ab_0
-  - pickleshare=0.7.5=py_1003
-  - pillow=10.3.0=py312hdcec9eb_0
-  - pip=24.0=pyhd8ed1ab_0
-  - pixman=0.43.2=h59595ed_0
-  - pkg-config=0.29.2=h36c2ea0_1008
-  - platformdirs=4.2.2=pyhd8ed1ab_0
-  - pooch=1.8.2=pyhd8ed1ab_0
-  - proj=9.3.1=h1d62c97_0
-  - prompt-toolkit=3.0.47=pyha770c72_0
-  - psutil=5.9.8=py312h98912ed_0
-  - pthread-stubs=0.4=h36c2ea0_1001
-  - ptscotch=7.0.4=h23d43cc_5
-  - ptyprocess=0.7.0=pyhd3deb0d_0
-  - pugixml=1.14=h59595ed_0
-  - pulseaudio-client=17.0=hb77b528_0
-  - pure_eval=0.2.2=pyhd8ed1ab_0
-  - pycparser=2.22=pyhd8ed1ab_0
-  - pygments=2.18.0=pyhd8ed1ab_0
-  - pyparsing=3.1.2=pyhd8ed1ab_0
-  - pysocks=1.7.1=pyha2e5f31_6
-  - python=3.12.3=hab00c5b_0_cpython
-  - python-dateutil=2.9.0=pyhd8ed1ab_0
-  - python_abi=3.12=4_cp312
-  - pyvista=0.43.10=pyhd8ed1ab_0
-  - pyzmq=26.0.3=py312h8fd38d8_0
-  - qt6-main=6.7.1=h2471661_2
-  - readline=8.2=h8228510_1
-  - requests=2.32.3=pyhd8ed1ab_0
-  - scalapack=2.2.0=h417d24c_2
-  - scooby=0.10.0=pyhd8ed1ab_0
-  - scotch=7.0.4=h23d43cc_5
-  - setuptools=70.1.0=pyhd8ed1ab_0
-  - six=1.16.0=pyh6c4a22f_0
-  - slepc=3.21.1=real_h97ad6bc_302
-  - slepc4py=3.21.1=py312h68801f8_102
-  - snappy=1.2.0=hdb0a2a9_1
-  - sqlite=3.46.0=h6d4b2fc_0
-  - stack_data=0.6.2=pyhd8ed1ab_0
-  - suitesparse=7.7.0=hf4753ba_1
-  - superlu=5.2.2=h00795ac_0
-  - superlu_dist=9.0.0=h3feb4ed_1
-  - svt-av1=2.1.0=hac33072_0
-  - sysroot_linux-64=2.12=he073ed8_17
-  - tbb=2021.12.0=h297d8ca_1
-  - tbb-devel=2021.12.0=h7c56ddd_1
-  - tk=8.6.13=noxft_h4845f30_101
-  - tornado=6.4.1=py312h9a8786e_0
-  - traitlets=5.14.3=pyhd8ed1ab_0
-  - typing_extensions=4.12.2=pyha770c72_0
-  - urllib3=2.2.2=pyhd8ed1ab_0
-  - utfcpp=4.0.5=ha770c72_0
-  - vtk=9.3.0=qt_py312h1234567_200
-  - vtk-base=9.3.0=qt_py312h1234567_200
-  - vtk-io-ffmpeg=9.3.0=qt_py312h1234567_200
-  - wayland=1.23.0=h5291e77_0
-  - wcwidth=0.2.13=pyhd8ed1ab_0
-  - wheel=0.43.0=pyhd8ed1ab_1
-  - wslink=2.1.1=pyhd8ed1ab_0
-  - x264=1!164.3095=h166bdaf_2
-  - x265=3.5=h924138e_3
-  - xcb-util=0.4.0=hd590300_1
-  - xcb-util-cursor=0.1.4=hd590300_1
-  - xcb-util-image=0.4.0=h8ee46fc_1
-  - xcb-util-keysyms=0.4.0=h8ee46fc_1
-  - xcb-util-renderutil=0.3.9=hd590300_1
-  - xcb-util-wm=0.4.1=h8ee46fc_1
-  - xkeyboard-config=2.42=h4ab18f5_0
-  - xorg-fixesproto=5.0=h7f98852_1002
-  - xorg-kbproto=1.0.7=h7f98852_1002
-  - xorg-libice=1.1.1=hd590300_0
-  - xorg-libsm=1.2.4=h7391055_0
-  - xorg-libx11=1.8.9=h8ee46fc_0
-  - xorg-libxau=1.0.11=hd590300_0
-  - xorg-libxdmcp=1.1.3=h7f98852_0
-  - xorg-libxext=1.3.4=h0b41bf4_2
-  - xorg-libxfixes=5.0.3=h7f98852_1004
-  - xorg-libxrender=0.9.11=hd590300_0
-  - xorg-libxt=1.3.0=hd590300_1
-  - xorg-renderproto=0.11.1=h7f98852_1002
-  - xorg-xextproto=7.3.0=h0b41bf4_1003
-  - xorg-xproto=7.0.31=h7f98852_1007
-  - xz=5.2.6=h166bdaf_0
-  - yaml=0.2.5=h7f98852_2
-  - yarl=1.9.4=py312h98912ed_0
-  - zeromq=4.3.5=h75354e8_4
-  - zfp=0.5.5=h9c3ff4c_8
-  - zipp=3.19.2=pyhd8ed1ab_0
-  - zlib=1.2.13=h4ab18f5_6
-  - zlib-ng=2.0.7=h0b41bf4_0
-  - zstd=1.5.6=ha6fb4c9_0
-  - pip:
-      - pandas==2.2.2
-      - pytz==2024.1
-      - scipy==1.14.0
-      - tzdata==2024.1
-prefix: /home/janha/anaconda3/envs/fenicsx
+aiohttp @ file:///home/conda/feedstock_root/build_artifacts/aiohttp_1713964820544/work
+aiosignal @ file:///home/conda/feedstock_root/build_artifacts/aiosignal_1667935791922/work
+asttokens @ file:///home/conda/feedstock_root/build_artifacts/asttokens_1698341106958/work
+attrs @ file:///home/conda/feedstock_root/build_artifacts/attrs_1704011227531/work
+Brotli @ file:///home/conda/feedstock_root/build_artifacts/brotli-split_1695989787169/work
+certifi @ file:///home/conda/feedstock_root/build_artifacts/certifi_1720457958366/work/certifi
+cffi @ file:///home/conda/feedstock_root/build_artifacts/cffi_1696001721842/work
+charset-normalizer @ file:///home/conda/feedstock_root/build_artifacts/charset-normalizer_1698833585322/work
+comm @ file:///home/conda/feedstock_root/build_artifacts/comm_1710320294760/work
+contourpy @ file:///home/conda/feedstock_root/build_artifacts/contourpy_1712429918028/work
+cycler @ file:///home/conda/feedstock_root/build_artifacts/cycler_1696677705766/work
+debugpy @ file:///home/conda/feedstock_root/build_artifacts/debugpy_1722923741382/work
+decorator @ file:///home/conda/feedstock_root/build_artifacts/decorator_1641555617451/work
+distlib==0.3.8
+exceptiongroup @ file:///home/conda/feedstock_root/build_artifacts/exceptiongroup_1720869315914/work
+executing @ file:///home/conda/feedstock_root/build_artifacts/executing_1698579936712/work
+fenics-basix @ file:///home/conda/feedstock_root/build_artifacts/fenics-basix-meta_1716296365445/work/python
+fenics-dolfinx @ file:///home/conda/feedstock_root/build_artifacts/fenics-dolfinx-split_1715597584729/work/python
+fenics-ffcx @ file:///home/conda/feedstock_root/build_artifacts/fenics-ffcx-ufcx_1714058897751/work
+fenics-ufl @ file:///home/conda/feedstock_root/build_artifacts/fenics-ufl_1714027187309/work
+filelock==3.16.0
+fonttools @ file:///home/conda/feedstock_root/build_artifacts/fonttools_1717209197958/work
+frozenlist @ file:///home/conda/feedstock_root/build_artifacts/frozenlist_1702645449276/work
+idna @ file:///home/conda/feedstock_root/build_artifacts/idna_1713279365350/work
+importlib_metadata @ file:///home/conda/feedstock_root/build_artifacts/importlib-metadata_1724187233579/work
+iniconfig==2.0.0
+ipykernel @ file:///home/conda/feedstock_root/build_artifacts/ipykernel_1719845459717/work
+ipython @ file:///home/conda/feedstock_root/build_artifacts/ipython_1719582526268/work
+jedi @ file:///home/conda/feedstock_root/build_artifacts/jedi_1696326070614/work
+jupyter_client @ file:///home/conda/feedstock_root/build_artifacts/jupyter_client_1716472197302/work
+jupyter_core @ file:///home/conda/feedstock_root/build_artifacts/jupyter_core_1710257406420/work
+kiwisolver @ file:///home/conda/feedstock_root/build_artifacts/kiwisolver_1695379925569/work
+loguru @ file:///home/conda/feedstock_root/build_artifacts/loguru_1695547342418/work
+matplotlib @ file:///home/conda/feedstock_root/build_artifacts/matplotlib-suite_1715976243782/work
+matplotlib-inline @ file:///home/conda/feedstock_root/build_artifacts/matplotlib-inline_1713250518406/work
+mdit-py-plugins==0.4.2
+mpi4py @ file:///home/conda/feedstock_root/build_artifacts/mpi4py_1716628219734/work
+msgpack @ file:///home/conda/feedstock_root/build_artifacts/msgpack-python_1715670634347/work
+multidict @ file:///home/conda/feedstock_root/build_artifacts/multidict_1707040729125/work
+munkres==1.1.4
+myst-parser==4.0.0
+nest_asyncio @ file:///home/conda/feedstock_root/build_artifacts/nest-asyncio_1705850609492/work
+numpy @ file:///home/conda/feedstock_root/build_artifacts/numpy_1707225359967/work/dist/numpy-1.26.4-cp312-cp312-linux_x86_64.whl#sha256=031b7d6b2e5e604d9e21fc21be713ebf28ce133ec872dce6de006742d5e49bab
+packaging @ file:///home/conda/feedstock_root/build_artifacts/packaging_1718189413536/work
+parso @ file:///home/conda/feedstock_root/build_artifacts/parso_1712320355065/work
+petsc4py @ file:///home/conda/feedstock_root/build_artifacts/petsc4py_1718362402291/work
+pexpect @ file:///home/conda/feedstock_root/build_artifacts/pexpect_1706113125309/work
+pickleshare @ file:///home/conda/feedstock_root/build_artifacts/pickleshare_1602536217715/work
+pillow @ file:///home/conda/feedstock_root/build_artifacts/pillow_1718833743537/work
+platformdirs @ file:///home/conda/feedstock_root/build_artifacts/platformdirs_1715777629804/work
+pluggy==1.5.0
+pooch @ file:///home/conda/feedstock_root/build_artifacts/pooch_1717777836653/work
+prompt_toolkit @ file:///home/conda/feedstock_root/build_artifacts/prompt-toolkit_1718047967974/work
+psutil @ file:///home/conda/feedstock_root/build_artifacts/psutil_1719274590885/work
+ptyprocess @ file:///home/conda/feedstock_root/build_artifacts/ptyprocess_1609419310487/work/dist/ptyprocess-0.7.0-py2.py3-none-any.whl
+pure_eval @ file:///home/conda/feedstock_root/build_artifacts/pure_eval_1721585709575/work
+pycparser @ file:///home/conda/feedstock_root/build_artifacts/pycparser_1711811537435/work
+Pygments @ file:///home/conda/feedstock_root/build_artifacts/pygments_1714846767233/work
+pyparsing @ file:///home/conda/feedstock_root/build_artifacts/pyparsing_1709721012883/work
+PySocks @ file:///home/conda/feedstock_root/build_artifacts/pysocks_1661604839144/work
+pytest==8.3.3
+python-dateutil @ file:///home/conda/feedstock_root/build_artifacts/python-dateutil_1709299778482/work
+pyvista @ file:///home/conda/feedstock_root/build_artifacts/pyvista_1718665741203/work
+pyzmq @ file:///home/conda/feedstock_root/build_artifacts/pyzmq_1724088076888/work
+requests @ file:///home/conda/feedstock_root/build_artifacts/requests_1717057054362/work
+scipy==1.14.1
+scooby @ file:///home/conda/feedstock_root/build_artifacts/scooby_1714897440316/work
+setuptools==70.1.0
+six @ file:///home/conda/feedstock_root/build_artifacts/six_1620240208055/work
+slepc4py @ file:///home/conda/feedstock_root/build_artifacts/slepc4py_1718913137780/work
+sphinx-copybutton==0.5.2
+sphinx-rtd-theme==2.0.0
+sphinxcontrib-jquery==4.1
+stack-data @ file:///home/conda/feedstock_root/build_artifacts/stack_data_1669632077133/work
+tornado @ file:///home/conda/feedstock_root/build_artifacts/tornado_1717722879218/work
+traitlets @ file:///home/conda/feedstock_root/build_artifacts/traitlets_1713535121073/work
+typing_extensions @ file:///home/conda/feedstock_root/build_artifacts/typing_extensions_1717802530399/work
+urllib3 @ file:///home/conda/feedstock_root/build_artifacts/urllib3_1718653059220/work
+virtualenv==20.26.4
+vtk==9.3.0
+wcwidth @ file:///home/conda/feedstock_root/build_artifacts/wcwidth_1704731205417/work
+wheel==0.43.0
+wslink @ file:///home/conda/feedstock_root/build_artifacts/wslink_1718909657344/work
+yarl @ file:///home/conda/feedstock_root/build_artifacts/yarl_1705508298164/work
+zipp @ file:///home/conda/feedstock_root/build_artifacts/zipp_1723591248676/work
diff --git a/examples/ReproducableCode/__pycache__/TernaryElectrolyte.cpython-312.pyc b/examples/ReproducableCode/__pycache__/TernaryElectrolyte.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..a80ab9785688a27caa0273ead29eea11d02a698f
Binary files /dev/null and b/examples/ReproducableCode/__pycache__/TernaryElectrolyte.cpython-312.pyc differ
diff --git a/examples/ReproducableCode/__pycache__/__init__.cpython-312.pyc b/examples/ReproducableCode/__pycache__/__init__.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..ab8cda1f924e7ee83dade32b42603274c1ab7693
Binary files /dev/null and b/examples/ReproducableCode/__pycache__/__init__.cpython-312.pyc differ
diff --git a/examples/__pycache__/__init__.cpython-312.pyc b/examples/__pycache__/__init__.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..95d74194526c61fbd40b89c98c402b5e5d7edd55
Binary files /dev/null and b/examples/__pycache__/__init__.cpython-312.pyc differ
diff --git a/src/__pycache__/Eq04.cpython-312.pyc b/src/__pycache__/Eq04.cpython-312.pyc
index 9cc7e7b2191c75121f413d3823a6721a84e10309..372244f1b18580495c00297122b6ba8e8c7e8296 100644
Binary files a/src/__pycache__/Eq04.cpython-312.pyc and b/src/__pycache__/Eq04.cpython-312.pyc differ
diff --git a/src/__pycache__/__init__.cpython-312.pyc b/src/__pycache__/__init__.cpython-312.pyc
index a4c2f1275d9e6807ccc319935c53511bffbf67a2..19664c3690d8a063333fd239f1546b1f8295d536 100644
Binary files a/src/__pycache__/__init__.cpython-312.pyc and b/src/__pycache__/__init__.cpython-312.pyc differ
diff --git a/tests/TestData/DLKap.npz b/tests/TestData/DLKap.npz
new file mode 100644
index 0000000000000000000000000000000000000000..dd1008368504e29fc2eb5889716b3989ec93c15f
Binary files /dev/null and b/tests/TestData/DLKap.npz differ
diff --git a/tests/TestData/ElectrolyticDiode.npz b/tests/TestData/ElectrolyticDiode.npz
new file mode 100644
index 0000000000000000000000000000000000000000..83cca5ebb0221af55b3322fcc07d907c1f9c1e82
Binary files /dev/null and b/tests/TestData/ElectrolyticDiode.npz differ
diff --git a/tests/TestData/Eq02.npz b/tests/TestData/Eq02.npz
new file mode 100644
index 0000000000000000000000000000000000000000..23a62d36b8a52046dc49b7701adcfe35c9534f65
Binary files /dev/null and b/tests/TestData/Eq02.npz differ
diff --git a/tests/TestData/Eq04.npz b/tests/TestData/Eq04.npz
new file mode 100644
index 0000000000000000000000000000000000000000..591c07db40d299c6ec7ea34266ee049d84cbaa8b
Binary files /dev/null and b/tests/TestData/Eq04.npz differ
diff --git a/tests/TestData/EqN.npz b/tests/TestData/EqN.npz
new file mode 100644
index 0000000000000000000000000000000000000000..aac5a727d530f047c3f6ddd9b4d3568368d2c64d
Binary files /dev/null and b/tests/TestData/EqN.npz differ
diff --git a/examples/__init__.py b/tests/__init__.py
similarity index 100%
rename from examples/__init__.py
rename to tests/__init__.py
diff --git a/tests/__pycache__/TestEq04.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/TestEq04.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..91bcbb2c55bb8114a2c2beb58e7dc24c3d87b50a
Binary files /dev/null and b/tests/__pycache__/TestEq04.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/__init__.cpython-312.pyc b/tests/__pycache__/__init__.cpython-312.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..089662da0b56caba69f0501d02dac24686b4e462
Binary files /dev/null and b/tests/__pycache__/__init__.cpython-312.pyc differ
diff --git a/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..aebed3727349670ee4178aee4a98d1060dad63ca
Binary files /dev/null and b/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..df145f680fec694c6c644d8d8fdd81d57546a1f8
Binary files /dev/null and b/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..389b1d06261f49dbfb4ffe9d76c7efeda9f2ad47
Binary files /dev/null and b/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..6316e86eb0c6a2b14ab1657d500b84b1b8dfdab3
Binary files /dev/null and b/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..3f329ce080859a5eda5164959c5001f4000d54fa
Binary files /dev/null and b/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/test_ElectrolyticDiode.py b/tests/test_ElectrolyticDiode.py
new file mode 100644
index 0000000000000000000000000000000000000000..b6aedddf45dcf5bfd3ee0d01a285ee9c33421304
--- /dev/null
+++ b/tests/test_ElectrolyticDiode.py
@@ -0,0 +1,94 @@
+''' 
+Tests the ElectrolyticDiode implementation in src.ElectrolyticDiode.py
+'''
+
+
+from src.ElectrolyticDiode import ElectrolyticDiode
+import numpy as np
+
+# Define parameter to use in the test
+phi_bias = 10
+g_phi = 350
+y_fixed = 0.01
+z_A = -1.0
+z_C = 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+number_cells = [20, 128]
+Lx = 0.02
+Ly = 0.1
+x0 = np.array([0, 0])
+x1 = np.array([Lx, Ly])
+refinement_style = 'uniform'
+solvation = 5
+PoissonBoltzmann = False
+rtol = 1e-3 
+relax_param = 0.15 # 0.1
+max_iter = 15_000
+    
+
+data_comparison = np.load('tests/TestData/ElectrolyticDiode.npz')
+
+def test_ForwardBias():
+    y_A_ForwardBias, y_C_ForwardBias, phi_ForwardBias, p_ForwardBias, x_ForwardBias = ElectrolyticDiode('ForwardBias', phi_bias, g_phi, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x_ForwardBias'], x_ForwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'ForwardBias grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_ForwardBias'], y_A_ForwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A_ForwardBias not calculated correctly'
+    assert np.allclose(data_comparison['y_C_ForwardBias'], y_C_ForwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C_ForwardBias not calculated correctly'
+    assert np.allclose(data_comparison['phi_ForwardBias'], phi_ForwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi_ForwardBias not calculated correctly'
+    assert np.allclose(data_comparison['p_ForwardBias'], p_ForwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p_ForwardBias not calculated correctly'
+    
+def test_NoBias():
+    y_A_NoBias, y_C_NoBias, phi_NoBias, p_NoBias, x_NoBias = ElectrolyticDiode('NoBias', phi_bias, g_phi, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x_NoBias'], x_NoBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'NoBias grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_NoBias'], y_A_NoBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A_NoBias not calculated correctly'
+    assert np.allclose(data_comparison['y_C_NoBias'], y_C_NoBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C_NoBias not calculated correctly'
+    assert np.allclose(data_comparison['phi_NoBias'], phi_NoBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi_NoBias not calculated correctly'
+    assert np.allclose(data_comparison['p_NoBias'], p_NoBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p_NoBias not calculated correctly'
+    
+def test_BackwardBias():
+    y_A_BackwardBias, y_C_BackwardBias, phi_BackwardBias, p_BackwardBias, x_BackwardBias = ElectrolyticDiode('BackwardBias', phi_bias, g_phi, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x_BackwardBias'], x_BackwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'BackwardBias grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_BackwardBias'], y_A_BackwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A_BackwardBias not calculated correctly'
+    assert np.allclose(data_comparison['y_C_BackwardBias'], y_C_BackwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C_BackwardBias not calculated correctly'
+    assert np.allclose(data_comparison['phi_BackwardBias'], phi_BackwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi_BackwardBias not calculated correctly'
+    assert np.allclose(data_comparison['p_BackwardBias'], p_BackwardBias,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p_BackwardBias not calculated correctly'
\ No newline at end of file
diff --git a/tests/test_Eq02.py b/tests/test_Eq02.py
new file mode 100644
index 0000000000000000000000000000000000000000..6356f24d54da28be63be14f3f5af6ab7e70c69f2
--- /dev/null
+++ b/tests/test_Eq02.py
@@ -0,0 +1,94 @@
+''' 
+Tests the Eq02 implementation in src.Eq02.py
+'''
+
+
+from src.Eq02 import solve_System_2eq
+import numpy as np
+
+# Define parameter to use in the test
+phi_left = 10.0
+phi_right = 0.0
+p_right = 0.0
+y_A_R = 0.01
+y_C_R = 0.01
+z_A = -1.0
+z_C = 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+solvation = 0
+number_cells = 128
+relax_param = .1
+rtol = 1e-4
+max_iter = 500
+
+data_comparison = np.load('tests/TestData/Eq02.npz')
+
+def test_standard():
+    y_A, y_C, phi, p, x = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x'], x,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Uniform grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A'], y_A,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A not calculated correctly'
+    assert np.allclose(data_comparison['y_C'], y_C,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C not calculated correctly'
+    assert np.allclose(data_comparison['phi'], phi,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi not calculated correctly'
+    assert np.allclose(data_comparison['p'], p,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p not calculated correctly'
+
+
+def test_solvation():
+    solvation = 5
+    y_A_solvation, y_C_solvation, phi_solvation, p_solvation, x_solvation = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='log', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with logarithmic refinement
+    assert np.allclose(data_comparison['x_solvation'], x_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'log grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_solvation'], y_A_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A with solvation not calculated correctly'
+    assert np.allclose(data_comparison['y_C_solvation'], y_C_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C with solvation not calculated correctly'
+    assert np.allclose(data_comparison['phi_solvation'], phi_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi with solvation not calculated correctly'
+    assert np.allclose(data_comparison['p_solvation'], p_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p with solvation not calculated correctly'
+
+def test_compressibility():
+    solvation = 0
+    K = 5_000
+    y_A_compressible, y_C_compressible, phi_compressible, p_compressible, x_compressible = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='hard_log', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with hard-logarithmic refinement
+    assert np.allclose(data_comparison['x_compressible'], x_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'hard-log grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_compressible'], y_A_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A with compressible not calculated correctly'
+    assert np.allclose(data_comparison['y_C_compressible'], y_C_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C with compressible not calculated correctly'
+    assert np.allclose(data_comparison['phi_compressible'], phi_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi with compressible not calculated correctly'
+    assert np.allclose(data_comparison['p_compressible'], p_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p with compressible not calculated correctly'
+    
\ No newline at end of file
diff --git a/tests/test_Eq04.py b/tests/test_Eq04.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d2652e8c8313f206af14c0540b47ca5c469f1b4
--- /dev/null
+++ b/tests/test_Eq04.py
@@ -0,0 +1,93 @@
+''' 
+Tests the Eq04 implementation in src.Eq04.py
+'''
+
+from src.Eq04 import solve_System_4eq
+import numpy as np
+
+# Define parameter to use in the test
+phi_left = 10.0
+phi_right = 0.0
+p_right = 0.0
+y_A_R = 0.01
+y_C_R = 0.01
+z_A = -1.0
+z_C = 1.0
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+solvation = 0
+number_cells = 128
+relax_param = .1
+rtol = 1e-4
+max_iter = 500
+
+data_comparison = np.load('tests/TestData/Eq04.npz')
+
+def test_standard():
+    y_A, y_C, phi, p, x = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x'], x,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Uniform grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A'], y_A,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A not calculated correctly'
+    assert np.allclose(data_comparison['y_C'], y_C,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C not calculated correctly'
+    assert np.allclose(data_comparison['phi'], phi,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi not calculated correctly'
+    assert np.allclose(data_comparison['p'], p,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p not calculated correctly'
+
+
+def test_solvation():
+    solvation = 5
+    y_A_solvation, y_C_solvation, phi_solvation, p_solvation, x_solvation = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='log', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with logarithmic refinement
+    assert np.allclose(data_comparison['x_solvation'], x_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'log grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_solvation'], y_A_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A with solvation not calculated correctly'
+    assert np.allclose(data_comparison['y_C_solvation'], y_C_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C with solvation not calculated correctly'
+    assert np.allclose(data_comparison['phi_solvation'], phi_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi with solvation not calculated correctly'
+    assert np.allclose(data_comparison['p_solvation'], p_solvation,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p with solvation not calculated correctly'
+
+def test_compressibility():
+    solvation = 0
+    K = 5_000
+    y_A_compressible, y_C_compressible, phi_compressible, p_compressible, x_compressible = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='hard_log', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with hard-logarithmic refinement
+    assert np.allclose(data_comparison['x_compressible'], x_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'hard-log grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_A_compressible'], y_A_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_A with compressible not calculated correctly'
+    assert np.allclose(data_comparison['y_C_compressible'], y_C_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_C with compressible not calculated correctly'
+    assert np.allclose(data_comparison['phi_compressible'], phi_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi with compressible not calculated correctly'
+    assert np.allclose(data_comparison['p_compressible'], p_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p with compressible not calculated correctly'
+    
\ No newline at end of file
diff --git a/tests/test_EqN.py b/tests/test_EqN.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc68dfb8aa2a79026540993821494d150ae21dfe
--- /dev/null
+++ b/tests/test_EqN.py
@@ -0,0 +1,66 @@
+''' 
+Tests the EqN implementation in src.EqN.py
+'''
+
+
+from src.EqN import solve_System_Neq
+import numpy as np
+
+# Define parameter to use in the test
+phi_left = 10.0
+phi_right = 0.0
+p_right = 0.0
+y_R = [3/6, 1/6, 1/6]
+z_alpha = [-1.0, 1.0, 2.0]
+K = 'incompressible'
+Lambda2 = 8.553e-6
+a2 = 7.5412e-4
+solvation = 0
+number_cells = 128
+relax_param = .1
+rtol = 1e-4
+max_iter = 500
+
+data_comparison = np.load('tests/TestData/EqN.npz')
+
+def test_1():
+    y_alpha, phi, p, x = solve_System_Neq(phi_left, phi_right, p_right, z_alpha, y_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x'], x,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Uniform grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_alpha'], y_alpha,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_alpha not calculated correctly'
+    assert np.allclose(data_comparison['phi'], phi,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi not calculated correctly'
+    assert np.allclose(data_comparison['p'], p,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p not calculated correctly'
+    
+def test_2():
+    y_R = [1/23, 1/23, 1/23, 1/23, 10/23]
+    z_alpha = [-4.0, -3.0, -2.0, -1.0, 1.0]
+    relax_param = 0.01
+    max_iter = 10_000
+
+    y_alpha_2, phi_2, p_2, x_2 = solve_System_Neq(phi_left, phi_right, p_right, z_alpha, y_R, K, Lambda2, a2, number_cells, solvation=solvation,  relax_param=relax_param, x0=0, x1=1, refinement_style='log', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+    # test grid generation with no refinement
+    assert np.allclose(data_comparison['x_2'], x_2,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Log grid not generated correctly'
+    # test solution vector
+    assert np.allclose(data_comparison['y_alpha_2'], y_alpha_2,\
+                        rtol=1e-15, atol=1e-15),\
+                        'y_alpha_2 not calculated correctly'
+    assert np.allclose(data_comparison['phi_2'], phi_2,\
+                        rtol=1e-15, atol=1e-15),\
+                        'phi_2 not calculated correctly'
+    assert np.allclose(data_comparison['p_2'], p_2,\
+                        rtol=1e-15, atol=1e-15),\
+                        'p_2 not calculated correctly'
+    
\ No newline at end of file
diff --git a/tests/test_Helper_DoubleLayerCapacity.py b/tests/test_Helper_DoubleLayerCapacity.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f16e3ee275608c36cd1a18266d9538d15d8055c
--- /dev/null
+++ b/tests/test_Helper_DoubleLayerCapacity.py
@@ -0,0 +1,195 @@
+''' 
+Tests the Helpfer functions implemented in src.Helper_DoubleLayerCapacity.py
+'''
+
+
+from src.Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_num_dim, Q_DL_dimless_ana, Q_DL_dim_ana, C_DL_dimless_ana, C_DL_dim_ana
+from src.Eq04 import solve_System_4eq
+from src.Eq02 import solve_System_2eq
+import numpy as np
+
+# 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-9
+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_incompresible = 'incompressible'
+K_compressible = 20_000
+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 = 128
+relax_param = 0.03
+p_right = 0
+rtol = 1e-4 # ! Change back to 1e-8
+
+
+# phi^L domain
+Vol_start = 0.05 # ! Change back to 0
+Volt_end = 0.75
+n_Volts = 10#0
+
+phi_left_dim = np.linspace(Vol_start, Volt_end, n_Volts)
+phi_left = phi_left_dim * e0/(k*T)
+
+data_comparison = np.load('tests/TestData/DLKap.npz')
+
+
+
+# Incompressible tests
+def test_incompressible():
+    # Calculate the numerical solution
+    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_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K_incompresible, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=rtol, 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_)
+
+    # Calculate the numerical charge and capacity with finite differences and numerical integration
+    Q_num = []
+    Q_num_dim_ = []
+    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_incompresible), x_num[j]))
+        Q_num_dim_.append(Q_num_dim(y_A_num[j], y_C_num[j], n(p_num[j], K_incompresible), x_num[j], z_A, z_C, nR_m, e0, LR))
+    Q_num = np.array(Q_num)
+    Q_num_dim_ = np.array(Q_num_dim_)
+
+    dx_ = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
+    C_dl_num = C_dl(Q_num, phi_left)
+
+    dx_dim = phi_left_dim[1] - phi_left_dim[0] 
+    C_dl_num_dim_ = C_dl(Q_num_dim_, phi_left_dim)
+
+
+    # Calculate the analytical charge and capacity
+    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_incompresible, Lambda2, a2, kappa)
+    Q_ana_dim_ = Q_DL_dim_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, K_incompresible, Lambda2, a2, nR_m, e0, LR, kappa)
+    C_dl_ana = C_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, Phi_pot_center(phi_left), phi_right, p_right, K_incompresible, Lambda2, a2, kappa)
+    C_dl_ana_dim = C_DL_dim_ana(y_R, y_R, 1-2*y_R, z_A, z_C, Phi_pot_center(phi_left), phi_right, p_right, K_incompresible, Lambda2, a2, nR_m, e0, LR, k, T, kappa)
+
+    # Test the numerical charge and capacity
+    assert np.allclose(data_comparison['Q_num'], Q_num,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Numerical charge not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['Q_num_dim_'], Q_num_dim_,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Numerical charge not calculated correctly for dimensionless'
+    assert np.allclose(data_comparison['C_dl_num'], C_dl_num,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Numerical capacity not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['C_dl_num_dim_'], C_dl_num_dim_,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Numerical capacity not calculated correctly for dimensionless'
+    
+    # Test the analytical charge and capacity
+    assert np.allclose(data_comparison['Q_ana'], Q_ana,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Analytical charge not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['Q_ana_dim_'], Q_ana_dim_,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Analytical charge not calculated correctly for dimensionless'
+    assert np.allclose(data_comparison['C_dl_ana'], C_dl_ana,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Analytical capacity not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['C_dl_ana_dim'], C_dl_ana_dim,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Incompressible: Analytical capacity not calculated correctly for dimensionless'
+    
+# Compressible tests
+def test_compressible():
+    # Calculate the numerical solution
+    y_A_num_compressible, y_C_num_compressible, y_S_num_compressible, phi_num_compressible, p_num_compressible, x_num_compressible = [], [], [], [], [], []
+    for i, phi_bcs in enumerate(phi_left):
+        y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K_compressible, 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_compressible.append(y_A_)
+        y_C_num_compressible.append(y_C_)
+        y_S_num_compressible.append(y_S_)
+        phi_num_compressible.append(phi_)
+        p_num_compressible.append(p_)
+        x_num_compressible.append(x_)
+
+    # Calculate the numerical charge and capacity with finite differences and numerical integration
+    Q_num_compressible = []
+    Q_num_dim_compressible = []
+    for j in range(len(phi_left)):
+        Q_num_compressible.append(Q_num_(y_A_num_compressible[j], y_C_num_compressible[j], n(p_num_compressible[j], K_compressible), x_num_compressible[j]))
+        Q_num_dim_compressible.append(Q_num_dim(y_A_num_compressible[j], y_C_num_compressible[j], n(p_num_compressible[j], K_compressible), x_num_compressible[j], z_A, z_C, nR_m, e0, LR))
+    Q_num_compressible = np.array(Q_num_compressible)
+    Q_num_dim_compressible = np.array(Q_num_dim_compressible)
+
+    dx__compressible = phi_left[1] - phi_left[0] # [1/V], Assumption: phi^L is uniformly distributed
+    C_dl_num_compressible = C_dl(Q_num_compressible, phi_left)
+
+    dx_dim_compressible = phi_left_dim[1] - phi_left_dim[0] 
+    C_dl_num_dim_compressible = C_dl(Q_num_dim_compressible, phi_left_dim)
+
+    # Calculate the analytical charge and capacity
+    Q_ana_compressible = Q_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, K_compressible, Lambda2, a2, kappa)
+    Q_ana_dim__compressible = Q_DL_dim_ana(y_R, y_R, 1-2*y_R, z_A, z_C, phi_left, phi_right, p_right, K_compressible, Lambda2, a2, nR_m, e0, LR, kappa)
+    C_dl_ana_compressible = C_DL_dimless_ana(y_R, y_R, 1-2*y_R, z_A, z_C, Phi_pot_center(phi_left), phi_right, p_right, K_compressible, Lambda2, a2, kappa)
+    C_dl_ana_dim_compressible = C_DL_dim_ana(y_R, y_R, 1-2*y_R, z_A, z_C, Phi_pot_center(phi_left), phi_right, p_right, K_compressible, Lambda2, a2, nR_m, e0, LR, k, T, kappa)
+
+    # Test the numerical charge and capacity
+    assert np.allclose(data_comparison['Q_num_compressible'], Q_num_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Numerical charge not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['Q_num_dim_compressible'], Q_num_dim_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Numerical charge not calculated correctly for dimensionless'
+    assert np.allclose(data_comparison['C_dl_num_compressible'], C_dl_num_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Numerical capacity not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['C_dl_num_dim_compressible'], C_dl_num_dim_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Numerical capacity not calculated correctly for dimensionless'
+    
+    # Test the analytical charge and capacity
+    assert np.allclose(data_comparison['Q_ana_compressible'], Q_ana_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Analytical charge not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['Q_ana_dim__compressible'], Q_ana_dim__compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Analytical charge not calculated correctly for dimensionless'
+    assert np.allclose(data_comparison['C_dl_ana_compressible'], C_dl_ana_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Analytical capacity not calculated correctly for dimensions'
+    assert np.allclose(data_comparison['C_dl_ana_dim_compressible'], C_dl_ana_dim_compressible,\
+                        rtol=1e-15, atol=1e-15),\
+                        'Compressible: Analytical capacity not calculated correctly for dimensionless'
+    
+# Test the functions
+def test_Phi_pot_center():
+    assert np.allclose(Phi_pot_center(phi_left), data_comparison['Phi_pot_center'],\
+                        rtol=1e-15, atol=1e-15),\
+                        'Phi_pot_center not calculated correctly'
+    
+def test_dx():
+    assert np.allclose(dx(phi_left), data_comparison['dx'],\
+                        rtol=1e-15, atol=1e-15),\
+                        'dx not calculated correctly'
+def test_n():
+    assert np.allclose(n(data_comparison['p_num_compressible'][0], K_compressible), data_comparison['n'], rtol=1e-15, atol=1e-15),\
+                        'n not calculated correctly'
+    
\ No newline at end of file