diff --git a/Project.toml b/Project.toml
index a50d1058a..3b8b5abd4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,7 +29,7 @@ Latexify = "0.15"
 ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
 SafeTestsets = "0.0.1"
-SciMLBase = "1"
+SciMLBase = "1.90, 1.91"
 StaticArrays = "1"
 SymbolicUtils = "1"
 Symbolics = "5"
diff --git a/docs/src/advection_schemes.md b/docs/src/advection_schemes.md
index 97704e673..264cc965e 100644
--- a/docs/src/advection_schemes.md
+++ b/docs/src/advection_schemes.md
@@ -5,7 +5,7 @@ Used as a keyword argument `advection_scheme` to `MOLFiniteDifference`.
 ```julia
 UpwindScheme(approx_order = 1)
 ```
-Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` keyword argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme).
+Changes the direction of the stencil based on the sign of the coefficient of the odd order derivative to be discretized. Scheme order can be increased by changing the `approx_order` argument. For more information, see [Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme).
 
 ## WENO Scheme of Jiang and Shu
 ```julia
@@ -21,4 +21,6 @@ Problems with first order derivatives which multiply one another will need to us
 
 Supports only first order derivatives, other odd order derivatives are unsupported with this scheme. At present, does not support Nonuniform grids, though this is a planned feature.
 
-Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/)
\ No newline at end of file
+Specified on pages 8-9 of [this document](https://repository.library.brown.edu/studio/item/bdr:297524/PDF/)
+
+## FunctionalScheme
\ No newline at end of file
diff --git a/docs/src/api/utils.md b/docs/src/api/utils.md
index d8d944dac..8c886f111 100644
--- a/docs/src/api/utils.md
+++ b/docs/src/api/utils.md
@@ -6,6 +6,6 @@ CurrentModule = MethodOfLines
 
 ```@docs
 unitindex
-params
+ivs
 Idx
 ```
\ No newline at end of file
diff --git a/docs/src/devnotes.md b/docs/src/devnotes.md
index 54f252420..7bb5698f0 100644
--- a/docs/src/devnotes.md
+++ b/docs/src/devnotes.md
@@ -48,7 +48,7 @@ end
 # Note that indexmap is used along with the function `Idx` to create an equivalent index for the discrete form of `u`,
 # which may have a different number of dimensions to `II`
 function generate_central_difference_rules(II::CartesianIndex, s::DiscreteSpace, terms::Vector{<:Term}, indexmap::Dict)
-    rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in params(u, s)] for u in depvars]
+    rules = [[@rule Differential(x)(u) => second_order_central_difference(Idx(II, s, u, indexmap), s, u, x) for x in ivs(u, s)] for u in depvars]
 
     rules = reduce(vcat, rules)
 
