diff --git a/Project.toml b/Project.toml
index 24916a6..21b6fd5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,18 +5,21 @@ version = "0.1.0"
 
 [deps]
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [compat]
 DomainSets = "0.5, 0.6"
-Interpolations = "0.12"
+Interpolations = "0.14"
 Markdown = "1"
-ModelingToolkit = "8"
+ModelingToolkit = "8.48"
 OrdinaryDiffEq = "6"
+SciMLBase = "1.88"
 
 [extras]
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
diff --git a/README.md b/README.md
index f860a3f..8e97b37 100644
--- a/README.md
+++ b/README.md
@@ -6,7 +6,7 @@ These can then be solved with the help of the various discretizer packages of Sc
 - [MethodOfLines.jl](https://www.github.com/SciML/MethodOfLines.jl)
 - [NeuralPDE.jl](https://www.github.com/SciML/NeuralPDE.jl)
 
-It can be used for benchmarking, verification, research in to discretization methods, and any other ideas you might have.
+It can be used for benchmarking, verification, showing off your work, research in to discretization methods, and any other ideas you might have.
 
 If you have a well posed system, please add it! Any and all PDE systems are welcome, even if they cannot currently be solved by discretizer packages.
 Please include a short abstract where possible, explaining where the system arises to aid future readers and large language models.
diff --git a/lib/brusselator.jl b/lib/brusselator.jl
index 7a15cc4..919b5a0 100644
--- a/lib/brusselator.jl
+++ b/lib/brusselator.jl
@@ -138,7 +138,7 @@ function brusselator_2d()
     tags = ["2D", "Brusselator", "Periodic", "Nonlinear", "Reaction", "Diffusion"]
 
     @named bruss = PDESystem(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)],
-                             analytic_func = analytic_func)
+                             analytic_func = analytic_func, metadata = tags)
 
     bruss
 end
