diff --git a/.gitignore b/.gitignore
index c4b4b9719..2db1a2f8d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,4 @@ dev/*
 # Build artifacts for creating documentation generated by the Documenter package
 docs/build/
 docs/site/
+.DS_Store
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 96397a032..464e8a9cf 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MethodOfLines"
 uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 authors = ["Alex Jones, <alex.jones@juliahub.com>"]
-version = "0.10.0"
+version = "0.10.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -31,7 +31,7 @@ Interpolations = "0.14"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
-PDEBase = "0.1.4"
+PDEBase = "0.1.6"
 PrecompileTools = "1"
 RuntimeGeneratedFunctions = "0.5.9"
 SafeTestsets = "0.0.1"
diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl
index 7a818230e..5bd8e6107 100644
--- a/src/MethodOfLines.jl
+++ b/src/MethodOfLines.jl
@@ -50,6 +50,7 @@ import Base.ndims
 # Interface
 include("interface/grid_types.jl")
 include("interface/scheme_types.jl")
+include("interface/callbacks.jl")
 include("interface/disc_strategy_types.jl")
 include("interface/MOLFiniteDifference.jl")
 
@@ -72,6 +73,7 @@ include("discretization/schemes/upwind_difference/upwind_diff_weights.jl")
 include("discretization/schemes/half_offset_weights.jl")
 include("discretization/schemes/extrapolation_weights.jl")
 include("discretization/differential_discretizer.jl")
+include("discretization/schemes/callbacks/callback_rules.jl")
 
 # System Parsing
 include("system_parsing/pde_system_transformation.jl")
@@ -104,6 +106,6 @@ include("precompile.jl")
 
 # Export
 export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace
-export UpwindScheme, WENOScheme, FunctionalScheme
+export UpwindScheme, WENOScheme, FunctionalScheme, MOLDiscCallback
 
 end
diff --git a/src/discretization/differential_discretizer.jl b/src/discretization/differential_discretizer.jl
index 5a6250c76..d23615d0e 100644
--- a/src/discretization/differential_discretizer.jl
+++ b/src/discretization/differential_discretizer.jl
@@ -8,6 +8,7 @@ struct DifferentialDiscretizer{T, D1, S} <: AbstractDifferentialDiscretizer
     interpmap::Dict{Num,DerivativeOperator}
     orders::Dict{Num,Vector{Int}}
     boundary::Dict{Num,DerivativeOperator}
+    callbacks
 end
 
 function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, discretization::MOLFiniteDifference, orders)
@@ -15,6 +16,7 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di
     bcs = pdesys.bcs isa Vector ? pdesys.bcs : [pdesys.bcs]
     approx_order = discretization.approx_order
     advection_scheme = discretization.advection_scheme
+    callbacks = discretization.callbacks
     upwind_order = advection_scheme isa UpwindScheme ? advection_scheme.order : 1
 
     differentialmap = Array{Pair{Num,DerivativeOperator},1}()
@@ -63,5 +65,5 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di
         push!(boundary, x => BoundaryInterpolatorExtrapolator(max(6, approx_order), dx))
     end
 
-    return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary))
+    return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary), callbacks)
 end
diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl
index 8d2ad1b5d..ac11911ec 100644
--- a/src/discretization/generate_finite_difference_rules.jl
+++ b/src/discretization/generate_finite_difference_rules.jl
@@ -80,6 +80,8 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace,
         integration_rules = []
     end
 
+    cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+
     integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms)))
-    return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules)
+    return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules)
 end
diff --git a/src/discretization/schemes/callbacks/callback_rules.jl b/src/discretization/schemes/callbacks/callback_rules.jl
new file mode 100644
index 000000000..cc88217e5
--- /dev/null
+++ b/src/discretization/schemes/callbacks/callback_rules.jl
@@ -0,0 +1,5 @@
+function generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+    cbs = derivweights.callbacks
+    cb_rules = [cb.sym => cb(s) for cb in cbs]
+    return cb_rules
+end
\ No newline at end of file
diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
index b8e99144b..8325c8794 100644
--- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
+++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
@@ -116,7 +116,7 @@ end
 
     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, indexmap, bcmap, 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, indexmapp, bcmap, 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, indexmap, bcmap, 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))(*(~~a, $(Differential(x))(u), ~~d))), ~~c), ~e) => /(*(~b..., ~c..., cartesian_nonlinear_laplacian(*(~a..., ~d...), Idx(II, s, u, indexmap), derivweights, s, indexmap, bcmap, depvars, x, u)), replacevals(~e, s, u, depvars, II, indexmap)) for x in ivs(u, s)]) for u in depvars], init = []))
 
diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl
index a4a7d7173..968e105e2 100644
--- a/src/interface/MOLFiniteDifference.jl
+++ b/src/interface/MOLFiniteDifference.jl
@@ -33,11 +33,12 @@ struct MOLFiniteDifference{G,D} <: AbstractEquationSystemDiscretization
     use_ODAE::Bool
     disc_strategy::D
     useIR::Bool
+    callbacks
     kwargs
 end
 
 # Constructors. If no order is specified, both upwind and centered differences will be 2nd order
-function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, kwargs...)
+function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, callbacks = [], kwargs...)
     if upwind_order !== nothing
         @warn "`upwind_order` no longer does anything, and will be removed in a future release. See the docs for the current interface."
     end
@@ -50,7 +51,7 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche
 
     dxs = dxs isa Dict ? dxs : Dict(dxs)
 
-    return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, kwargs)
+    return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, callbacks, kwargs)
 end
 
 PDEBase.get_time(disc::MOLFiniteDifference) = disc.time
diff --git a/src/interface/callbacks.jl b/src/interface/callbacks.jl
new file mode 100644
index 000000000..5c5fc427e
--- /dev/null
+++ b/src/interface/callbacks.jl
@@ -0,0 +1,13 @@
+struct MOLDiscCallback
+    f::Function
+    p
+    sym
+end
+
+function MOLDiscCallback(f, disc_ps)
+    symh = hash(f)
+    sym = unwrap(first(@parameters(Symbol("MOLDiscCallback_$symh"))))
+    return MOLDiscCallback(f, disc_ps, sym)
+end
+
+(M::MOLDiscCallback)(s) = M.f(s, M.p)
\ No newline at end of file
diff --git a/src/interface/solution/common.jl b/src/interface/solution/common.jl
index 96e16abd2..a862a0d67 100644
--- a/src/interface/solution/common.jl
+++ b/src/interface/solution/common.jl
@@ -19,13 +19,28 @@ function (sol::SciMLBase.PDESolution{T,N,S,D})(args...;
                 @assert i !== nothing "Independent variable $(arg_iv) in dependent variable $(dv) not found in the solution."
                 i
             end
-
+            if any(length.(sol.ivdomain) .== 1)
+                is = drop_singleton_inds(is, sol.ivdomain)
+            end
             sol.interp[dv](args[is]...)
         end
     end
+    if any(length.(sol.ivdomain) .== 1)
+        good_inds = get_good_inds(sol.ivdomain)
+        return sol.interp[dv](args[good_inds]...)
+    end
     return sol.interp[dv](args...)
 end
 
+function drop_singleton_inds(is, ivdomain)
+    bad_is = findall(a->length(a) == 1, ivdomain)
+    filter(a->!(a in bad_is), is)
+end
+
+function get_good_inds(ivdomain)
+    findall(a->length(a) > 1, ivdomain)
+end
+
 Base.@propagate_inbounds function Base.getindex(A::SciMLBase.PDESolution{T,N,S,D},
     sym) where {T,N,S,D<:MOLMetadata}
     iv = nothing
diff --git a/src/interface/solution/solution_utils.jl b/src/interface/solution/solution_utils.jl
index 2215caa5f..0bab30a75 100644
--- a/src/interface/solution/solution_utils.jl
+++ b/src/interface/solution/solution_utils.jl
@@ -12,10 +12,14 @@ function build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys, replaced_vars)
         kreplaced = get(replaced_vars, k, nothing)
         if all(length.(nodes) .> 1)
             interp = interpolate(nodes, umap[k], Gridded(Linear()))
-        else
-            i = findfirst(a -> a <= 1, length.(nodes))
-            @warn "Solution has length 1 in dimension $(ivs[i]). Interpolation will not be possible for variable $k. Solution return code is $(sol.retcode)."
+        elseif all(length.(nodes) .== 1)
+            @warn "Solution has length 1 in all dimensions. Interpolation will not be possible. Solution return code is $(sol.retcode)."
             interp = nothing
+        # If one or more ivs have only length 1 but not all
+        else
+            is_todrop = tuple(findall(a->a==1, length.(nodes))...)
+            is = findall(a->a > 1, length.(nodes))
+            interp = interpolate(nodes[is], dropdims(umap[k], dims=is_todrop), Gridded(Linear()))
         end
         if isnothing(kreplaced)
             [k => interp]
@@ -25,4 +29,5 @@ function build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys, replaced_vars)
     end |> Dict
 end
 
+
 sym_to_index(sym, syms) = findfirst(isequal(sym), syms)
diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl
index 8e3374952..b1239c845 100644
--- a/src/interface/solution/timeindep.jl
+++ b/src/interface/solution/timeindep.jl
@@ -11,7 +11,7 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada
     umap = mapreduce(vcat, dvs) do u
         let discu = discretespace.discvars[u]
             solu = map(CartesianIndices(discu)) do I
-                i = sym_to_index(discu[I], odesys.states)
+                i = sym_to_index(discu[I], get_states(odesys))
                 # Handle Observed
                 if i !== nothing
                     sol.u[i]
diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl
index 681b00866..ca722b151 100644
--- a/src/system_parsing/pde_system_transformation.jl
+++ b/src/system_parsing/pde_system_transformation.jl
@@ -35,6 +35,9 @@ function PDEBase.transform_pde_system!(v::PDEBase.VariableMap, boundarymap, sys:
 end
 
 function PDEBase.should_transform(pdesys::PDESystem, disc::MOLFiniteDifference, boundarymap)
+    if !disc.should_transform
+        return false
+    end
     if has_interfaces(boundarymap)
         @warn "The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature."
         return false
diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl
new file mode 100644
index 000000000..e5a6b89ea
--- /dev/null
+++ b/test/components/callbacks.jl
@@ -0,0 +1,38 @@
+using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, OrdinaryDiffEq
+
+@testset "Discrete callback" begin
+    @parameters x, t
+    @variables u(..)
+    Dt = Differential(t)
+    Dxx = Differential(x)^2
+    t_min = 0.0 
+    t_max = 2.0
+    x_min = 0.0
+    x_max = 2.0
+
+    a = 0.1*pi
+
+    eq = Dt(u(t, x)) ~ a*Dxx(u(t, x))
+
+    bcs = [u(t_min, x) ~ sinpi(x), u(t, x_min) ~ sinpi(t+x), u(t, x_max) ~ sinpi(t+x)]
+
+    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
+
+    @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])
+
+    cb = MOLDiscCallback((s, p) -> prod(p), [0.1, pi])
+    a_cb = cb.sym
+    cbeq = Dt(u(t, x)) ~ a_cb*Dxx(u(t, x))
+
+    @named cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)])
+
+    disc = MOLFiniteDifference([x => 0.1], t; approx_order=2, callbacks = [cb])
+
+    prob = discretize(pdesys, disc)
+    cbprob = discretize(cbpdesys, disc)
+
+    sol1 = solve(prob, Tsit5(), saveat=0.1)
+    sol2 = solve(cbprob, Tsit5(), saveat=0.1) 
+
+    @test sol1[u(t, x)] ≈ sol2[u(t, x)]
+end
\ No newline at end of file
diff --git a/test/pde_systems/MOL_1D_Linear_Diffusion.jl b/test/pde_systems/MOL_1D_Linear_Diffusion.jl
index f4ff722b4..09f01484b 100644
--- a/test/pde_systems/MOL_1D_Linear_Diffusion.jl
+++ b/test/pde_systems/MOL_1D_Linear_Diffusion.jl
@@ -42,7 +42,7 @@ using ModelingToolkit: Differential
 
     for disc in [discretization, discretization_edge, discretization_centered, discretization_approx_order4]
         # Convert the PDE problem into an ODE problem
-        prob = discrebranch tize(pdesys, disc)
+        prob = discretize(pdesys, disc)
 
         # Solve ODE problem      # Solve ODE problem
         sol = solve(prob, Tsit5(), saveat=0.1)
diff --git a/test/runtests.jl b/test/runtests.jl
index 7acef7ba8..263692a8f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,7 +30,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         @time @safetestset "MOLFiniteDifference Interface: Advanced Nonlinear Diffusion" begin
             include("pde_systems/nonlinear_laplacian_advanced.jl")
         end
-        end
     end
 
     if GROUP == "All" || GROUP == "Sol_Interface"
@@ -72,6 +71,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         @time @safetestset "ODEFunction" begin
             include("components/ODEFunction_test.jl")
         end
+        @time @safetestset "Discrete Callbacks" begin
+            include("components/callbacks.jl")
+        end
         #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
     end