diff --git a/docs/src/tutorials/PIDE.md b/docs/src/tutorials/PIDE.md
index 18c946da4..1dcd16b3b 100644
--- a/docs/src/tutorials/PIDE.md
+++ b/docs/src/tutorials/PIDE.md
@@ -56,7 +56,7 @@ To have an integral over the whole domain, be sure to wrap the integral in an au
 Due to a limitation, the whole domain integral needs to have the same arguments as the integrand, but is constant in x. To use it in an equation one dimension lower, use a boundary value like integral(t, 0.0)
 
 ```@example integrals2
-using MethodOfLines, ModelingToolkit, DomainSets
+using MethodOfLines, ModelingToolkit, DomainSets, OrdinaryDiffEq, Plots
 
 @parameters t, x
 @variables integrand(..) integral(..)
@@ -90,5 +90,6 @@ xdisc = sol[x]
 tdisc = sol[t]
 
 integralsol = sol[integral(t, x)]
-exact = [asf(t_) for t_ in tdisc, x_ in xdisc]
+
+heatmap(integralsol)
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/params.md b/docs/src/tutorials/params.md
index bd044c243..882c375e0 100644
--- a/docs/src/tutorials/params.md
+++ b/docs/src/tutorials/params.md
@@ -2,7 +2,7 @@
 
 We can also build up more complicated systems with multiple dependent variables and parameters as follows
 
-```@example params1
+```@example ivs1
 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
 
 @parameters t x
@@ -50,7 +50,7 @@ gif(anim, "plot.gif",fps=30)
 
 The system does not need to be re-discretized every time we want to plot with different parameters, the system can be remade with new parameters with `remake`. See the `ModelingToolkit.jl` [docs](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/#Algebraic-relations-and-structural-simplification) for more ways to manipulate a `prob` post discretization.
 
-```@example params2
+```@example ivs2
 using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
 
 @parameters t x
diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
index 38f0ccab1..1312db733 100644
--- a/src/MOL_discretization.jl
+++ b/src/MOL_discretization.jl
@@ -4,8 +4,8 @@ function interface_errors(depvars, indvars, discretization)
     for x in indvars
         @assert haskey(discretization.dxs, Num(x)) || haskey(discretization.dxs, x) "Variable $x has no step size"
     end
-    if !(typeof(discretization.advection_scheme) ∈ [UpwindScheme, WENOScheme])
-        throw(ArgumentError("Only `UpwindScheme()` and `WENOScheme()` are supported advection schemes."))
+    if !any(s -> discretization.advection_scheme isa s,  [UpwindScheme, FunctionalScheme])
+        throw(ArgumentError("Only `UpwindScheme()` and `FunctionalScheme()` are supported advection schemes. Got $(typeof(discretization.advection_scheme))."))
     end
     if !(typeof(discretization.disc_strategy) ∈ [ScalarizedDiscretization])
         throw(ArgumentError("Only `ScalarizedDiscretization()` are supported discretization strategies."))
@@ -107,7 +107,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         eqvar = interiormap.var[pde]
 
         # * Assumes that all variables in the equation have same dimensionality except edgevals
-        args = params(eqvar, s)
+        args = ivs(eqvar, s)
         indexmap = Dict([args[i] => i for i in 1:length(args)])
         if disc_strategy isa ScalarizedDiscretization
             # Generate the equations for the interior points
@@ -147,7 +147,8 @@ function SciMLBase.discretize(pdesys::PDESystem,discretization::MethodOfLines.MO
                     return prob
                 else
                     f = ODEFunction(pdesys, discretization, analytic=analytic, discretization.kwargs..., kwargs...)
-                    return ODEProblem(f, prob.u0, prob.tspan; discretization.kwargs..., kwargs...)
+
+                    return ODEProblem(f, prob.u0, prob.tspan, prob.p; discretization.kwargs..., kwargs...)
                 end
             end
         end
diff --git a/src/MOL_utils.jl b/src/MOL_utils.jl
index c81418190..2fcc11deb 100644
--- a/src/MOL_utils.jl
+++ b/src/MOL_utils.jl
@@ -64,13 +64,13 @@ end
 
 function get_gridloc(u, s)
     if isequal(operation(u), getindex)
+        # Remember arguments of getindex have u(t) first
         return _get_gridloc(s, arguments(u)...)
     else
         return (operation(u), [])
     end
 end
 
-
 function generate_function_from_gridlocs(analyticmap, gridlocs, s)
     is_t_first_map = Dict(map(s.ū) do u
         operation(u) => (findfirst(x -> isequal(s.time, x), arguments(u)) == 1)
@@ -84,14 +84,14 @@ function generate_function_from_gridlocs(analyticmap, gridlocs, s)
         is_t_first = is_t_first_map[uop]
         _f = analyticmap[opsmap[uop]]
         if is_t_first
-            return t -> _f(t, x̄...)
+            return (p, t) -> _f(p, t, x̄...)
         else
-            return t -> _f(x̄..., t)
+            return (p, t) -> _f(p, x̄..., t)
         end
     end
 
     f = (u0, p, t) -> map(fs_) do f_
-        f_(t)
+        f_(p, t)
     end
 
     return f
@@ -100,7 +100,7 @@ end
 function newindex(u_, II, s, indexmap)
     u = depvar(u_, s)
     args_ = remove(arguments(u_), s.time)
-    args = params(u, s)
+    args = ivs(u, s)
     is = map(enumerate(args_)) do (j, x)
         if haskey(indexmap, x)
             II[indexmap[x]]
@@ -135,8 +135,14 @@ end
 function chebyspace(N, dom)
     interval = dom.domain
     a, b = DomainSets.infimum(interval), DomainSets.supremum(interval)
-    x = reverse([(a + b) / 2 + (b - a) / 2 * cos(π * (2k - 1) / (2N)) for k in 1:N])
+    x = reverse([(a + b) / 2 + (b - a) / 2 * cospi((2k - 1) / (2N)) for k in 1:N])
     x[1] = a
     x[end] = b
     return dom.variables => x
 end
+
+@inline function sym_dot(a, b)
+    mapreduce((+), zip(a, b)) do (a_, b_)
+        a_ * b_
+    end
+end
diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl
index 1afac10c8..5b8296401 100644
--- a/src/MethodOfLines.jl
+++ b/src/MethodOfLines.jl
@@ -3,7 +3,8 @@ using LinearAlgebra
 using SciMLBase
 using DiffEqBase
 using ModelingToolkit
-using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states
+using ModelingToolkit: operation, istree, arguments, variable, get_metadata, get_states,
+parameters, defaults, varmap_to_vars
 using SymbolicUtils, Symbolics
 using Symbolics: unwrap, solve_for, expand_derivatives, diff2term, setname, rename, similarterm
 using SymbolicUtils: operation, arguments
@@ -60,6 +61,7 @@ include("system_parsing/pde_system_transformation.jl")
 include("discretization/interface_boundary.jl")
 
 # Schemes
+include("discretization/schemes/function_scheme/function_scheme.jl")
 include("discretization/schemes/centered_difference/centered_difference.jl")
 include("discretization/schemes/upwind_difference/upwind_difference.jl")
 include("discretization/schemes/half_offset_centred_difference.jl")
@@ -79,6 +81,6 @@ include("scalar_discretization.jl")
 include("MOL_discretization.jl")
 
 export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace
-export UpwindScheme, WENOScheme
+export UpwindScheme, WENOScheme, FunctionalScheme
 
 end
diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl
index 6f75312d3..5f97f7088 100644
--- a/src/discretization/discretize_vars.jl
+++ b/src/discretization/discretize_vars.jl
@@ -136,12 +136,12 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
         end
         if t === nothing
             uaxes = collect(axes(grid[x])[1] for x in arguments(u))
-            u => collect(first(@variables $sym[uaxes...]))
+            u => unwrap.(collect(first(@variables $sym[uaxes...])))
         elseif isequal(SymbolicUtils.arguments(u), [t])
-            u => fill(u, ()) #Create a 0-dimensional array
+            u => fill(safe_unwrap(u), ()) #Create a 0-dimensional array
         else
             uaxes = collect(axes(grid[x])[1] for x in remove(arguments(u), t))
-            u => collect(first(@variables $sym(t)[uaxes...]))
+            u => unwrap.(collect(first(@variables $sym(t)[uaxes...])))
         end
     end
 
@@ -149,28 +149,30 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
 end
 
 function Base.getproperty(s::DiscreteSpace, p::Symbol)
-    if p in [:ū, :x̄, :time, :args, :x2i, :i2x]
+    if p in [:ū, :x̄, :ps, :time, :args, :x2i, :i2x]
         getfield(s.vars, p)
     else
         getfield(s, p)
     end
 end
 
+params(s::DiscreteSpace) = s.ps
+
 get_grid_type(::DiscreteSpace{N,M,G}) where {N,M,G} = G
 
 prepare_dx(dx::Integer, xdomain, ::CenterAlignedGrid) = (xdomain[2] - xdomain[1])/(dx - 1)
 prepare_dx(dx::Integer, xdomain, ::EdgeAlignedGrid) = (xdomain[2] - xdomain[1])/dx
 prepare_dx(dx, xdomain, ::AbstractGrid) = dx
 
-nparams(::DiscreteSpace{N,M}) where {N,M} = N
+nivs(::DiscreteSpace{N,M}) where {N,M} = N
 nvars(::DiscreteSpace{N,M}) where {N,M} = M
 
 """
-    params(u, s::DiscreteSpace)
+    ivs(u, s::DiscreteSpace)
 
 Filter out the time variable and get the spatial variables of `u` in `s`.
 """
-params(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time)
+ivs(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time)
 Base.ndims(u, s::DiscreteSpace) = ndims(s.discvars[depvar(u, s)])
 
 Base.length(s::DiscreteSpace, x) = length(s.grid[x])
@@ -185,9 +187,9 @@ of `II` that corresponds to only the spatial arguments of `u`.
 """
 function Idx(II::CartesianIndex, s::DiscreteSpace, u, indexmap)
     # We need to construct a new index as indices may be of different size
-    length(params(u, s)) == 0 && return CartesianIndex()
-    !all(x -> haskey(indexmap, x), params(u, s)) && return II
-    is = [II[indexmap[x]] for x in params(u, s)]
+    length(ivs(u, s)) == 0 && return CartesianIndex()
+    !all(x -> haskey(indexmap, x), ivs(u, s)) && return II
+    is = [II[indexmap[x]] for x in ivs(u, s)]
 
     II = CartesianIndex(is...)
     return II
@@ -198,7 +200,7 @@ A function that returns what to replace independent variables with in boundary e
 """
 @inline function axiesvals(s::DiscreteSpace{N,M,G}, u_, x_, I) where {N,M,G}
     u = depvar(u_, s)
-    map(params(u, s)) do x
+    map(ivs(u, s)) do x
         if isequal(x, x_)
             x => (I[x2i(s, u, x)] == 1 ? first(s.axies[x]) : last(s.axies[x]))
         else
@@ -207,8 +209,8 @@ A function that returns what to replace independent variables with in boundary e
     end
 end
 
-gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in params(u, s)], s.Igrid[u])
-gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in params(u, s)]
+gridvals(s::DiscreteSpace{N}, u) where {N} = ndims(u, s) == 0 ? [] : map(y -> [x => s.grid[x][y.I[x2i(s, u, x)]] for x in ivs(u, s)], s.Igrid[u])
+gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 ? [] : [x => s.grid[x][I[x2i(s, u, x)]] for x in ivs(u, s)]
 
 
 varmaps(s::DiscreteSpace, depvars, II, indexmap) = [u => s.discvars[u][Idx(II, s, u, indexmap)] for u in depvars]
diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl
index 769826c03..43ff62dd8 100644
--- a/src/discretization/generate_bc_eqs.jl
+++ b/src/discretization/generate_bc_eqs.jl
@@ -1,5 +1,5 @@
 @inline function generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary::AbstractTruncatingBoundary)
-    args = params(depvar(boundary.u, s), s)
+    args = ivs(depvar(boundary.u, s), s)
     indexmap = Dict([args[i]=>i for i in 1:length(args)])
     push!(bceqs, generate_bc_eqs(s, boundaryvalfuncs, boundary, interiormap, indexmap))
 end
@@ -48,7 +48,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights
     # replace u(t,0) with u₁, etc
 
     u = depvar(u_, s)
-    args = params(u, s)
+    args = ivs(u, s)
     j = findfirst(isequal(x_), args)
     IIold = II
     # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
@@ -99,7 +99,7 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights
 
     # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary
     u = depvar(u_, s)
-    args = params(u, s)
+    args = ivs(u, s)
     j = findfirst(isequal(x_), args)
     IIold = II
     # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower
@@ -155,7 +155,7 @@ Pads the boundaries with extrapolation equations, extrapolated with 6th order la
 Reuses `central_difference` as this already dispatches the correct stencil, given a `DerivativeOperator` which contains the correct weights.
 """
 function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
-    args = params(u, s)
+    args = ivs(u, s)
     length(args) == 0 && return
 
     lowerextents, upperextents = interiormap.stencil_extents[pde]
@@ -167,7 +167,12 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
     for (j, x) in enumerate(args)
         ninterp = lowerextents[j] - vlower[j]
         I1 = unitindex(length(args), j)
+        bs = bcmap[operation(u)][x]
+        haslower, hasupper = haslowerupper(bs, x)
         while ninterp >= vlower[j]
+            if haslower
+                break
+            end
             for Il in (edge(interiormap, s, u, j, true) .+ (ninterp * I1,))
                 expr = central_difference(derivweights.boundary[x], Il, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc)
                 push!(eqmap[Il], expr)
@@ -176,6 +181,9 @@ function generate_extrap_eqs!(eqs, pde, u, s, derivweights, interiormap, bcmap)
         end
         ninterp = upperextents[j] - vupper[j]
         while ninterp >= vupper[j]
+            if hasupper
+                break
+            end
             for Iu in (edge(interiormap, s, u, j, false) .- (ninterp * I1,))
                 expr = central_difference(derivweights.boundary[x], Iu, s, filter_interfaces(bcmap[operation(u)][x]), (j, x), u, ufunc)
                 push!(eqmap[Iu], expr)
@@ -201,7 +209,7 @@ end
 
 @inline function generate_corner_eqs!(bceqs, s, interiormap, N, u)
     interior = interiormap.I[interiormap.pde[u]]
-    ndims(u, s) == 0 && return
+    N == 0 || ndims(u, s) == 0 && return
     sd(i, j) = selectdim(interior, j, i)
     domain = setdiff(s.Igrid[u], interior)
     II1 = unitindices(N)
diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl
index 379be5a89..e443b83b3 100644
--- a/src/discretization/generate_finite_difference_rules.jl
+++ b/src/discretization/generate_finite_difference_rules.jl
@@ -54,10 +54,12 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace,
         # Advection rules
         if derivweights.advection_scheme isa UpwindScheme
             advection_rules = generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
-        elseif derivweights.advection_scheme isa WENOScheme
-            advection_rules = generate_WENO_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+        elseif derivweights.advection_scheme isa FunctionalScheme
+            advection_rules = generate_advection_rules(derivweights.advection_scheme, II, s,
+                                                       depvars, derivweights, bmap, indexmap, terms)
             advection_rules = vcat(advection_rules,
-                                generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms; skip = [1]))
+                                generate_winding_rules(II, s, depvars, derivweights, bmap,
+                                                       indexmap, terms; skip = [1]))
         else
             error("Unsupported advection scheme $(derivweights.advection_scheme) encountered.")
         end
diff --git a/src/discretization/generate_ic_defaults.jl b/src/discretization/generate_ic_defaults.jl
index 8b0e10fe6..d922bfa62 100644
--- a/src/discretization/generate_ic_defaults.jl
+++ b/src/discretization/generate_ic_defaults.jl
@@ -6,7 +6,7 @@ function generate_ic_defaults(tconds, s, ::ScalarizedDiscretization)
                 throw(ArgumentError("Upper boundary condition $(ic.eq) on time variable is not supported, please use a change of variables `t => -τ` to make this an initial condition."))
             end
 
-            args = params(depvar(ic.u, s), s)
+            args = ivs(depvar(ic.u, s), s)
             indexmap = Dict([args[i]=>i for i in 1:length(args)])
             D = ic.order == 0 ? identity : (Differential(t)^ic.order)
             defaultvars = D.(s.discvars[depvar(ic.u, s)])
diff --git a/src/discretization/schemes/WENO/WENO.jl b/src/discretization/schemes/WENO/WENO.jl
index 17ba3fd9c..990cda6a4 100644
--- a/src/discretization/schemes/WENO/WENO.jl
+++ b/src/discretization/schemes/WENO/WENO.jl
@@ -1,37 +1,16 @@
-
 """
 Implements the WENO scheme of Jiang and Shu.
-
 Specified in https://repository.library.brown.edu/studio/item/bdr:297524/PDF/ (Page 8-9)
-
 Implementation *heavily* inspired by https://github.com/ranocha/HyperbolicDiffEq.jl/blob/84c2d882e0c8956457c7d662bf7f18e3c27cfa3d/src/finite_volumes/weno_jiang_shu.jl by H. Ranocha.
 """
-function weno(II::CartesianIndex, s::DiscreteSpace, wenoscheme::WENOScheme, bs, jx, u, dx::Number)
-    j, x = jx
-    ε = wenoscheme.epsilon
-
-    I1 = unitindex(ndims(u, s), j)
+function weno_f(u, p, t, x, dx)
+    ε = p[1]
 
-    udisc = s.discvars[u]
-
-    Im2 = bwrap(II - 2I1, bs, s, jx)
-    Im1 = bwrap(II - I1, bs, s, jx)
-    Ip1 = bwrap(II + I1, bs, s, jx)
-    Ip2 = bwrap(II + 2I1, bs, s, jx)
-    is = map(I -> I[j], [Im2, Im1, Ip1, Ip2])
-    for i in is
-        if i < 1
-            return nothing
-        elseif i > length(s, x)
-            return nothing
-        end
-    end
-
-    u_m2 = udisc[Im2]
-    u_m1 = udisc[Im1]
-    u_0 = udisc[II]
-    u_p1 = udisc[Ip1]
-    u_p2 = udisc[Ip2]
+    u_m2 = u[1]
+    u_m1 = u[2]
+    u_0 = u[3]
+    u_p1 = u[4]
+    u_p2 = u[5]
 
     γm1 = 1 / 10
     γm2 = 3 / 5
@@ -74,16 +53,15 @@ function weno(II::CartesianIndex, s::DiscreteSpace, wenoscheme::WENOScheme, bs,
     hp = wp1 * hp1 + wp2 * hp2 + wp3 * hp3
     hm = wm1 * hm1 + wm2 * hm2 + wm3 * hm3
 
-    return recursive_unwrap((hp - hm) / dx)
-end
-
-function weno(II::CartesianIndex, s::DiscreteSpace, b, jx, u, dx::AbstractVector)
-    @assert false "WENO scheme not implemented for nonuniform grids."
+    return (hp - hm) / dx
 end
 
 """
-This is a catch all ruleset, as such it does not use @rule.
+`WENOScheme` of Jiang and Shu
+## Keyword Arguments
+- `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensetive problems will benefit from a smaller value. It is defined as a functional scheme.
 """
-@inline function generate_WENO_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
-    return reduce(safe_vcat, [[(Differential(x))(u) => weno(Idx(II, s, u, indexmap), s, derivweights.advection_scheme, filter_interfaces(bcmap[operation(u)][x]), (x2i(s, u, x), x), u, s.dxs[x]) for x in params(u, s)] for u in depvars], init = [])
+function WENOScheme(epsilon = 1e-6)
+    boundary_f = [nothing, nothing]
+    return FunctionalScheme{5, 0}(weno_f, boundary_f, boundary_f, false, [epsilon], name = "WENO")
 end
diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl
index 05aab1dc9..b47d6dc62 100644
--- a/src/discretization/schemes/centered_difference/centered_difference.jl
+++ b/src/discretization/schemes/centered_difference/centered_difference.jl
@@ -23,7 +23,7 @@ function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u
         Itap = [bwrap(II + i * I1, bs, s, jx) for i in half_range(D.stencil_length)]
     end
     # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want.
-    return recursive_unwrap(dot(weights, ufunc(u, Itap, x)))
+    return sym_dot(weights, ufunc(u, Itap, x))
 end
 
 function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u, ufunc) where {T,N,Wind,DX<:AbstractVector}
@@ -48,7 +48,7 @@ function central_difference(D::DerivativeOperator{T,N,Wind,DX}, II, s, bs, jx, u
     end
     # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want.
 
-    return recursive_unwrap(dot(weights, ufunc(u, Itap, x)))
+    return sym_dot(weights, ufunc(u, Itap, x))
 end
 
 """