diff --git a/lib/linear_convection.jl b/lib/linear_convection.jl
new file mode 100644
index 0000000..9051a9e
--- /dev/null
+++ b/lib/linear_convection.jl
@@ -0,0 +1,465 @@
+# Linear Convection Systems
+"""
+# Linear Convection Equation
+
+1D linear convection equation with periodic boundary conditions.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+It can take any initial condition profile, as long as it is periodic.
+
+```math
+\\frac{\\partial u}{\\partial t} + c \\frac{\\partial u}{\\partial x} = 0
+```
+"""
+function linear_convection(f, ps, name = :linear_convection)
+    @variables u(..)
+    @parameters t x v
+    Dt = Differential(t)
+    Dx = Differential(x)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ -v * Dx(u(t, x))
+
+    bcs = [u(0, x) ~ f(x),
+        u(t, 0) ~ u(t, 1)]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 10.0),
+        x ∈ Interval(0.0, 1.0)]
+
+    # Analytic solution
+    u_exact = [u(t, x) ~ f(x - v * t)]
+
+    tags = ["1D", "Periodic", "Linear", "Advection"]
+
+    # PDE system
+    lin_conv = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => ps[1]],
+                         analytic = u_exact, metadata = tags, name = name)
+
+    lin_conv
+end
+
+# sinusoidal input
+convsin = linear_convection(x -> sin(2pi * x), [1], :convsin)
+push!(convsin.metadata, "Sinusoidal")
+push!(all_systems, convsin)
+
+#cosine input
+convcos = linear_convection(x -> cos(2pi * x), [-2], :convcos)
+push!(convcos.metadata, "Sinusoidal")
+push!(all_systems, convcos)
+
+# triangular input
+convtri = linear_convection(x -> 1.0 - abs(x - floor(x + 0.5)), [0.6], :convtri)
+push!(convtri.metadata, "Triangular")
+push!(all_systems, convtri)
+
+# square wave
+f = (x) -> IfElse.ifelse(x - floor(x) < 0.5, 1.0, -1.0)
+
+convsquare = linear_convection(f, [-1.1], :convsquare)
+push!(convsquare.metadata, "Square")
+push!(convsquare.metadata, "Discontinuous")
+push!(all_systems, convsquare)
+
+convsquare = linear_convection(f, [2.1], :convsquare2)
+push!(convsquare.metadata, "Square")
+push!(convsquare.metadata, "Discontinuous")
+push!(all_systems, convsquare)
+
+"""
+# Linear Convection Equation with Dirichlet Boundary Conditions 1
+
+1D linear convection equation with Dirichlet boundary conditions.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+v must be positive.
+
+It can take any initial condition profile, as long as it is windowed over the domain.
+
+```math
+\\frac{\\partial u}{\\partial t} + v \\frac{\\partial u}{\\partial x} = 0
+```
+
+"""
+function linear_convection_dirichlet1(f, ps, name = :linear_convection)
+    @variables u(..)
+    @parameters t x v
+    Dt = Differential(t)
+    Dx = Differential(x)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ -v * Dx(u(t, x))
+
+    bcs = [u(0, x) ~ f(x),
+        u(t, 0) ~ 0.0]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 2 / ps[1]),
+        x ∈ Interval(0.0, 1.0)]
+
+    # Analytic solution
+    u_exact = [u(t, x) ~ windowlower(f(x - v * t), x - v * t)]
+
+    tags = ["1D", "Dirichlet", "Linear", "Advection"]
+
+    # PDE system
+    lin_conv = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => ps[1]],
+                         analytic = u_exact, metadata = tags, name = name)
+
+    lin_conv
+end
+
+function windowlower(f, x, domain = (0.0, 1.0))
+    IfElse.ifelse(x <= domain[1], IfElse.ifelse(x > domain[2], f, 0.0), 0.0)
+end
+
+# sinusoidal input
+convsin = linear_convection_dirichlet1(x -> sin(2pi * x), [1], :dconvsin)
+push!(convsin.metadata, "Sinusoidal")
+push!(all_systems, convsin)
+
+convsin = linear_convection_dirichlet1(x -> sin(2pi * x), [0.5], :dconvsin2)
+push!(convsin.metadata, "Sinusoidal")
+push!(all_systems, convsin)
+
+#cosine input
+convcos = linear_convection_dirichlet1(x -> cos(2pi * x), [2], :dconvcos)
+push!(convcos.metadata, "Sinusoidal")
+push!(all_systems, convcos)
+
+# triangular input
+convtri = linear_convection_dirichlet1(x -> 1.0 - abs(x - floor(x + 0.5)), [0.6], :dconvtri)
+push!(convtri.metadata, "Triangular")
+push!(all_systems, convtri)
+
+# square wave
+f = (x) -> IfElse.ifelse(x - floor(x) < 0.5, 1.0, -1.0)
+
+convsquare = linear_convection_dirichlet1(f, [1.1], :dconvsquare)
+push!(convsquare.metadata, "Square")
+push!(convsquare.metadata, "Discontinuous")
+push!(all_systems, convsquare)
+
+"""
+# Linear Convection Equation with Dirichlet Boundary Conditions 2
+
+1D linear convection equation with Dirichlet boundary conditions.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+v must be positive.
+
+It can take any initial condition profile, as long as it is windowed over the domain.
+
+```math
+\\frac{\\partial u}{\\partial t} - v \\frac{\\partial u}{\\partial x} = 0
+```
+
+"""
+function linear_convection_dirichlet2(f, ps, name = :linear_convection)
+    @variables u(..)
+    @parameters t x v
+    Dt = Differential(t)
+    Dx = Differential(x)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ v * Dx(u(t, x))
+
+    bcs = [u(0, x) ~ f(x),
+        u(t, 2 / ps[1]) ~ 0.0]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 2 / ps[1]),
+        x ∈ Interval(0.0, 1.0)]
+
+    # Analytic solution
+    u_exact = [u(t, x) ~ windowlower(f(x - v * t), x - v * t)]
+
+    tags = ["1D", "Dirichlet", "Linear", "Advection"]
+
+    # PDE system
+    lin_conv = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => ps[1]],
+                         analytic = u_exact, metadata = tags, name = name)
+
+    lin_conv
+end
+
+function windowupper(f, x, domain = (0.0, 1.0))
+    IfElse.ifelse(x < domain[1], IfElse.ifelse(x >= domain[2], f, 0.0), 0.0)
+end
+
+# sinusoidal input
+convsin = linear_convection_dirichlet1(x -> sin(2pi * x), [1], :ddconvsin)
+push!(convsin.metadata, "Sinusoidal")
+push!(all_systems, convsin)
+
+convsin = linear_convection_dirichlet1(x -> sin(2pi * x), [0.5], :ddconvsin2)
+push!(convsin.metadata, "Sinusoidal")
+push!(all_systems, convsin)
+
+#cosine input
+convcos = linear_convection_dirichlet1(x -> cos(2pi * x), [0.7], :ddconvcos)
+push!(convcos.metadata, "Sinusoidal")
+push!(all_systems, convcos)
+
+# triangular input
+convtri = linear_convection_dirichlet1(x -> 1.0 - abs(x - floor(x + 0.5)), [0.1],
+                                       :ddconvtri)
+push!(convtri.metadata, "Triangular")
+push!(all_systems, convtri)
+
+# square wave
+sq = f = (x) -> IfElse.ifelse(x - floor(x) < 0.5, 1.0, -1.0)
+
+convsquare = linear_convection_dirichlet1(f, [3.1], :ddconvsquare)
+push!(convsquare.metadata, "Square")
+push!(convsquare.metadata, "Discontinuous")
+push!(all_systems, convsquare)
+
+"""
+# Linear Convection Equation with Dirichlet Boundary Conditions 3
+
+1D linear convection equation with Dirichlet boundary conditions. the left boundary is time dependent.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+v must be positive.
+
+It can take any initial condition profile, as long as it is windowed over the domain.
+
+```math
+\\frac{\\partial u}{\\partial t} + v \\frac{\\partial u}{\\partial x} = 0
+```
+
+"""
+function linear_convection_dirichlet3(f, h, ps, name = :linear_convection)
+    @variables u(..)
+    @parameters t x v
+    Dt = Differential(t)
+    Dx = Differential(x)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ -v * Dx(u(t, x))
+
+    bcs = [u(0, x) ~ f(x),
+        u(t, 0.0) ~ h(t)]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 10.0),
+        x ∈ Interval(0.0, 1.0)]
+
+    # Analytic solution
+    u_exact = [u(t, x) ~ IfElse.ifelse(x > v * t, f(x - v * t), h(t - x / v))]
+
+    tags = ["1D", "Dirichlet", "Linear", "Advection"]
+
+    # PDE system
+    lin_conv = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => ps[1]],
+                         analytic = u_exact, metadata = tags, name = name)
+
+    lin_conv
+end
+
+funcs = [x -> x, x -> x^2, x -> x^3, sinpi, cospi, x -> 1.0 - abs(x - floor(x + 0.5)), sq]
+function add_systems!(all_systems, funcs, sysconstructor, name)
+    i = 0
+    for f in funcs
+        for h in funcs
+            conv = sysconstructor(f, h, [rand()], Symbol(name, i))
+            push!(all_systems, conv)
+            i += 1
+            conv = sysconstructor(f, x -> -2 * h(x), [2 * rand()], Symbol(name, i))
+            push!(all_systems, conv)
+            i += 1
+            conv = sysconstructor(x -> -10 * f(x), h, [rand()], Symbol(name, i))
+            push!(all_systems, conv)
+            i += 1
+            conv = sysconstructor(x -> -6 * f(x), x -> -5 * h(x), [rand()],
+                                  Symbol(name, i))
+            push!(all_systems, conv)
+            i += 1
+        end
+    end
+end
+
+add_systems!(all_systems, funcs, linear_convection_dirichlet3, "funcconv")
+
+"""
+# Linear Convection Equation with Dirichlet Boundary Conditions 4
+
+1D linear convection equation with Dirichlet boundary conditions. the left boundary is time dependent.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+v must be positive.
+
+It can take any initial condition profile, as long as it is windowed over the domain.
+
+```math
+\\frac{\\partial u}{\\partial t} - v \\frac{\\partial u}{\\partial x} = 0
+```
+
+"""
+function linear_convection_dirichlet4(f, h, ps, name = :linear_convection)
+    @variables u(..)
+    @parameters t x v
+    Dt = Differential(t)
+    Dx = Differential(x)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ v * Dx(u(t, x))
+
+    bcs = [u(0, x) ~ f(x),
+        u(t, 1.0) ~ h(t)]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 10.0),
+        x ∈ Interval(0.0, 1.0)]
+
+    # Analytic solution
+    u_exact = [u(t, x) ~ IfElse.ifelse(x < v * t, f(x + v * t), h(t + x / v))]
+
+    tags = ["1D", "Dirichlet", "Linear", "Advection"]
+
+    # PDE system
+    lin_conv = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => ps[1]],
+                         analytic = u_exact, metadata = tags, name = name)
+
+    lin_conv
+end
+
+add_systems!(all_systems, funcs, linear_convection_dirichlet4, "funcconvneg")
+
+"""
+# Convection Diffusion Equation in 1D
+
+1D convection diffusion equation with Dirichlet boundary conditions.
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity and diffusing in a medium.
+
+It is initialized with a zero profile, and the boundary conditions are set to 1 at the right boundary and 0 at the left boundary.
+
+```math
+\\frac{\\partial u}{\\partial t} = \\kappa \\frac{\\partial^2 u}{\\partial z^2} + v \\frac{\\partial u}{\\partial z} +
+```
+
+Reference Solution found here: https://www.12000.org/my_notes/diffusion_convection_PDE/diffusion_convection_PDE.htm
+Author of reference solution: Nasser M. Abbasi
+
+"""
+function convection_diffusion(L, ps, name = :convection_diffusion)
+    @variables f(..)
+    @parameters z, t, k, v
+    Dt = Differential(t)
+    Dz = Differential(z)
+    Dzz = Differential(z)^2
+
+    # 1D PDE and boundary conditions
+    eq = Dt(f(t, z)) ~ k * Dzz(f(t, z)) + v * Dz(f(t, z))
+
+    f_0(z) = IfElse.ifelse(z == L, 1.0, 0.0)
+
+    bcs = [f(0, z) ~ f_0(z),
+        f(t, 0) ~ 0.0, f(t, L) ~ 1.0]
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 30.0),
+        z ∈ Interval(0.0, L)]
+
+    # Analytic/reference solution
+    λ(n) = (n * π / L)^2
+
+    maxiters = 1000
+
+    function A(ps, t, z)
+        k = ps[1]
+        v = ps[2]
+        exp(-((t * v^2) / (4k) + (v * z) / (2k)))
+    end
+
+    function u(ps, t, z)
+        k = ps[1]
+        v = ps[2]
+        a = z / L * exp((t * v^2) / (4k) + (v * L) / (2k))
+        s = mapreduce((+), 1:maxiters) do n
+            acu = (2 * (-1)^n * v^2 * exp(-k * λ(n) * t + (v * L) / (2k)) *
+                   (exp(k * λ(n) * t + (t * v^2) / (4k)) - 1)) /
+                  (n * π * (4 * k^2 * λ(n) + v^2))
+            acu += 2 / L * ((-1)^n) * exp(v * L / (2k)) * exp(-k * λ(n) * t) / sqrt(λ(n))
+            acu *= sin(λ(n) * z)
+
+            acu
+        end
+        a + s
+    end
+
+    ref = [f(t, z) => (ps, t, z) -> A(ps, t, z) * u(ps, t, z)]
+
+    tags = ["1D", "Dirichlet", "Linear", "Advection", "Diffusion", "Monotonic"]
+
+    # PDE system
+    convdiff = PDESystem(eq, bcs, domains, [t, z], [f(t, z)],
+                         [k => ps[1], v => ps[2]], analytic_func = ref, metadata = tags,
+                         name = name)
+
+    convdiff
+end
+
+# L = 1.0, k = 0.01, v = 1.0
+convdiff1 = convection_diffusion(1.0, [0.01, 1.0], :convdiff1)
+push!(all_systems, convdiff1)
+
+# L = 5.0, k = 1.0, v = -1.0
+convdiff2 = convection_diffusion(5.0, [1.0, -1.0], :convdiff2)
+push!(all_systems, convdiff2)
+
+# L = 10.0, k = 0.01, v = -3.0
+convdiff3 = convection_diffusion(10.0, [0.01, -3.0], :convdiff3)
+push!(all_systems, convdiff3)
+
+# L = 10.0, k = 0.5, v = 3.0
+convdiff4 = convection_diffusion(10.0, [0.5, 3.0], :convdiff4)
+push!(all_systems, convdiff4)
+
+"""
+# Transport Equation in 1D
+
+1D transport equation without boundary conditions
+
+This equation models the flow of a substance, where the substance is moving at a constant velocity.
+
+it is initialized with a sinusoid, and has a sinusoidal source term.
+
+```math
+\\frac{\\partial u}{\\partial t} + v 2\\frac{\\partial u}{\\partial z} = \\sin(z)
+```
+"""
+function trans_sin()
+    @variables u(..)
+    @parameters z, t
+    Dt = Differential(t)
+    Dz = Differential(z)
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, z)) + 2 * Dz(u(t, z)) ~ sin(z)
+
+    u_exact(t, z) = sin(z - 2t) + 0.5 * cos(z - 2t) - 0.5 * cos(z)
+
+    bcs = [u(0, z) ~ u_exact(0, z),
+        u(t, 0) ~ u_exact(t, 0),
+        (t, 2π) ~ u_exact(t, 2π)]
+
+    # Space and time domains
+
+    domains = [t ∈ Interval(0.0, 1.0),
+        z ∈ Interval(0.0, 2π)]
+
+    # Analytic/reference solution
+    ref = [u(t, z) ~ u_exact(t, z)]
+
+    tags = ["1D", "Transport", "Linear", "Sinusoidal", "Inhomogeneous", "Advection"]
+
+    # PDESystem
+
+    @named trans_sin = PDESystem(eqs, bcs, domains, [t, z], [u(t, z)], analytic = ref,
+                                 metadata = tags)
+
+    trans_sin
+end
diff --git a/lib/linear_diffusion.jl b/lib/linear_diffusion.jl
index 06a6fd2..5197c40 100644
--- a/lib/linear_diffusion.jl
+++ b/lib/linear_diffusion.jl
@@ -34,5 +34,46 @@ function heat_1d1()
 
     heat_1d1
 end
