diff --git a/Project.toml b/Project.toml
index fe81d8b..39528f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PoissonBoltzmannIPAM2025"
 uuid = "480e0da0-c3fa-11f0-becb-6b90fb73a51e"
-version = "0.1.0"
 authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>"]
+version = "0.1.0"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -13,6 +13,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuliaMPBSolver = "d8b18f01-5396-498d-b34d-247825c18ff0"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LessUnitful = "f29f6376-6e90-4d80-80c9-fb8ec61203d5"
@@ -42,6 +43,7 @@ ForwardDiff = "1"
 HypertextLiteral = "0.9"
 InteractiveUtils = "1.11.0"
 Interpolations = "0.16"
+JLD2 = "0.6.3"
 LaTeXStrings = "1.0.1"
 LessUnitful = "1"
 Markdown = "1.11.0"
diff --git a/draftnotebooks/demo-notebook.jl b/draftnotebooks/demo-notebook.jl
deleted file mode 100644
index f617898..0000000
--- a/draftnotebooks/demo-notebook.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-### A Pluto.jl notebook ###
-# v0.20.19
-
-using Markdown
-using InteractiveUtils
-
-# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
-begin
-    using Pkg
-    Pkg.activate(joinpath(@__DIR__, ".."))
-    using PlutoUI
-    using JuliaMPBSolver
-end
-
-# ╔═╡ 5b7da564-8c8e-45c9-8c72-534c5a13301d
-md"""
-# Demo notebook for JuliaMPBSolver
-"""
-
-# ╔═╡ 882dda23-63b9-4b1e-a04e-69071deff69a
-md"""
-This notebook is only relocatable together with the whole {PKGNAME} project.
-All packages used by this notebook need to be added to the environment in 
-$(dirname(joinpath(@__DIR__,".."))).
-"""
-
-# ╔═╡ a8e37976-5db2-485f-87aa-0cf7155e8e00
-JuliaMPBSolver.greet()
-
-# ╔═╡ Cell order:
-# ╟─5b7da564-8c8e-45c9-8c72-534c5a13301d
-# ╟─882dda23-63b9-4b1e-a04e-69071deff69a
-# ╠═60941eaa-1aea-11eb-1277-97b991548781
-# ╠═a8e37976-5db2-485f-87aa-0cf7155e8e00
diff --git a/scripts/simplecell-bsk.jl b/scripts/simplecell-bsk.jl
index 35cacf9..a3a74ea 100644
--- a/scripts/simplecell-bsk.jl
+++ b/scripts/simplecell-bsk.jl
@@ -7,6 +7,7 @@ using Interpolations
 using VoronoiFVM
 using PythonPlot
 using JuliaMPBSolver
+using JLD2
 
 const nel = 20.0 * JuliaMPBSolver.Units.el_surface_density # number of electrons/nm^2 at interfaces
 const M_bulk = 1 # (bulk) molarity at center of domain
@@ -51,120 +52,4 @@ solution, X, nv, ε_r =
     computed_parameters,
 )
 