@@ -60,5 +60,5 @@ This is a catch all ruleset, as such it does not use @rule. Any even ordered der
         let orders = derivweights.orders[x]
             orders[iseven.(orders)]
         end
-    )] for x in params(u, s)], init = []) for u in depvars], init = [])
+    )] for x in ivs(u, s)], init = []) for u in depvars], init = [])
 end
diff --git a/src/discretization/schemes/function_scheme/function_scheme.jl b/src/discretization/schemes/function_scheme/function_scheme.jl
new file mode 100644
index 000000000..757289e72
--- /dev/null
+++ b/src/discretization/schemes/function_scheme/function_scheme.jl
@@ -0,0 +1,63 @@
+function get_f_and_taps(F::FunctionalScheme, II, s, bs, jx, u)
+    j, x = jx
+    # unit index in direction of the derivative
+    I1 = unitindex(ndims(u, s), j)
+    # offset is important due to boundary proximity
+    haslower, hasupper = haslowerupper(bs, x)
+
+    lower_point_count = length(F.lower)
+    upper_point_count = length(F.upper)
+
+    if (II[j] <= lower_point_count) & !haslower
+        f = F.lower[II[j]]
+        offset = 1 - II[j]
+        Itap = [II + (i + offset) * I1 for i in 0:(F.boundary_points-1)]
+    elseif (II[j] > (length(s, x) - upper_point_count)) & !hasupper
+        f = F.upper[length(s, x)-II[j]+1]
+        offset = length(s, x) - II[j]
+        Itap = [II + (i + offset) * I1 for i in (-F.boundary_points+1):1:0]
+    else
+        f = F.interior
+        Itap = [bwrap(II + i * I1, bs, s, jx) for i in half_range(F.interior_points)]
+    end
+
+    return f, Itap
+end
+
+function function_scheme(F::FunctionalScheme, II, s, bs, jx, u, ufunc)
+    j, x = jx
+    ndims(u, s) == 0 && return 0
+
+    f, Itap = get_f_and_taps(F, II, s, bs, jx, u)
+    if isnothing(f)
+        error("Scheme $(F.name) applied to $u in direction of $x at point $II is not defined.")
+    end
+    # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want.
+    u_disc = ufunc(u, Itap, x)
+    ps = vcat(F.ps, params(s))
+    t = s.time
+    itap = map(I -> I[j], Itap)
+    discx = @view s.grid[x][itap]
+    dx = s.dxs[x]
+    if F.is_nonuniform
+        if dx isa AbstractVector
+            dx = @views dx[itap[1:end-1]]
+        end
+    elseif dx isa AbstractVector
+        error("Scheme $(F.name) not implemented for nonuniform dxs.")
+    end
+
+    return f(u_disc, ps, t, discx, dx)
+end
+
+@inline function generate_advection_rules(F::FunctionalScheme, II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
+    central_ufunc(u, I, x) = s.discvars[u][I]
+    return reduce(safe_vcat, [[(Differential(x))(u) =>
+                                        function_scheme(F,
+                                                        Idx(II, s, u, indexmap), s,
+                                                        filter_interfaces(bcmap[operation(u)][x]),
+                                                        (x2i(s, u, x), x), u,
+                                                        central_ufunc)
+                                for x in ivs(u, s)]
+                               for u in depvars], init=[])
+end
diff --git a/src/discretization/schemes/half_offset_centred_difference.jl b/src/discretization/schemes/half_offset_centred_difference.jl
index 2ae1a4357..8b6a00a8f 100644
--- a/src/discretization/schemes/half_offset_centred_difference.jl
+++ b/src/discretization/schemes/half_offset_centred_difference.jl
@@ -61,5 +61,5 @@ function half_offset_centered_difference(D, II, s, bs, jx, u, ufunc)
     ndims(u, s) == 0 && return Num(0)
     j, x = jx
     weights, Itap = get_half_offset_weights_and_stencil(D, II, s, bs, u, jx)
-    return recursive_unwrap(dot(weights, ufunc(u, Itap, x)))
+    return sym_dot(weights, ufunc(u, Itap, x))
 end
diff --git a/src/discretization/schemes/integral_expansion/integral_expansion.jl b/src/discretization/schemes/integral_expansion/integral_expansion.jl
index 8b0e5c82d..b6c60ec68 100644
--- a/src/discretization/schemes/integral_expansion/integral_expansion.jl
+++ b/src/discretization/schemes/integral_expansion/integral_expansion.jl
@@ -10,7 +10,7 @@ function _euler_integral(II, s, jx, u, ufunc, dx::Number) #where {T,N,Wind,DX<:N
     Itap = [II - I1, II]
     weights = [dx / 2, dx / 2]
 
-    return dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx)
+    return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx)
 end
 
 # Nonuniform dx