-
 push!(all_systems, heat_1d1())
+
+"""
+# The Heat Equation in 1D with Neumann Boundary Conditions.
+
+1D heat equation with Neumann boundary conditions.
+This models the temperature of a rod over time, where the ends are held at a constant temperature.
+
+It is initialized with a sinusoidal profile.
+The equation is given by:
+
+```math
+\\frac{\\partial u}{\\partial t} = D \\frac{\\partial^2 u}{\\partial x^2}
+```
+"""
+function heat_1d_neumann()
+    # Method of Manufactured Solutions
+
+    # Parameters, variables, and derivatives
+    @parameters t x
+    @variables u(..)
+    Dt = Differential(t)
+    Dx = Differential(x)
+    Dxx = Differential(x)^2
+
+    # 1D PDE and boundary conditions
+    eq = Dt(u(t, x)) ~ Dxx(u(t, x))
+    bcs = [u(0, x) ~ cos(x),
+        Dx(u(t, 0)) ~ 0,
+        Dx(u(t, Float64(pi))) ~ 0]
+
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 1.0),
+        x ∈ Interval(0.0, Float64(pi))]
+
+    analytic = [u(t, x) ~ exp(-t) * cos(x)]
+
+    tags = ["1D", "Neumann", "Linear", "Diffusion", "Heat"]
+    # PDE system
+    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], analytic = analytic,
+                              metadata = tags)
+end
+push!(all_systems, heat_1d_neumann())
diff --git a/lib/nonlinear_diffusion.jl b/lib/nonlinear_diffusion.jl
new file mode 100644
index 0000000..9d0ae14
--- /dev/null
+++ b/lib/nonlinear_diffusion.jl
@@ -0,0 +1,41 @@
+"""
+# Spherical Laplacian 1D
+
+The laplacian in spherical coordinates, with Dirichlet and Neumann boundary conditions.
+Symmetrical in azimuth and elevation, so the solution is a function of only one spatial
+variable.
+
+```math
+\\frac{\\partial^2 u}{\\partial t^2} - \\frac{4}{r^2} \\frac{\\partial}{\\partial r}(r^2 \\frac{\\partial}{\\partial r}(u) = 0
+```
+"""
+function spherical_laplacian()
+    # Parameters, variables, and derivatives
+    @parameters t r
+    @variables u(..)
+    Dt = Differential(t)
+    Dr = Differential(r)
+
+    # 1D PDE and boundary conditions
+
+    eq = Dt(u(t, r)) ~ 4 / r^2 * Dr(r^2 * Dr(u(t, r)))
+    bcs = [u(0, r) ~ sin(r) / r,
+        Dr(u(t, 0)) ~ 0,
+        u(t, 1) ~ exp(-4t) * sin(1)]
+
+    # Space and time domains
+    domains = [t ∈ Interval(0.0, 1.0),
+        r ∈ Interval(0.0, 1.0)]
+
+    u_exact = [u(t, r) ~ exp.(-4t) * sin.(r) ./ r]
+
+    tags = ["1D", "Dirichlet", "Neumann", "Spherical", "Diffusion"]
+
+    # PDE system
+    @named sph = PDESystem(eq, bcs, domains, [t, r], [u(t, r)], analytic = u_exact,
+                           metadata = tags)
+
+    sph
+end
+
+push!(all_systems, spherical_laplacian())
diff --git a/src/PDESystemLibrary.jl b/src/PDESystemLibrary.jl
index f6f4bd3..06bd07f 100644
--- a/src/PDESystemLibrary.jl
+++ b/src/PDESystemLibrary.jl
@@ -3,6 +3,10 @@ using ModelingToolkit, DomainSets
 using OrdinaryDiffEq
 using Interpolations
 