-function bee!(y, ϕ)
-    N = length(user_parameters.charge_numbers)
-    for i in 1:N
-        y[i] =
-            JuliaMPBSolver.Units.thermal_energy(user_parameters.temperature) *
-            log(c_bulk[i] / computed_parameters.bulk_solvent_concentration) /
-            JuliaMPBSolver.Units.F - user_parameters.charge_numbers[i] * ϕ
-    end
-    return nothing
-end
-
-function bee(sol)
-    n = size(sol, 2)
-    N = length(user_parameters.charge_numbers)
-    e = zeros(N, n)
-    y = zeros(N)
-    for i in 1:n
-        bee!(y, sol[1, i])
-        for j in 1:N
-            e[j, i] = y[j]
-        end
-    end
-    return e
-end
-
-function plotsol(sol; size = (600, 400))
-    PythonPlot.clf()
-    fig, ax = pyplot.subplots(2, 1)
-    ax1 = ax[0]
-    ax2 = ax[1]
-    ax1.grid()
-    ax1r = ax1.twinx()
-
-    c = JuliaMPBSolver.Postprocess.compute_concentrations(
-        sol[1, :],
-        user_parameters,
-        computed_parameters,
-    )
-    cm = c[1, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    cp = c[2, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    c0 = -(sum(c, dims = 1) .- c̄)
-    e = bee(sol)
-    ax1.set_title(
-        "ϕ∈$(round.(Float64.(extrema(sol[1, :])), sigdigits = 3)), ε_r ∈$(round.(Float64.(extrema(ε_r)), sigdigits = 3))",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        sol[1, :],
-        color = "green",
-        linewidth = 2,
-        label = "ϕ",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[1, :],
-        color = "blue",
-        linewidth = 2,
-        label = L"ψ^+",
-        linestyle = "dotted",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[2, :],
-        color = "red",
-        linewidth = 2,
-        linestyle = "dotted",
-        label = L"ψ^+",
-    )
-    ax1r.plot(
-        X[1:(end - 1)] / JuliaMPBSolver.Units.nm,
-        ε_r,
-        color = "pink",
-        linewidth = 3,
-        label = L"ε_r",
-    )
-    ax1.set_ylim(-10, 10)
-    ax1.set_xlabel("z/nm")
-    ax1.set_ylabel("ϕ/V")
-    ax1.legend(loc = (0.1, 0.1))
-    ax1r.legend(loc = (0.8, 0.1))
-    ax1r.set_ylim(0, 80)
-
-    ax2.grid()
-    ax2.set_title("M_avg=$(round.((cm, cp), sigdigits = 3))")
-    ax2.set_xlabel("z/nm")
-    ax2.set_ylabel("c/(mol/L)")
-    ax2.set_ylim(0, 60)
-
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[1, :] / (JuliaMPBSolver.Units.M),
-        color = "blue",
-        linewidth = 2,
-        label = L"c^-",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[2, :] / (JuliaMPBSolver.Units.M),
-        color = "red",
-        linewidth = 2,
-        label = L"c^+",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c0[1, :] / (JuliaMPBSolver.Units.M),
-        color = "green",
-        linewidth = 2,
-        label = L"c_{solvent}",
-    )
-    ax2.legend(loc = (0.4, 0.1))
-
-    tight_layout()
-    savefig("simplecell-bsk.jpg", dpi = 300)
-    return PythonPlot.gcf()
-end
-
-plotsol(solution)
+save("simplecell-bsk.jld2", "solution", solution, "X", X, "nv", nv, "ε_r", ε_r)
diff --git a/scripts/simplecell.jl b/scripts/simplecell.jl
index 0ca3552..934a250 100644
--- a/scripts/simplecell.jl
+++ b/scripts/simplecell.jl
@@ -7,6 +7,7 @@ using Interpolations
 using VoronoiFVM
 using PythonPlot
 using JuliaMPBSolver
+using JLD2
 
 const nel = 20.0 * JuliaMPBSolver.Units.el_surface_density # number of electrons/nm^2 at interfaces
 const M_bulk = 1 # (bulk) molarity at center of domain
@@ -52,120 +53,4 @@ solution, X, nv, ε_r =
     computed_parameters,
 )
 