@@ -25,7 +25,7 @@ function _euler_integral(II, s, jx, u, ufunc, dx::AbstractVector) #where {T,N,Wi
     Itap = [II - I1, II]
     weights = fill(dx[II[j]-1] / 2, 2)
 
-    return dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx)
+    return sym_dot(weights, ufunc(u, Itap, x)) + _euler_integral(II - I1, s, jx, u, ufunc, dx)
 end
 
 function euler_integral(II, s, jx, u, ufunc)
@@ -52,14 +52,14 @@ end
     ufunc(u, I, x) = s.discvars[u][I]
 
     eulerrules = reduce(safe_vcat, [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], Num(x)))(u) => euler_integral(Idx(II, s, u, indexmap), s, (x2i(s, u, x), x), u, ufunc)
-                                for x in params(u, s)]
+                                for x in ivs(u, s)]
                                for u in depvars], init = [])
     return eulerrules
 end
 
 function wd_integral_Idx(II::CartesianIndex, s::DiscreteSpace, u, x, indexmap)
     # We need to construct a new index as indices may be of different size
-    length(params(u, s)) == 0 && return CartesianIndex()
+    length(ivs(u, s)) == 0 && return CartesianIndex()
     # A hack using the boundary value re-indexing function to get an index that will work
     u_ = substitute(u, [x => s.axies[x][end]])
     II = newindex(u_, II, s, indexmap)
@@ -71,7 +71,7 @@ end
     wholedomainrules = reduce(safe_vcat,
                               [[Integral(x in DomainSets.ClosedInterval(s.vars.intervals[x][1], s.vars.intervals[x][2]))(u) =>
                                     whole_domain_integral(wd_integral_Idx(II, s, u, x, indexmap), s, (x2i(s, u, x), x), u, ufunc)
-                                for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), params(u, s))]
+                                for x in filter(x -> (!haskey(indexmap, x) | isequal(x, bvar)), ivs(u, s))]
                                for u in depvars],
                               init = [])
     return wholedomainrules
diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
index d97ca7886..7360c7acd 100644
--- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
+++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
@@ -51,34 +51,34 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace,
     interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, I, s, bs, u, jx) for I in outerstencil]
 
     # map variables to symbolically inerpolated/extrapolated expressions
-    map_vars_to_interpolated(stencil, weights) = [v => dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars]
+    map_vars_to_interpolated(stencil, weights) = [v => sym_dot(weights, s.discvars[v][interface_wrap(stencil)]) for v in depvars]
 
     # Map parameters to interpolated values. Using simplistic extrapolation/interpolation for now as grids are uniform
     #TODO: make this more efficient
-    map_params_to_interpolated(stencil, weights) = safe_vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(interface_wrap(stencil), (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])])
+    map_ivs_to_interpolated(stencil, weights) = safe_vcat([x => dot(weights, getindex.((s.grid[x],), getindex.(interface_wrap(stencil), (j,))))], [s.x̄[k] => s.grid[s.x̄[k]][II[k]] for k in setdiff(1:N, [j])])
 
     # Take the inner finite difference