+import SciMLBase
+
+using IfElse
+using IfElse: ifelse
 using Markdown
 using Random
 
@@ -12,12 +16,15 @@ all_systems = []
 
 include("../lib/burgers.jl")
 include("../lib/linear_diffusion.jl")
+include("../lib/linear_convection.jl")
+include("../lib/nonlinear_diffusion.jl")
 include("../lib/general_linear_system.jl")
 include("../lib/brusselator.jl")
 
-function get_pdesys_with_tags(tags...)
+function get_pdesys_with_tags(withtags; without = [], f = all)
     filter(all_systems) do ex
-        all(t -> t in ex.metadata, tags)
+        b = f(t -> t ∈ ex.metadata, withtags)
+        b && all(t -> t ∉ ex.metadata, withouttags)
     end
 end
 
diff --git a/test/mol_test.jl b/test/mol_test.jl
index 798e553..803f8b2 100644
--- a/test/mol_test.jl
+++ b/test/mol_test.jl
@@ -1,7 +1,7 @@
 using PDESystemLibrary
 PSL = PDESystemLibrary
 
-using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, NonlinearSolve
+using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, NonlinearSolve, Test
 
 N = 100
 
@@ -19,11 +19,13 @@ for ex in PSL.all_systems
             disc = MOLFiniteDifference(dxs)
             prob = discretize(ex, disc)
             sol = NonlinerSolve.solve(prob, NewtonRaphsom())
+            @test sol.retcode == SciMLBase.ReturnCode.Success
         else
             @parameters t
             disc = MOLFiniteDifference(dxs, t)
             prob = discretize(ex, disc)
             sol = solve(prob, FBDF())
+            @test sol.retcode == SciMLBase.ReturnCode.Success
         end
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1c2b39b..2fa4b48 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,8 +9,8 @@ if GROUP == "All" || GROUP == "MOL"
     @time @safetestset "Test against MethodOfLines.jl" begin include("mol_test.jl") end
 end
 
-# TODO: fix this when NeuralPDE.jl can be added to the test environment.
+# TODO: fix this when NeuralPDE.jl can be added to the test environment. (compat with Symbolics.jl 5)
 
-if GROUP == "All" || GROUP == "NeuralPDE"
-    @time @safetestset "Test against NeuralPDE.jl" begin include("neuralpde_test.jl") end
-end
+# if GROUP == "All" || GROUP == "NeuralPDE"
+#     @time @safetestset "Test against NeuralPDE.jl" begin include("neuralpde_test.jl") end
+# end