-function bee!(y, ϕ)
-    N = length(user_parameters.charge_numbers)
-    for i in 1:N
-        y[i] =
-            JuliaMPBSolver.Units.thermal_energy(user_parameters.temperature) *
-            log(c_bulk[i] / computed_parameters.bulk_solvent_concentration) /
-            JuliaMPBSolver.Units.F - user_parameters.charge_numbers[i] * ϕ
-    end
-    return nothing
-end
-
-function bee(sol)
-    n = size(sol, 2)
-    N = length(user_parameters.charge_numbers)
-    e = zeros(N, n)
-    y = zeros(N)
-    for i in 1:n
-        bee!(y, sol[1, i])
-        for j in 1:N
-            e[j, i] = y[j]
-        end
-    end
-    return e
-end
-
-function plotsol(sol; size = (600, 400))
-    PythonPlot.clf()
-    fig, ax = pyplot.subplots(2, 1)
-    ax1 = ax[0]
-    ax2 = ax[1]
-    ax1.grid()
-    ax1r = ax1.twinx()
-
-    c = JuliaMPBSolver.Postprocess.compute_concentrations(
-        sol[1, :],
-        user_parameters,
-        computed_parameters,
-    )
-    cm = c[1, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    cp = c[2, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    c0 = -(sum(c, dims = 1) .- c̄)
-    e = bee(sol)
-    ax1.set_title(
-        "ϕ∈$(round.(Float64.(extrema(sol[1, :])), sigdigits = 3)), ε_r ∈$(round.(Float64.(extrema(ε_r)), sigdigits = 3))",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        sol[1, :],
-        color = "green",
-        linewidth = 2,
-        label = "ϕ",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[1, :],
-        color = "blue",
-        linewidth = 2,
-        label = L"ψ^+",
-        linestyle = "dotted",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[2, :],
-        color = "red",
-        linewidth = 2,
-        linestyle = "dotted",
-        label = L"ψ^+",
-    )
-    ax1r.plot(
-        X[1:(end - 1)] / JuliaMPBSolver.Units.nm,
-        ε_r,
-        color = "pink",
-        linewidth = 3,
-        label = L"ε_r",
-    )
-    ax1.set_ylim(-10, 10)
-    ax1.set_xlabel("z/nm")
-    ax1.set_ylabel("ϕ/V")
-    ax1.legend(loc = (0.1, 0.1))
-    ax1r.legend(loc = (0.8, 0.1))
-    ax1r.set_ylim(0, 80)
-
-    ax2.grid()
-    ax2.set_title("M_avg=$(round.((cm, cp), sigdigits = 3))")
-    ax2.set_xlabel("z/nm")
-    ax2.set_ylabel("c/(mol/L)")
-    ax2.set_ylim(0, 60)
-
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[1, :] / (JuliaMPBSolver.Units.M),
-        color = "blue",
-        linewidth = 2,
-        label = L"c^-",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[2, :] / (JuliaMPBSolver.Units.M),
-        color = "red",
-        linewidth = 2,
-        label = L"c^+",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c0[1, :] / (JuliaMPBSolver.Units.M),
-        color = "green",
-        linewidth = 2,
-        label = L"c_{solvent}",
-    )
-    ax2.legend(loc = (0.4, 0.1))
-
-    tight_layout()
-    savefig("simplecell.jpg", dpi = 300)
-    return PythonPlot.gcf()
-end
-
-plotsol(solution)
+save("simplecell.jld2", "solution", solution, "X", X, "nv", nv, "ε_r", ε_r)
diff --git a/simplecell-bsk.jpg b/simplecell-bsk.jpg
deleted file mode 100644
index 70e265a..0000000
Binary files a/simplecell-bsk.jpg and /dev/null differ
diff --git a/simplecell.jpg b/simplecell.jpg
deleted file mode 100644
index d45df33..0000000
Binary files a/simplecell.jpg and /dev/null differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 0a49493..42abb73 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,9 @@
-using Pkg, JuliaMPBSolver, ExampleJuggler, Test, Markdown
+using Pkg
+using JuliaMPBSolver
+using ExampleJuggler
+using Test
+using Markdown
+using JLD2
 
 ExampleJuggler.verbose!(true)
 
@@ -11,8 +16,31 @@ notebooks = [
     "SymmetricCellSurfaceCharge.jl",
 ]
 
+scripts = [
+    "simplecell.jl",
+    "simplecell-bsk.jl",
+]
+
 @testset "Notebooks" begin
     @testscripts(joinpath(@__DIR__, "..", "notebooks"), notebooks)
 end
 
+@testset "Scripts" begin
+    @testscripts(joinpath(@__DIR__, "..", "scripts"), scripts)
+
+    for script in scripts
+        script_name = script[1:(end - 3)]
+        solution, X, nv, ε_r = load(script_name * ".jld2", "solution", "X", "nv", "ε_r")
+        solution_reference, X_reference, nv_reference, ε_r_reference = load(script_name * "-reference.jld2", "solution", "X", "nv", "ε_r")
+
+        @test solution ≈ solution_reference
+        @test X ≈ X_reference
+        @test nv ≈ nv_reference
+        @test ε_r ≈ ε_r_reference
+
+        rm(script_name * ".jld2")
+    end
+
+end
+
 Pkg.test("JuliaMPBSolver")
diff --git a/test/simplecell-bsk-reference.jld2 b/test/simplecell-bsk-reference.jld2
new file mode 100644
index 0000000..543451a
Binary files /dev/null and b/test/simplecell-bsk-reference.jld2 differ
diff --git a/test/simplecell-reference.jld2 b/test/simplecell-reference.jld2
new file mode 100644
index 0000000..452db0b
Binary files /dev/null and b/test/simplecell-reference.jld2 differ