-    inner_difference = [dot(inner_weights, s.discvars[u][interface_wrap(inner_stencil)]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil]
+    inner_difference = [sym_dot(inner_weights, s.discvars[u][interface_wrap(inner_stencil)]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil]
 
     # Symbolically interpolate the multiplying expression
 
 
     interpolated_expr = map(interp_weights_and_stencil) do (weights, stencil)
-        substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_params_to_interpolated(stencil, weights))
+        substitute(substitute(expr, map_vars_to_interpolated(stencil, weights)), map_ivs_to_interpolated(stencil, weights))
     end
 
     # multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference
-    return dot(outerweights, inner_difference .* interpolated_expr)
+    return sym_dot(outerweights, inner_difference .* interpolated_expr)
 end
 
 @inline function generate_nonlinlap_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
-    rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in params(u, s)]) for u in depvars], init = [])
+    rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in ivs(u, s)]) for u in depvars], init = [])
 
-    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = []))
+    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in ivs(u, s)]) for u in depvars], init = []))
 
-    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u) / ~a)) => cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = []))
+    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u) / ~a)) => cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in ivs(u, s)]) for u in depvars], init = []))
 
-    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(~b..., ~c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in params(u, s)]) for u in depvars], init = []))
+    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(~b..., ~c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in ivs(u, s)]) for u in depvars], init = []))
 
     nonlinlap_rules = []
     for t in terms
diff --git a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl
index 013a752e6..9fcb666f1 100644
--- a/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl
+++ b/src/discretization/schemes/spherical_laplacian/spherical_laplacian.jl
@@ -36,13 +36,13 @@ end
 
 @inline function generate_spherical_diffusion_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
     rules = reduce(safe_vcat, [vec([@rule *(~~a, 1 / (r^2), ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) => *(~a..., spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u), ~b...)
-                               for r in params(u, s)]) for u in depvars], init = [])
+                               for r in ivs(u, s)]) for u in depvars], init = [])
 
     rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(*(~~a, $(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e)), ~~b), (r^2)) => *(~a..., ~b..., spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u))
-                                           for r in params(u, s)]) for u in depvars], init = []))
+                                           for r in ivs(u, s)]) for u in depvars], init = []))
 
     rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule /(($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), (r^2)) => spherical_diffusion(*(~c..., ~d..., ~e..., Num(1)), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][r]), depvars, r, u)
-                                           for r in params(u, s)]) for u in depvars], init = []))
+                                           for r in ivs(u, s)]) for u in depvars], init = []))
 
     spherical_diffusion_rules = []
     for t in terms
diff --git a/src/discretization/schemes/upwind_difference/upwind_difference.jl b/src/discretization/schemes/upwind_difference/upwind_difference.jl
index 1feab1bdd..92626876c 100644
--- a/src/discretization/schemes/upwind_difference/upwind_difference.jl
+++ b/src/discretization/schemes/upwind_difference/upwind_difference.jl
@@ -65,7 +65,7 @@ function upwind_difference(d::Int, II::CartesianIndex, s::DiscreteSpace, bs, der
     #@show D.stencil_coefs, D.stencil_length, D.boundary_stencil_length, D.boundary_point_count
     # unit index in direction of the derivative
     weights, Itap = _upwind_difference(D, II, s, bs, ispositive, u, jx)
-    return recursive_unwrap(dot(weights, ufunc(u, Itap, x)))
+    return sym_dot(weights, ufunc(u, Itap, x))
 end
 
 function upwind_difference(expr, d::Int, II::CartesianIndex, s::DiscreteSpace, bs, depvars, derivweights, (j, x), u, central_ufunc, indexmap)
@@ -84,14 +84,14 @@ end
             let orders = derivweights.orders[x]
                 setdiff(orders[isodd.(orders)], skip)
             end
-        )] for x in params(u, s)], init = []) for u in depvars], init = []),
+        )] for x in ivs(u, s)], init = []) for u in depvars], init = []),
 
         #Catch division and multiplication, see issue #1
         reduce(safe_vcat, [reduce(safe_vcat, [[@rule /(*(~~a, $(Differential(x)^d)(u), ~~b), ~c) => upwind_difference(*(~a..., ~b...) / ~c, d, Idx(II, s, u, indexmap), s, filter_interfaces(bcmap[operation(u)][x]), depvars, derivweights, (x2i(s, u, x), x), u, wind_ufunc, indexmap) for d in (
             let orders = derivweights.orders[x]
                 setdiff(orders[isodd.(orders)], skip)
             end
-        )] for x in params(u, s)], init = []) for u in depvars], init = [])
+        )] for x in ivs(u, s)], init = []) for u in depvars], init = [])
     )
 
     wind_rules = []
@@ -106,7 +106,7 @@ end
     end
 
     return safe_vcat(wind_rules, vec(mapreduce(safe_vcat, depvars, init = []) do u
-        mapreduce(safe_vcat, params(u, s), init = []) do x
+        mapreduce(safe_vcat, ivs(u, s), init = []) do x
             j = x2i(s, u, x)
             let orders = setdiff(derivweights.orders[x], skip)
                 oddorders = orders[isodd.(orders)]
diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl
index 3d7c9035f..e7efbcaf9 100644
--- a/src/interface/scheme_types.jl
+++ b/src/interface/scheme_types.jl
@@ -2,28 +2,100 @@ abstract type AbstractScheme end
 
 struct UpwindScheme <: AbstractScheme
     order
-    function UpwindScheme(approx_order = 1)
+    function UpwindScheme(approx_order=1)
         return new(approx_order)
     end
 end
 
-extent(scheme::UpwindScheme, dorder) = dorder+scheme.order-1
+extent(scheme::UpwindScheme, dorder) = dorder + scheme.order - 1
+
+# Functional Schemes
 
 """
-`WENOScheme` of Jiang and Shu
-## Keyword Arguments
-- `epsilon`: A quantity used to prevent vanishing denominators in the scheme, defaults to `1e-6`. More sensetive problems will benefit from a smaller value.
+# FunctionalScheme
+```
+F = FunctionalScheme{interior_points, boundary_points}(interior, lower, upper, is_nonuniform, parameters; name)
+```
+A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary.
+
+`lower` and `upper` should be vectors of functions. In general, `upper` and `lower` must be at least `floor(interior_points/2)` long. Where you have no good approximation for a derivative at the boundary, you can use `nothing` as a placeholder. MethodOfLines will then attempt to use an extrapolation here where nessesary. Be warned that this can lead to instability.
+
+The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if `boundary_points = 3`, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps.
+
+The functions making up the scheme take the following inputs:
+
+Functions must be of the form `f(u, p, t, deriv_iv, d_iv)`.
+
+Where the function would branch on a value of u, p, t, or d_iv, use IfElse.ifelse instead of standard conditionals.
+
+For the interior, `u` takes a vector of dependent variable values in the direction of the derivative
+of length `interior_points`. `interior_points` must be odd, as this function defines the derivative at the center of the input points.
+
+For the lower and upper boundaries, `u` takes a vector of dependent variable values of length `boundary_points`. This will be the `boundary_points` number of points closest to the lower and upper boundary respectively.
+`p` will take all parameter values in the order specified in the PDESystem, with the scheme's parameters prepended to the list.
+
+`deriv_iv` takes a vector of independent variable values of the same support as for `u`, for the independent variable in the direction of the derivative.
+
+If `is_nonuniform` is false, `d_iv` will take a scalar value of the stepsize between the points used to call the function in `u` and `deriv_iv`.
+
+If `is_nonuniform` is true, the scheme must be able to accept `d_iv` as a vector of stepsizes between the points used to call the function in `u` and `deriv_iv`, therefore of length `length(u)-1`. A method should also be defined for the case where `d_iv` is a scalar, in which case the stepsizes are assumed to be uniform.
 """
-struct WENOScheme <: AbstractScheme
-    epsilon
-    function WENOScheme(epsilon = 1e-6)
-        new(epsilon)
-    end
+struct FunctionalScheme{F, V1, V2, V3} <: AbstractScheme
+    """
+    The function that defines the scheme on the interior.
+    """
+    interior::F
+    """
+    The vector of functions that defines the scheme near the lower boundary.
+    """
+    lower::V1
+    """
+    The vector functions that defines the scheme near the upper boundary.
+    """
+    upper::V2
+    """
+    The number of interior points that the interior scheme takes as input
+    """
+    interior_points::Int
+    """
+    The number of boundary points that the lower and upper schemes take as input.
+    """
+    boundary_points::Int
+    """
+    Whether this scheme takes grid steps as input for the nonuniform case.
+    """
+    is_nonuniform::Bool
+    """
+    parameters for the scheme, should be a vector of numbers.
+    """
+    ps::V3
+    """
+    The name of the scheme.
+    """
+    name::String
+end
+
+function FunctionalScheme{ips, bps}(interior, lower, upper, is_nonuniform = false, ps = []; name = "FunctionalScheme") where {ips, bps}
+    @assert ips % 2 == 1 "interior_points must be odd."
+    FunctionalScheme{typeof(interior), typeof(lower),
+                     typeof(upper), typeof(ps)}(interior, lower, upper,
+                                                ips, bps, is_nonuniform, ps, name)
+end
+
+function extent(scheme::FunctionalScheme, dorder)
+    @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes."
+    lower = length(findall(isnothing, scheme.lower))
+    upper = length(findall(isnothing, scheme.upper))
+    @assert lower == upper "Scheme must have symmetric extent; same number of placeholders in lower and upper boundary functions."
+    return lower
 end
 
-function extent(::WENOScheme, dorder)
-    @assert dorder == 1 "Only first order spatial derivatives are implemented for WENO schemes."
-    return 2
+function lower_extent(scheme::FunctionalScheme, dorder)
+    @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes."
+    return length(findall(isnothing, scheme.lower))
 end
 
-# Note: This type and its subtypes will become important later with the stencil interfaces as we will need to dispatch on derivative order and approximation order
+function upper_extent(scheme::FunctionalScheme, dorder)
+    @assert dorder == 1 "Only first order spatial derivatives are implemented for functional schemes."
+    return length(findall(isnothing, scheme.upper))
+end
diff --git a/src/interface/solution/timedep.jl b/src/interface/solution/timedep.jl
index 8bad97a41..f0c731ed6 100644
--- a/src/interface/solution/timedep.jl
+++ b/src/interface/solution/timedep.jl
@@ -50,10 +50,10 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T},
         return SciMLBase.PDETimeSeriesSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata),
             typeof(sol),typeof(sol.errors),typeof(sol.t),typeof(ivgrid),
             typeof(ivs),typeof(pdesys.dvs),typeof(sol.prob),typeof(sol.alg),
-            typeof(interp)}(umap, sol, sol.errors, sol.t, ivgrid, ivs,
+            typeof(interp), typeof(sol.stats)}(umap, sol, sol.errors, sol.t, ivgrid, ivs,
             pdesys.dvs, metadata, sol.prob, sol.alg,
             interp, sol.dense, sol.tslocation,
-            sol.retcode)
+            sol.retcode, sol.stats)
     catch e
         rethrow(e)
         return sol, e
diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl
index fbd5056c0..3fce94195 100644
--- a/src/interface/solution/timeindep.jl
+++ b/src/interface/solution/timeindep.jl
@@ -30,7 +30,7 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada
 
     return SciMLBase.PDENoTimeSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata),
         typeof(sol),typeof(ivgrid),typeof(ivs),typeof(pdesys.dvs),typeof(sol.prob),typeof(sol.alg),
-        typeof(interp)}(umap, sol, ivgrid, ivs,
+        typeof(interp), typeof(sol.stats)}(umap, sol, ivgrid, ivs,
         pdesys.dvs, metadata, sol.prob, sol.alg,
-        interp, sol.retcode)
+        interp, sol.retcode, sol.stats)
 end
diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl
index 2ecb1f840..bc5e7d207 100644
--- a/src/scalar_discretization.jl
+++ b/src/scalar_discretization.jl
@@ -47,18 +47,23 @@ function generate_system(alleqs, bceqs, ics, discvars, defaults, ps, tspan, meta
     alleqs = vcat(alleqs, unique(bceqs))
     alldepvarsdisc = vec(reduce(vcat, vec(unique(reduce(vcat, vec.(values(discvars)))))))
     # Finalize
+    # if haskey(metadata.disc.kwargs, :checks)
+    #     checks = metadata.disc.kwargs[:checks]
+    # else
+         checks = true
+    # end
     try
         if t === nothing
             # At the time of writing, NonlinearProblems require that the system of equations be in this form:
             # 0 ~ ...
             # Thus, before creating a NonlinearSystem we normalize the equations s.t. the lhs is zero.
             eqs = map(eq -> 0 ~ eq.rhs - eq.lhs, alleqs)
-            sys = NonlinearSystem(eqs, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata)
+            sys = NonlinearSystem(eqs, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata, checks = checks)
             return sys, nothing
         else
             # * In the end we have reduced the problem to a system of equations in terms of Dt that can be solved by an ODE solver.
 
-            sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata)
+            sys = ODESystem(alleqs, t, alldepvarsdisc, ps, defaults=defaults, name=name, metadata=metadata, checks = checks)
             return sys, tspan
         end
     catch e
diff --git a/src/system_parsing/bcs/parse_boundaries.jl b/src/system_parsing/bcs/parse_boundaries.jl
index db451dfb3..dbaac9b17 100644
--- a/src/system_parsing/bcs/parse_boundaries.jl
+++ b/src/system_parsing/bcs/parse_boundaries.jl
@@ -179,9 +179,9 @@ function generate_boundary_matching_rules(v, orders)
     upperboundary(x) = v.intervals[x][2]
 
     # Rules to match boundary conditions on the lower boundaries
-    lower = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, lowerboundary(x)) for x in all_params(u, v)]) for u in v.ū])
+    lower = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, lowerboundary(x)) for x in all_ivs(u, v)]) for u in v.ū])
 
-    upper = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, upperboundary(x)) for x in all_params(u, v)]) for u in v.ū])
+    upper = Dict([operation(u) => Dict([x => _boundary_rules(v, orders, u, x, upperboundary(x)) for x in all_ivs(u, v)]) for u in v.ū])
 
     return (lower, upper)
 end
diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl
index 1fdfccf7a..49790f600 100644
--- a/src/system_parsing/interior_map.jl
+++ b/src/system_parsing/interior_map.jl
@@ -59,7 +59,7 @@ end
 function generate_interior(lower, upper, u, s, disc::MOLFiniteDifference{G,D}) where {G, D<:ScalarizedDiscretization}
     args = remove(arguments(u), s.time)
 
-    ret = s.Igrid[u][[(1+lower[x2i(s, u, x)]:length(s.grid[x])-upper[x2i(s, u, x)]) for x in args]...]
+    ret = collect(s.Igrid[u][[(1+lower[x2i(s, u, x)]:length(s.grid[x])-upper[x2i(s, u, x)]) for x in args]...])
     return ret
 end
 
diff --git a/src/system_parsing/variable_map.jl b/src/system_parsing/variable_map.jl
index 7ec4a0646..8e3fd3ed3 100644
--- a/src/system_parsing/variable_map.jl
+++ b/src/system_parsing/variable_map.jl
@@ -1,6 +1,7 @@
 struct VariableMap
     ū
     x̄
+    ps
     time
     intervals
     args
@@ -9,8 +10,14 @@ struct VariableMap
     i2x
 end
 
-function VariableMap(eqs, depvars, domain, time)
+function VariableMap(eqs, depvars, domain, time, ps = [])
+    if ps isa SciMLBase.NullParameters
+        ps = []
+    end
     time = safe_unwrap(time)
+    ps = map(ps) do p
+        safe_unwrap(p.first)
+    end
     depvar_ops = get_ops(depvars)
     # Get all dependent variables in the correct type
     alldepvars = get_all_depvars(eqs, depvar_ops)
@@ -27,11 +34,11 @@ function VariableMap(eqs, depvars, domain, time)
     args = [operation(u) => arguments(u) for u in ū]
     x̄2dim = [x̄[i] => i for i in 1:nspace]
     dim2x̄ = [i => x̄[i] for i in 1:nspace]
-    return VariableMap(ū, x̄, time, Dict(intervals), Dict(args), depvar_ops, Dict(x̄2dim), Dict(dim2x̄))
+    return VariableMap(ū, x̄, ps, time, Dict(intervals), Dict(args), depvar_ops, Dict(x̄2dim), Dict(dim2x̄))
 end
 
-VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, disc.time)
-VariableMap(pdesys::PDESystem, t) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, t)
+VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, disc.time, pdesys.ps)
+VariableMap(pdesys::PDESystem, t) = VariableMap(pdesys.eqs, pdesys.dvs, pdesys.domain, t, pdesys.ps)
 
 function update_varmap!(v, newdv)
     push!(v.ū, newdv)
@@ -40,13 +47,13 @@ function update_varmap!(v, newdv)
 end
 
 
-params(u, v::VariableMap) = remove(v.args[operation(u)], v.time)
+ivs(u, v::VariableMap) = remove(v.args[operation(u)], v.time)
 
-Base.ndims(u, v::VariableMap) = length(params(u, v))
+Base.ndims(u, v::VariableMap) = length(ivs(u, v))
 
 all_ivs(v::VariableMap) = v.time === nothing ? v.x̄ : v.x̄ ∪ [v.time]
 
-all_params(u, v::VariableMap) = v.args[operation(u)]
+all_ivs(u, v::VariableMap) = v.args[operation(u)]
 
 depvar(u, v::VariableMap) = operation(u)(v.args[operation(u)]...)
 
@@ -54,7 +61,7 @@ x2i(v::VariableMap, u, x) = findfirst(isequal(x), remove(v.args[operation(u)], v
 
 @inline function axiesvals(v::VariableMap, u_, x_, I)
     u = depvar(u_, v)
-    map(params(u, v)) do x
+    map(ivs(u, v)) do x
         x => (I[x2i(v, u, x)] == 1 ? v.intervals[x][1] : v.intervals[x][2])
     end
 end
diff --git a/test/components/ODEFunction_test.jl b/test/components/ODEFunction_test.jl
index ecb701126..b0cf165d8 100644
--- a/test/components/ODEFunction_test.jl
+++ b/test/components/ODEFunction_test.jl
@@ -1,6 +1,6 @@
 using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, Test
 
-@parameters x t
+@parameters x t a
 @variables u(..)
 Dx = Differential(x)
 Dt = Differential(t)
@@ -9,22 +9,22 @@ x_max = 1.0
 t_min = 0.0
 t_max = 6.0
 
-analytic_u(t, x) = x / (t + 1)
+analytic_u(p, t, x) = x / (t + p[1])
 
-eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x))
+eq = Dt(u(t, x)) ~ -a*u(t, x) * Dx(u(t, x))
 
 bcs = [u(0, x) ~ x,
-    u(t, x_min) ~ analytic_u(t, x_min),
-    u(t, x_max) ~ analytic_u(t, x_max)]
+    u(t, x_min) ~ analytic_u([1], t, x_min),
+    u(t, x_max) ~ analytic_u([1], t, x_max)]
 
 domains = [t ∈ Interval(t_min, t_max),
     x ∈ Interval(x_min, x_max)]
 
-@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])
+@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [a => 1.0], analytic_func = [u(t, x) => analytic_u])
 
 disc = MOLFiniteDifference([x => 30], t, advection_scheme=WENOScheme())
 
-prob = discretize(pdesys, disc; analytic = [u(t, x) => analytic_u])
+prob = discretize(pdesys, disc; analytic = pdesys.analytic_func)
 
 sol = solve(prob, FBDF())
 
diff --git a/test/components/finite_diff_schemes.jl b/test/components/finite_diff_schemes.jl
index 256d62b14..66fb5f5e7 100644
--- a/test/components/finite_diff_schemes.jl
+++ b/test/components/finite_diff_schemes.jl
@@ -1,6 +1,6 @@
 using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, LinearAlgebra
 
-@parameters x, t 
+@parameters x, t
 @variables u(..)
 
 Dx(d) = Differential(x)^d
@@ -10,7 +10,7 @@ Dt = Differential(t)
 t_min= 0.
 t_max = 2.0
 x_min = 0.
-x_max = 20.0 
+x_max = 20.0
 
 dx = 1.0
 
@@ -28,7 +28,7 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
 
             @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])
 
-            # Test centered order 
+            # Test centered order
             disc = MOLFiniteDifference([x=>dx], t; approx_order=a)
 
             depvar_ops = map(x->operation(x.val),pdesys.depvars)
@@ -44,7 +44,7 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
             s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)
 
             derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)
-            
+
             #@show pde.rhs, operation(pde.rhs), arguments(pde.rhs)
             for II in s.Igrid[s.ū[1]][2:end-1]
                 #II = s.Igrid[end-1]
@@ -67,16 +67,16 @@ domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
 end
 
 @testset "Test 01: Nonlinear Diffusion discretization" begin
-    
+
     pde = Dt(u(t,x)) ~ Dx(1)(u(t,x))
     bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)]
 
     @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])
 
-    # Test centered order 
-    
+    # Test centered order
+
     depvar_ops = map(x->operation(x.val),pdesys.depvars)
-    
+
     depvars_lhs = MethodOfLines.get_depvars(pde.lhs, depvar_ops)
     depvars_rhs = MethodOfLines.get_depvars(pde.rhs, depvar_ops)
     depvars = collect(depvars_lhs ∪ depvars_rhs)
@@ -84,15 +84,15 @@ end
     # ignore if the only argument is [t]
     indvars = first(Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars))))
     x̄ = first(Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars))))
-    
+
     for order in [2]
         disc = MOLFiniteDifference([x=>dx], t; approx_order=order)
         s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)
 
         derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)
-        
+
         ufunc(u, I, x) = s.discvars[u][I]
-        #TODO Test Interpolation of params
+        #TODO Test Interpolation of ivs
         # Test simple case
         for II in s.Igrid[s.ū[1]][2:end-1]
             expr = MethodOfLines.cartesian_nonlinear_laplacian((1~1).lhs, II, derivweights, s, depvars, x, u(t,x))
@@ -103,14 +103,14 @@ end
 end
 
 @testset "Test 02: Spherical Diffusion discretization" begin
-    
+
     pde  = Dt(u(t,x)) ~ 1/x^2 * Dx(1)(x^2 * Dx(1)(u(t,x)))
 
     bcs = [u(0,x) ~ cos(x), u(t,0) ~ exp(-t), u(t,Float64(π)) ~ -exp(-t)]
 
     @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x)])
 
-    # Test centered order 
+    # Test centered order
     disc = MOLFiniteDifference([x=>dx], t; approx_order=2)
 
     depvar_ops = map(x->operation(x.val),pdesys.depvars)
@@ -126,11 +126,11 @@ end
     s = MethodOfLines.DiscreteSpace(domains, depvars, x̄, disc)
 
     derivweights = MethodOfLines.DifferentialDiscretizer(pdesys, s, disc)
-    
+
     for II in s.Igrid[s.ū[1]][2:end-1]
-        #TODO Test Interpolation of params
+        #TODO Test Interpolation of ivs
         expr = MethodOfLines.spherical_diffusion((1~1).lhs, II, derivweights, s, depvars, x, u(t,x))
         #@show II, expr
     end
 
-end
\ No newline at end of file
+end
diff --git a/test/components/function_scheme.jl b/test/components/function_scheme.jl
new file mode 100644
index 000000000..fc5c1facb
--- /dev/null
+++ b/test/components/function_scheme.jl
@@ -0,0 +1,11 @@
+@testset "User defined scheme, 3 interior points" begin
+    f_interior = (u, p, t, x, dx) -> IfElse.ifelse(u[2] < 0, (u[2] - u[1]) / dx, (u[3] - u[2]) / dx)
+    f_lower = (u, p, t, x, dx) -> (u[2] - u[1]) / dx
+    f_upper = (u, p, t, x, dx) -> (u[2] - u[1]) / dx
+
+    scheme = FunctionScheme{3, 2}(f_interior, f_lower, f_upper, false)
+
+    disc = MOLFiniteDifference([x => dx], t, advection_scheme=scheme)
+
+
+end
diff --git a/test/pde_systems/MOL_1D_HigherOrder.jl b/test/pde_systems/MOL_1D_HigherOrder.jl
index 83fb6f479..a21aeab5d 100644
--- a/test/pde_systems/MOL_1D_HigherOrder.jl
+++ b/test/pde_systems/MOL_1D_HigherOrder.jl
@@ -114,7 +114,7 @@ end
 
     sol = solve(prob,FBDF(),saveat=0.1)
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     xs = sol[x]
     ts = sol[t]
diff --git a/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl b/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl
index 2c569372a..bf925fae7 100644
--- a/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl
+++ b/test/pde_systems/MOL_1D_NonLinear_Diffusion.jl
@@ -46,7 +46,7 @@ using ModelingToolkit: Differential
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -161,7 +161,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -222,7 +222,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -337,7 +337,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32()) # TODO: check warnings
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -394,7 +394,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32()) # TODO: check warnings
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
 
     # Test against exact solution
@@ -453,7 +453,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
-    @test_broken sol.retcode == :Success
+    @test_broken sol.retcode == SciMLBase.ReturnCode.Success
     fail # don't run the rest of the test
     # Test against exact solution
     x_disc = sol[x]
@@ -510,7 +510,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
-    @test_broken sol.retcode == :Success
+    @test_broken sol.retcode == SciMLBase.ReturnCode.Success
     fail #don't run later tests
     # Test against exact solution
     x_disc = sol[x]
diff --git a/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl b/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl
index db4959296..1f75c76ce 100644
--- a/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl
+++ b/test/pde_systems/MOL_1D_NonLinear_Diffusion_NonUniform.jl
@@ -52,7 +52,7 @@ using StableRNGs
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -173,7 +173,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -236,7 +236,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -353,7 +353,7 @@ end
     # Solution of the ODE system
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -414,7 +414,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 
     # Test against exact solution
     x_disc = sol[x]
@@ -476,7 +476,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test_broken sol.retcode == :Success
+    @test_broken sol.retcode == SciMLBase.ReturnCode.Success
     fail
     # Test against exact solution
     x_disc = sol[x]
@@ -538,7 +538,7 @@ end
     using OrdinaryDiffEq
     sol = solve(prob, Rosenbrock32())
 
-    @test_broken sol.retcode == :Success
+    @test_broken sol.retcode == SciMLBase.ReturnCode.Success
     fail
     # Test against exact solution
     x_disc = sol[x]
diff --git a/test/pde_systems/MOL_Mixed_Deriv.jl b/test/pde_systems/MOL_Mixed_Deriv.jl
index f37dc4e09..dde8d1a29 100644
--- a/test/pde_systems/MOL_Mixed_Deriv.jl
+++ b/test/pde_systems/MOL_Mixed_Deriv.jl
@@ -10,7 +10,7 @@ using ModelingToolkit: Differential
     Dxx = Differential(x)^2
     Dtx = Differential(t)*Differential(x)
 
-    asf(t, x) = sinpi(2*(t + (1+sqrt(5))*x/2)) + cospi(2*(t + (1-sqrt(5))*x/2))
+    assf(t, x) = sinpi(2*(t + (1+sqrt(5))*x/2)) + cospi(2*(t + (1-sqrt(5))*x/2))
     aDtf(t, x) = 2pi*cospi(2*(t + (1+sqrt(5))*x/2)) - 2pi*sinpi(2*(t + (1-sqrt(5))*x/2))
     # Where asf(t, x) ~ 0, NonlinearSolved
     xmin = -0.1118033987645
@@ -18,7 +18,7 @@ using ModelingToolkit: Differential
 
     eq  = [Dtt(u(t, x)) + Dtx(u(t, x)) - Dxx(u(t, x)) ~ Dxx(x)]
 
-    bcs = [u(1e-9, x) ~ asf(1e-9, x),
+    bcs = [u(1e-9, x) ~ assf(1e-9, x),
            Dt(u(1e-9, x)) ~ aDtf(1e-9, x),
            u(t, xmin) ~ 0,
            u(t, xmax) ~ 0]
@@ -39,7 +39,7 @@ using ModelingToolkit: Differential
     tdisc = sol[t]
     usol = sol[u(t,x)]
 
-    asol = [asf(t, x) for t in tdisc, x in xdisc]
+    asol = [assf(t, x) for t in tdisc, x in xdisc]
     @test_broken usol ≈ asol atol = 1e-3
 end
 
@@ -67,5 +67,5 @@ end
 
     prob = discretize(pdesys, discretization, advection_scheme = WENOScheme())
     sol = solve(prob, FBDF(), saveat = 0.1);
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 end
diff --git a/test/pde_systems/MOLtest2.jl b/test/pde_systems/MOLtest2.jl
index e87496fdb..cef11ec27 100644
--- a/test/pde_systems/MOLtest2.jl
+++ b/test/pde_systems/MOLtest2.jl
@@ -156,7 +156,7 @@ end
 
     solc = hcat(solc1[:, :], solc2[:, 2:end])
 
-    @test sol.retcode == :Success
+    @test sol.retcode == SciMLBase.ReturnCode.Success
 end
 
 @testset "Another boundaries appearing in equations case" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 65c0cab7e..a78770afd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,13 +7,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
 # Start Test Script
 
 @time begin
-
-    if GROUP == "All" || GROUP == "Diffusion_NU"
-        @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion, Non-Uniform" begin
-            include("pde_systems/MOL_1D_Linear_Diffusion_NonUniform.jl")
-        end
-    end
-
     if GROUP == "All" || GROUP == "Components"
         @time @safetestset "MOLFiniteDifference Utils" begin
             include("utils_test.jl")
@@ -33,6 +26,19 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
     end
 
+    if GROUP == "All" || GROUP == "Convection_WENO"
+        @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin
+            include("pde_systems/MOL_1D_Linear_Convection_WENO.jl")
+        end
+    end
+
+    if GROUP == "All" || GROUP == "Diffusion_NU"
+        @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion, Non-Uniform" begin
+            include("pde_systems/MOL_1D_Linear_Diffusion_NonUniform.jl")
+        end
+    end
+
+
     if GROUP == "All" || GROUP == "Integrals"
         @time @safetestset "MOLFiniteDifference Interface: Integrals" begin
             include("pde_systems/MOL_1D_Integration.jl")
@@ -100,12 +106,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         end
     end
 
-    if GROUP == "All" || GROUP == "Convection_WENO"
-        @time @safetestset "MOLFiniteDifference Interface: Linear Convection, WENO Scheme." begin
-            include("pde_systems/MOL_1D_Linear_Convection_WENO.jl")
-        end
-    end
-
     if GROUP == "All" || GROUP == "Burgers"
         @time @safetestset "MOLFiniteDifference Interface: 2D Burger's Equation" begin
             include("pde_systems/burgers_eq.jl")
@@ -124,9 +124,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         end
     end
 
-    # if GROUP == "All" || GROUP == "Mixed_Derivatives"
-    #     @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin
-    #         include("pde_systems/MOL_Mixed_Deriv.jl")
-    #     end
-    # end
+    if GROUP == "All" || GROUP == "Mixed_Derivatives"
+        @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin
+            include("pde_systems/MOL_Mixed_Deriv.jl")
+        end
+    end
 end