diff --git a/Project.toml b/Project.toml
index 17303c5ad..07c2581cc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ DomainSets = "0.5"
 ForwardDiff = "0.10"
 LazyArrays = "0.17, 0.18, 0.19, 0.20, 0.21"
 LazyBandedMatrices = "0.3, 0.4, 0.5"
-ModelingToolkit = "4, 5"
+ModelingToolkit = "4 - 5.18"
 NNlib = "0.6, 0.7"
 NonlinearSolve = "0.3.7"
 RuntimeGeneratedFunctions = "0.4, 0.5"
diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl
index 0eccbecae..0465a017e 100644
--- a/src/derivative_operators/concretization.jl
+++ b/src/derivative_operators/concretization.jl
@@ -163,190 +163,190 @@ function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
 end
 
 # Multidimensional BC operators
-_concretize(Q::MultiDimDirectionalBC, M) = _concretize(Q.BCs, M)
+# _concretize(Q::MultiDimDirectionalBC, M) = _concretize(Q.BCs, M)
 
 function _concretize(Q::AbstractArray{T,N}, M) where {T,N}
     return (stencil.(Q, fill(M,size(Q))), affine.(Q))
 end
 
-function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
-    @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
-    blip = zeros(Int64, N)
-    blip[D] = 2
-    s_pad = s.+ blip # extend s in the right direction
-    Q = _concretize(Q.BCs, s[D])
-    ē = unit_indices(N)
-    QL = zeros(T, prod(s_pad), prod(s))
-    Qb = zeros(T, prod(s_pad))
-    ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
-    ranges[D] = ranges[D] .+ 1
-
-    interior = CartesianIndices(Tuple(ranges))
-    I1 = CartesianIndex(Tuple(ones(Int64, N)))
-    for I in interior
-        i = cartesian_to_linear(I, s_pad)
-        j = cartesian_to_linear(I-ē[D], s)
-        QL[i,j] = one(T)
-    end
-    ranges[D] = 1
-    lower = CartesianIndices((Tuple(ranges)))
-    ranges[D] = s_pad[D]
-    upper = CartesianIndices((Tuple(ranges)))
-    for K in CartesianIndices(upper)
-        I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
-        il = cartesian_to_linear(lower[K], s_pad)
-        iu = cartesian_to_linear(upper[K], s_pad)
-        Qb[il] = Q[2][I][1]
-        Qb[iu] = Q[2][I][2]
-        for k in 0:s[D]-1
-            j = cartesian_to_linear(lower[K] + k*ē[D], s)
-            QL[il, j] = Q[1][I][1][k+1]
-            QL[iu, j] = Q[1][I][2][k+1]
-        end
-    end
-
-    return (QL, Qb)
-end
-
-"""
-This is confusing, but it does work.
-"""
-function LinearAlgebra.Array(Q::ComposedMultiDimBC{T, B, N,M} , s::NTuple{N,G}) where {T, B, N, M, G<:Int}
-    for d in 1:N
-        @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
-    end
-    s_pad = s.+2
-    Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC
-
-    QL = zeros(T, prod(s_pad), prod(s))
-    Qb = zeros(T, prod(s_pad))
-
-    ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
-    interior = CartesianIndices(Tuple(ranges))
-
-    ē = unit_indices(N) #setup unit indices in each direction
-    I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
-    for I in interior #loop over interior
-        i = cartesian_to_linear(I, s_pad) #find the index on the padded side
-        j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
-        QL[i,j] = one(T)  #create a padded identity matrix
-    end
-    for dim in 1:N #Loop over boundaries
-        r_ = deepcopy(ranges)
-        r_[dim] = 1
-        lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
-        r_[dim] = s_pad[dim]
-        upper = CartesianIndices((Tuple(r_)))
-        for K in CartesianIndices(upper) #for every element of the boundaries
-            I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
-            il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
-            iu = cartesian_to_linear(upper[K], s_pad) # ditto
-            Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
-            Qb[iu] = Q[dim][2][I][2] #ditto with upper index
-            for k in 1:s[dim] #loop over the direction orthogonal to the boundary
-                j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
-                QL[il, j] = Q[dim][1][I][1][k]
-                QL[iu, j] = Q[dim][1][I][2][k]
-            end
-        end
-    end
-
-    return (QL, Qb)
-end
-
-"""
-See comments on the `Array` method for this type for an idea of what is going on.
-"""
-function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
-    @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
-    blip = zeros(Int64, N)
-    blip[D] = 2
-    s_pad = s.+ blip # extend s in the right direction
-    Q = _concretize(Q.BCs, s[D])
-    ē = unit_indices(N)
-    QL = spzeros(T, prod(s_pad), prod(s))
-    Qb = spzeros(T, prod(s_pad))
-    ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
-    ranges[D] = ranges[D] .+ 1
-
-    interior = CartesianIndices(Tuple(ranges))
-    I1 = CartesianIndex(Tuple(ones(Int64, N)))
-    for I in interior
-        i = cartesian_to_linear(I, s_pad)
-        j = cartesian_to_linear(I-ē[D], s)
-        QL[i,j] = one(T)
-    end
-    ranges[D] = 1
-    lower = CartesianIndices((Tuple(ranges)))
-    ranges[D] = s_pad[D]
-    upper = CartesianIndices((Tuple(ranges)))
-    for K in CartesianIndices(upper)
-        I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
-        il = cartesian_to_linear(lower[K], s_pad)
-        iu = cartesian_to_linear(upper[K], s_pad)
-        Qb[il] = Q[2][I][1]
-        Qb[iu] = Q[2][I][2]
-        for k in 0:s[D]-1
-            j = cartesian_to_linear(lower[K] + k*ē[D], s)
-            QL[il, j] = Q[1][I][1][k+1]
-            QL[iu, j] = Q[1][I][2][k+1]
-        end
-    end
-
-    return (QL, Qb)
-end
-
-
-function SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC{T, B, N,M}, s::NTuple{N,G}) where {T, B, N, M, G<:Int}
-    for d in 1:N
-        @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
-    end
-    s_pad = s.+2
-    Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC
-
-    QL = spzeros(T, prod(s_pad), prod(s))
-    Qb = spzeros(T, prod(s_pad))
-
-    ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
-    interior = CartesianIndices(Tuple(ranges))
-
-    ē = unit_indices(N) #setup unit indices in each direction
-    I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
-    for I in interior #loop over interior
-        i = cartesian_to_linear(I, s_pad) #find the index on the padded side
-        j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
-        QL[i,j] = one(T)  #create a padded identity matrix
-    end
-    for dim in 1:N #Loop over boundaries
-        r_ = deepcopy(ranges)
-        r_[dim] = 1
-        lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
-        r_[dim] = s_pad[dim]
-        upper = CartesianIndices((Tuple(r_)))
-        for K in CartesianIndices(upper) #for every element of the boundaries
-            I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
-            il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
-            iu = cartesian_to_linear(upper[K], s_pad) # ditto
-            Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
-            Qb[iu] = Q[dim][2][I][2] #ditto with upper index
-            for k in 1:s[dim] #loop over the direction orthogonal to the boundary
-                j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
-                QL[il, j] = Q[dim][1][I][1][k]
-                QL[iu, j] = Q[dim][1][I][2][k]
-            end
-        end
-    end
-
-    return (QL, Qb)
-end
-
-SparseArrays.sparse(Q::MultiDimDirectionalBC, s) = SparseMatrixCSC(Q, s)
-SparseArrays.sparse(Q::ComposedMultiDimBC, s) = SparseMatrixCSC(Q, s)
-
-
-function BandedMatrices.BandedMatrix(Q:: MultiDimensionalBC, M) where {T, B, D,N,K}
-    throw("Banded Matrix cocnretization not yet supported for MultiDimensionalBCs")
-end
+# function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
+#     @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
+#     blip = zeros(Int64, N)
+#     blip[D] = 2
+#     s_pad = s.+ blip # extend s in the right direction
+#     Q = _concretize(Q.BCs, s[D])
+#     ē = unit_indices(N)
+#     QL = zeros(T, prod(s_pad), prod(s))
+#     Qb = zeros(T, prod(s_pad))
+#     ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
+#     ranges[D] = ranges[D] .+ 1
+
+#     interior = CartesianIndices(Tuple(ranges))
+#     I1 = CartesianIndex(Tuple(ones(Int64, N)))
+#     for I in interior
+#         i = cartesian_to_linear(I, s_pad)
+#         j = cartesian_to_linear(I-ē[D], s)
+#         QL[i,j] = one(T)
+#     end
+#     ranges[D] = 1
+#     lower = CartesianIndices((Tuple(ranges)))
+#     ranges[D] = s_pad[D]
+#     upper = CartesianIndices((Tuple(ranges)))
+#     for K in CartesianIndices(upper)
+#         I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
+#         il = cartesian_to_linear(lower[K], s_pad)
+#         iu = cartesian_to_linear(upper[K], s_pad)
+#         Qb[il] = Q[2][I][1]
+#         Qb[iu] = Q[2][I][2]
+#         for k in 0:s[D]-1
+#             j = cartesian_to_linear(lower[K] + k*ē[D], s)
+#             QL[il, j] = Q[1][I][1][k+1]
+#             QL[iu, j] = Q[1][I][2][k+1]
+#         end
+#     end
+
+#     return (QL, Qb)
+# end
+
+# """
+# This is confusing, but it does work.
+# """
+# function LinearAlgebra.Array(Q::ComposedMultiDimBC{T, B, N,M} , s::NTuple{N,G}) where {T, B, N, M, G<:Int}
+#     for d in 1:N
+#         @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
+#     end
+#     s_pad = s.+2
+#     Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC
+
+#     QL = zeros(T, prod(s_pad), prod(s))
+#     Qb = zeros(T, prod(s_pad))
+
+#     ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
+#     interior = CartesianIndices(Tuple(ranges))
+
+#     ē = unit_indices(N) #setup unit indices in each direction
+#     I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
+#     for I in interior #loop over interior
+#         i = cartesian_to_linear(I, s_pad) #find the index on the padded side
+#         j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
+#         QL[i,j] = one(T)  #create a padded identity matrix
+#     end
+#     for dim in 1:N #Loop over boundaries
+#         r_ = deepcopy(ranges)
+#         r_[dim] = 1
+#         lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
+#         r_[dim] = s_pad[dim]
+#         upper = CartesianIndices((Tuple(r_)))
+#         for K in CartesianIndices(upper) #for every element of the boundaries
+#             I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
+#             il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
+#             iu = cartesian_to_linear(upper[K], s_pad) # ditto
+#             Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
+#             Qb[iu] = Q[dim][2][I][2] #ditto with upper index
+#             for k in 1:s[dim] #loop over the direction orthogonal to the boundary
+#                 j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
+#                 QL[il, j] = Q[dim][1][I][1][k]
+#                 QL[iu, j] = Q[dim][1][I][2][k]
+#             end
+#         end
+#     end
+
+#     return (QL, Qb)
+# end
+
+# """
+# See comments on the `Array` method for this type for an idea of what is going on.
+# """
+# function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
+#     @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
+#     blip = zeros(Int64, N)
+#     blip[D] = 2
+#     s_pad = s.+ blip # extend s in the right direction
+#     Q = _concretize(Q.BCs, s[D])
+#     ē = unit_indices(N)
+#     QL = spzeros(T, prod(s_pad), prod(s))
+#     Qb = spzeros(T, prod(s_pad))
+#     ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
+#     ranges[D] = ranges[D] .+ 1
+
+#     interior = CartesianIndices(Tuple(ranges))
+#     I1 = CartesianIndex(Tuple(ones(Int64, N)))
+#     for I in interior
+#         i = cartesian_to_linear(I, s_pad)
+#         j = cartesian_to_linear(I-ē[D], s)
+#         QL[i,j] = one(T)
+#     end
+#     ranges[D] = 1
+#     lower = CartesianIndices((Tuple(ranges)))
+#     ranges[D] = s_pad[D]
+#     upper = CartesianIndices((Tuple(ranges)))
+#     for K in CartesianIndices(upper)
+#         I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
+#         il = cartesian_to_linear(lower[K], s_pad)
+#         iu = cartesian_to_linear(upper[K], s_pad)
+#         Qb[il] = Q[2][I][1]
+#         Qb[iu] = Q[2][I][2]
+#         for k in 0:s[D]-1
+#             j = cartesian_to_linear(lower[K] + k*ē[D], s)
+#             QL[il, j] = Q[1][I][1][k+1]
+#             QL[iu, j] = Q[1][I][2][k+1]
+#         end
+#     end
+
+#     return (QL, Qb)
+# end
+
+
+# function SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC{T, B, N,M}, s::NTuple{N,G}) where {T, B, N, M, G<:Int}
+#     for d in 1:N
+#         @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
+#     end
+#     s_pad = s.+2
+#     Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC
+
+#     QL = spzeros(T, prod(s_pad), prod(s))
+#     Qb = spzeros(T, prod(s_pad))
+
+#     ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
+#     interior = CartesianIndices(Tuple(ranges))
+
+#     ē = unit_indices(N) #setup unit indices in each direction
+#     I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
+#     for I in interior #loop over interior
+#         i = cartesian_to_linear(I, s_pad) #find the index on the padded side
+#         j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
+#         QL[i,j] = one(T)  #create a padded identity matrix
+#     end
+#     for dim in 1:N #Loop over boundaries
+#         r_ = deepcopy(ranges)
+#         r_[dim] = 1
+#         lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
+#         r_[dim] = s_pad[dim]
+#         upper = CartesianIndices((Tuple(r_)))
+#         for K in CartesianIndices(upper) #for every element of the boundaries
+#             I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
+#             il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
+#             iu = cartesian_to_linear(upper[K], s_pad) # ditto
+#             Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
+#             Qb[iu] = Q[dim][2][I][2] #ditto with upper index
+#             for k in 1:s[dim] #loop over the direction orthogonal to the boundary
+#                 j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
+#                 QL[il, j] = Q[dim][1][I][1][k]
+#                 QL[iu, j] = Q[dim][1][I][2][k]
+#             end
+#         end
+#     end
+
+#     return (QL, Qb)
+# end
+
+# SparseArrays.sparse(Q::MultiDimDirectionalBC, s) = SparseMatrixCSC(Q, s)
+# SparseArrays.sparse(Q::ComposedMultiDimBC, s) = SparseMatrixCSC(Q, s)
+
+
+# function BandedMatrices.BandedMatrix(Q:: MultiDimensionalBC, M) where {T, B, D,N,K}
+#     throw("Banded Matrix cocnretization not yet supported for MultiDimensionalBCs")
+# end
 
 ################################################################################
 # Higher-Dimensional DerivativeOperator Concretizations
diff --git a/src/derivative_operators/multi_dim_bc_operators.jl b/src/derivative_operators/multi_dim_bc_operators.jl
index a9da48adb..3683bd8ea 100644
--- a/src/derivative_operators/multi_dim_bc_operators.jl
+++ b/src/derivative_operators/multi_dim_bc_operators.jl
@@ -1,8 +1,27 @@
 
 abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
 
+struct MultiDimDirectionalBC{T<:Real,S, D, N} <: MultiDimensionalBC{T, N}
+    stls :: S
+    l :: T
+    r :: T
+    coeff :: Number
+end
+
+struct ComposedMultiDimBC{T, S, N, M} <: MultiDimensionalBC{T, N}
+    BCs :: S
+    coeff :: Number
+end
+
+struct MultiDimBC{dim} end
 
+MultiDimBC{dim}(N::Int, approx_order::Int, dx::T; coeff::Number = 1) where {dim,T} = finite_stencil(approx_order,dx,dim,N,coeff)
+
+"""
+slice_rmul lets you multiply each vector like strip of an array `u` with a linear operator `A`, sliced along dimension `dim`
+"""
 @noinline function _slice_rmul!(u_temp::AbstractArray{T,N}, A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,N}
+    l = length(A)
     for J in post
         for I in pre
             u_temp[I, :, J] = A*u[I, :, J]
@@ -11,9 +30,6 @@ abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
     u_temp
 end
 
-"""
-slice_rmul lets you multiply each vector like strip of an array `u` with a linear operator `A`, sliced along dimension `dim`
-"""
 function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
     @assert N != 1
     u_temp = zero(u)
@@ -23,169 +39,281 @@ function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim:
     return u_temp
 end
 
-@noinline function _slice_rmul!(lower::AbstractArray, upper::AbstractArray, A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,B,N,M}
+@noinline function _slice_rmul!(lower::AbstractArray, upper::AbstractArray, A::MultiDimDirectionalBC{T, S, D, N}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,S,D,N}
     for J in post
         for I in pre
-
-            tmp = A[I,J]*u[I, :, J]
-            lower[I,J], upper[I,J] = tmp.l, tmp.r
+            tmp1 = 0
+            tmp2 = 0
+            len = Int(length(A.stls)/2)
+            for i in 1:len
+                tmp1 += A.stls[i]*u[I, i, J]
+                tmp2 += A.stls[i+len]*u[I,size(u,dim)-len+i,J]
+            end
+            tmp1 = A.coeff*(-tmp1*(A.l) + u[I,1,J])
+            tmp2 = A.coeff*(tmp2*(A.r) + u[I,size(u,dim),J])
+            lower[I,J], upper[I,J] = tmp1, tmp2
         end
     end
     return (lower, upper)
 end
 
-function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) where {T, B, N,M}
+function slice_rmul(A::MultiDimDirectionalBC{T, S, D, N}, u::AbstractArray{T,N}, dim::Int) where {T, S, D, N}
     @assert N != 1
-    @assert M == N-1
     lower = zeros(T,perpsize(u,dim))
     upper = zeros(T,perpsize(u,dim))
-
     _slice_rmul!(lower, upper, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
 
     return (lower, upper)
 end
 
-
-struct MultiDimDirectionalBC{T<:Number, B<:AtomicBC{T}, D, N, M} <: MultiDimensionalBC{T, N}
-    BCs::Array{B,M} #dimension M=N-1 array of BCs to extend dimension D
-end
-
-struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
-    BCs::Vector{Array{B, M}} # Why isn't this a vector of MultiDimBCs?
-end
-
-
 """
-A multidimensional BC, supporting arbitrary BCs at each boundary point.
-To construct an arbitrary BC, pass an Array of BCs with dimension `N-1`, if `N` is the dimensionality of your domain `u`
-with a size of `size(u)[setdiff(1:N, dim)]`, where dim is the dimension orthogonal to the boundary that you want to extend.
-
-It is also possible to call
-    Q_dim = MultiDimBC(YourBC, size(u), dim)
-to use YourBC for the whole boundary orthogonal to that dimension.
-
-Further, it is possible to call
-Qx, Qy, Qz... = MultiDimBC(YourBC, size(u))
-to use YourBC for the whole boundary for all dimensions. Valid for any number of dimensions greater than 1.
-However, this is only valid for Robin/General type BCs (including neummann/dirichlet) when the grid steps are equal in each dimension - including uniform grid case.
-
-In the case where you want to extend the same Robin/GeneralBC to the whole boundary with a non-uniform grid, please use
-    Qx, Qy, Qz... = RobinBC(l, r, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
-or
-    Qx, Qy, Qz... = GeneralBC(αl, αr, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
-
-There are also constructors for NeumannBC, DirichletBC and Dirichlet0BC. Simply replace `dx` in the call with the tuple dxyz... as above, and append `size(u)`` to the argument signature.
-The order is a required argument in this case,
-
-where dx, dy, and dz are vectors of grid steps.
-
-For Neumann0BC, please use
-    Qx, Qy, Qz... = Neumann0BC(T::Type, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
-where T is the element type of the domain to be extended.
+finite_stencil produces 2 stencils for computing derivates at the start & end of each vector along the dimension of padding.
+It internally produces MultiDimDirectionalBC structure with the produced stencils and lower, upper BV step size.
+Finally MultiDimDirectionalBC uses extrapolation for predicting BVs.
 """
-struct MultiDimBC{N} end
-
-
-MultiDimBC{dim}(BC::Array{B,N}) where {N, B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
-#s should be size of the domain
-MultiDimBC{dim}(BC::B, s) where  {B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)]))
-
-#Extra constructor to make a set of BC operators that extend an atomic BC Operator to the whole domain
-#Only valid in the uniform grid case!
-MultiDimBC(BC::B, s) where {B<:AtomicBC} = Tuple([MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
-
-# Additional constructors for cases when the BC is the same for all boundaries
-PeriodicBC{dim}(T,s) where dim = MultiDimBC{dim}(PeriodicBC(T), s)
-PeriodicBC(T,s) = MultiDimBC(PeriodicBC(T), s)
-
-NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where {T,dim} = RobinBC{dim}((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order, s)
-NeumannBC(α::NTuple{2,T}, dxyz, order, s) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dxyz, order, s)
-
-DirichletBC{dim}(αl::T, αr::T, s) where {T,dim} = RobinBC{dim}((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2.0, s)
-DirichletBC(αl::T, αr::T, s) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), fill(one(T), length(s)), 2.0, s)
-
-Dirichlet0BC{dim}(T::Type, s) where {dim} = DirichletBC{dim}(zero(T), zero(T), s)
-Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
-
-Neumann0BC{dim}(T::Type, dx, order, s) where {dim} = NeumannBC{dim}((zero(T), zero(T)), dx, order, s)
-Neumann0BC(T::Type, dxyz, order, s) = NeumannBC((zero(T), zero(T)), dxyz, order, s)
-
-RobinBC{dim}(l::NTuple{3,T}, r::NTuple{3,T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(RobinBC(l, r, dx, order), s)
-RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), perpindex(s,dim))) for dim in 1:length(s)])
-
-GeneralBC{dim}(αl::AbstractVector{T}, αr::AbstractVector{T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(GeneralBC(αl, αr, dx, order), s)
-GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order),perpindex(s,dim))) for dim in 1:length(s)])
+function finite_stencil(approx_order,dx::Union{T,AbstractArray{T}},D,N,coeff) where{T}
+    stencil_length = 1 + approx_order
+
+    if typeof(dx) <: Array
+        x = [0.0, cumsum(dx)...]
+        _upwind_coefs = convert(SVector{stencil_length, T}, calculate_weights(1, x[2], x[2:1+stencil_length]))
+        _downwind_coefs = convert(SVector{stencil_length,T}, calculate_weights(1, x[end-1], x[end-stencil_length:end-1]))
+        stencil_coefs = [_upwind_coefs ; _downwind_coefs]
+    else
+        dummy_x = 0.0 : approx_order
+        _upwind_coefs = convert(SVector{stencil_length, T}, (1/dx) * calculate_weights(1, 0.0, dummy_x))
+        _downwind_coefs = convert(SVector{stencil_length, T}, (1/dx) * calculate_weights(1, Float64(approx_order), dummy_x))
+        stencil_coefs = [_upwind_coefs ; _downwind_coefs]
+    end
+    MultiDimDirectionalBC{T,typeof(stencil_coefs), D, N}(stencil_coefs,dx[begin],dx[end],coeff)
+end
 
 
 """
 Q = compose(BCs...)
-
 -------------------------------------------------------------------------------------
-
 Example:
 Q = compose(Qx, Qy, Qz) # 3D domain
 Q = compose(Qx, Qy) # 2D Domain
-
 Creates a ComposedMultiDimBC operator, Q, that extends every boundary when applied to a `u` with a compatible size and number of dimensions.
-
 Qx Qy and Qz can be passed in any order, as long as there is exactly one BC operator that extends each dimension.
 """
 function compose(BCs::MultiDimDirectionalBC...)
-    T = gettype(BCs[1])
+    T = eltype(BCs[1].stls)
     N = ndims(BCs[1])
-    Ds = getaxis.(BCs)
-    (length(BCs) == N) || throw(ArgumentError("There must be enough BCs to cover every dimension - check that the number of MultiDimBCs == N"))
-    for D in Ds
-        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension"))
+    dims = getaxis.(BCs)
+    for D in dims
+        length(setdiff(dims, D)) == (length(dims)-1) || throw(ArgumentError("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension"))
     end
-    BCs = BCs[sortperm([Ds...])]
+    BCs = BCs[sortperm([dims...])]
 
-    ComposedMultiDimBC{T, Union{eltype.(BC.BCs for BC in BCs)...}, N,N-1}([condition.BCs for condition in BCs])
+    ComposedMultiDimBC{T, Union{eltype.(BC for BC in BCs)...}, N,N-1}(BCs,1)
 end
 
-Base.:+(BCs::MultiDimDirectionalBC...) = compose(BCs...)
-
 """
 Qx, Qy,... = decompose(Q::ComposedMultiDimBC{T,N,M})
-
 -------------------------------------------------------------------------------------
-
 Decomposes a ComposedMultiDimBC in to components that extend along each dimension individually
 """
-decompose(Q::ComposedMultiDimBC) = Tuple([MultiDimBC(Q.BCs[i], i) for i in 1:ndims(Q)])
+decompose(Q::ComposedMultiDimBC) = Tuple([Q.BCs[i] for i in 1:ndims(Q)])
 
-getaxis(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = D
-getboundarytype(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = B
+Base.:+(BCs::MultiDimDirectionalBC...) = compose(BCs...)
 
-Base.ndims(Q::MultiDimensionalBC{T,N}) where {T,N} = N
+getaxis(Q::MultiDimDirectionalBC{T, S, D, N}) where {T, S, D, N} = D
 
-Base.:*(BC::AtomicBC, u::AbstractArray) = MultiDimBC{1}(BC, size(u)) * u
+Base.ndims(Q::MultiDimensionalBC{T,N}) where {T,N} = N
 
-function Base.:*(Q::MultiDimDirectionalBC{T, B, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
-    @assert perpsize(u, D) == size(Q.BCs) "Size of the BCs array in the MultiDimBC is incorrect, needs to be $(perpsize(u,D)) to extend dimension $D, got $(size(Q.BCs))"
-    lower, upper = slice_rmul(Q.BCs, u, D)
-    return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
+function Base.:*(Q::MultiDimDirectionalBC{T, S, D, N}, u::AbstractArray{T, N}) where {T, S, D, N}
+    lower, upper = slice_rmul(Q, u, D)
+    return BoundaryPaddedArray{T, D, N, N-1, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, Q.coeff*u)
 end
 
-function Base.:*(Q::MultiDimDirectionalBC{T, PeriodicBC{T}, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
-    lower = selectdim(u, D, 1)
-    upper = selectdim(u, D, size(u,D))
-    return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
+function Base.:*(Q::ComposedMultiDimBC{T, S, N, K}, u::AbstractArray{T, N}) where {T, S, N, K}
+    Ds = getaxis.(Q.BCs)
+    for i in 1:length(Q.BCs)
+        lower,upper = slice_rmul(Q.BCs[i], u, Ds[i])
+        u = Array(BoundaryPaddedArray{T, Ds[i], N, N-1, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u))
+    end 
+    return Q.coeff*u 
 end
 
+*(c::Number, Q::MultiDimDirectionalBC{T, S, D, N}) where {T,S,D,N} = MultiDimDirectionalBC{T,S, D, N}(Q.stls,Q.l,Q.r,c*Q.coeff)
+*(c::Number, Q::ComposedMultiDimBC{T, S, N, K}) where {T,S,N,K} = ComposedMultiDimBC{T,S, N, K}(Q.BCs,c*Q.coeff)
+
+# abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
+
+
+# @noinline function _slice_rmul!(u_temp::AbstractArray{T,N}, A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,N}
+#     for J in post
+#         for I in pre
+#             u_temp[I, :, J] = A*u[I, :, J]
+#         end
+#     end
+#     u_temp
+# end
 
-function Base.:*(Q::ComposedMultiDimBC{T, B, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
-    for dim in 1:N
-        @assert perpsize(u, dim) == size(Q.BCs[dim]) "Size of the BCs array for dimension $dim in the MultiDimBC is incorrect, needs to be $(perpsize(u,dim)), got $(size(Q.BCs[dim]))"
-    end
-    out = slice_rmul.(Q.BCs, fill(u, N), 1:N)
-    lower = [A[1] for A in out]
-    upper = [A[2] for A in out]
-    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
-end
-
-function Base.:*(Q::ComposedMultiDimBC{T, PeriodicBC{T}, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
-    lower = [selectdim(u, d, 1) for d in 1:N]
-    upper = [selectdim(u, d, size(u, d)) for d in 1:N]
-    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
-end
+# """
+# slice_rmul lets you multiply each vector like strip of an array `u` with a linear operator `A`, sliced along dimension `dim`
+# """
+# function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
+#     @assert N != 1
+#     u_temp = zero(u)
+
+#     _slice_rmul!(u_temp, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
+
+#     return u_temp
+# end
+
+# @noinline function _slice_rmul!(lower::AbstractArray, upper::AbstractArray, A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,B,N,M}
+#     for J in post
+#         for I in pre
+
+#             tmp = A[I,J]*u[I, :, J]
+#             lower[I,J], upper[I,J] = tmp.l, tmp.r
+#         end
+#     end
+#     return (lower, upper)
+# end
+
+# function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) where {T, B, N,M}
+#     @assert N != 1
+#     @assert M == N-1
+#     lower = zeros(T,perpsize(u,dim))
+#     upper = zeros(T,perpsize(u,dim))
+
+#     _slice_rmul!(lower, upper, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
+
+#     return (lower, upper)
+# end
+
+
+# struct MultiDimDirectionalBC{T<:Number, B<:AtomicBC{T}, D, N, M} <: MultiDimensionalBC{T, N}
+#     BCs::Array{B,M} #dimension M=N-1 array of BCs to extend dimension D
+# end
+
+# struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
+#     BCs::Vector{Array{B, M}} # Why isn't this a vector of MultiDimBCs?
+# end
+
+
+# """
+# A multidimensional BC, supporting arbitrary BCs at each boundary point.
+# To construct an arbitrary BC, pass an Array of BCs with dimension `N-1`, if `N` is the dimensionality of your domain `u`
+# with a size of `size(u)[setdiff(1:N, dim)]`, where dim is the dimension orthogonal to the boundary that you want to extend.
+# It is also possible to call
+#     Q_dim = MultiDimBC(YourBC, size(u), dim)
+# to use YourBC for the whole boundary orthogonal to that dimension.
+# Further, it is possible to call
+# Qx, Qy, Qz... = MultiDimBC(YourBC, size(u))
+# to use YourBC for the whole boundary for all dimensions. Valid for any number of dimensions greater than 1.
+# However, this is only valid for Robin/General type BCs (including neummann/dirichlet) when the grid steps are equal in each dimension - including uniform grid case.
+# In the case where you want to extend the same Robin/GeneralBC to the whole boundary with a non-uniform grid, please use
+#     Qx, Qy, Qz... = RobinBC(l, r, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
+# or
+#     Qx, Qy, Qz... = GeneralBC(αl, αr, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
+# There are also constructors for NeumannBC, DirichletBC and Dirichlet0BC. Simply replace `dx` in the call with the tuple dxyz... as above, and append `size(u)`` to the argument signature.
+# The order is a required argument in this case,
+# where dx, dy, and dz are vectors of grid steps.
+# For Neumann0BC, please use
+#     Qx, Qy, Qz... = Neumann0BC(T::Type, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
+# where T is the element type of the domain to be extended.
+# """
+# struct MultiDimBC{N} end
+
+
+# MultiDimBC{dim}(BC::Array{B,N}) where {N, B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
+# #s should be size of the domain
+# MultiDimBC{dim}(BC::B, s) where  {B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)]))
+
+# #Extra constructor to make a set of BC operators that extend an atomic BC Operator to the whole domain
+# #Only valid in the uniform grid case!
+# MultiDimBC(BC::B, s) where {B<:AtomicBC} = Tuple([MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
+
+# # Additional constructors for cases when the BC is the same for all boundaries
+# PeriodicBC{dim}(T,s) where dim = MultiDimBC{dim}(PeriodicBC(T), s)
+# PeriodicBC(T,s) = MultiDimBC(PeriodicBC(T), s)
+
+# NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where {T,dim} = RobinBC{dim}((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order, s)
+# NeumannBC(α::NTuple{2,T}, dxyz, order, s) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dxyz, order, s)
+
+# DirichletBC{dim}(αl::T, αr::T, s) where {T,dim} = RobinBC{dim}((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2.0, s)
+# DirichletBC(αl::T, αr::T, s) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), fill(one(T), length(s)), 2.0, s)
+
+# Dirichlet0BC{dim}(T::Type, s) where {dim} = DirichletBC{dim}(zero(T), zero(T), s)
+# Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
+
+# Neumann0BC{dim}(T::Type, dx, order, s) where {dim} = NeumannBC{dim}((zero(T), zero(T)), dx, order, s)
+# Neumann0BC(T::Type, dxyz, order, s) = NeumannBC((zero(T), zero(T)), dxyz, order, s)
+
+# RobinBC{dim}(l::NTuple{3,T}, r::NTuple{3,T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(RobinBC(l, r, dx, order), s)
+# RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), perpindex(s,dim))) for dim in 1:length(s)])
+
+# GeneralBC{dim}(αl::AbstractVector{T}, αr::AbstractVector{T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(GeneralBC(αl, αr, dx, order), s)
+# GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order),perpindex(s,dim))) for dim in 1:length(s)])
+
+
+# """
+# Q = compose(BCs...)
+# -------------------------------------------------------------------------------------
+# Example:
+# Q = compose(Qx, Qy, Qz) # 3D domain
+# Q = compose(Qx, Qy) # 2D Domain
+# Creates a ComposedMultiDimBC operator, Q, that extends every boundary when applied to a `u` with a compatible size and number of dimensions.
+# Qx Qy and Qz can be passed in any order, as long as there is exactly one BC operator that extends each dimension.
+# """
+# function compose(BCs::MultiDimDirectionalBC...)
+#     T = gettype(BCs[1])
+#     N = ndims(BCs[1])
+#     Ds = getaxis.(BCs)
+#     (length(BCs) == N) || throw(ArgumentError("There must be enough BCs to cover every dimension - check that the number of MultiDimBCs == N"))
+#     for D in Ds
+#         length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension"))
+#     end
+#     BCs = BCs[sortperm([Ds...])]
+
+#     ComposedMultiDimBC{T, Union{eltype.(BC.BCs for BC in BCs)...}, N,N-1}([condition.BCs for condition in BCs])
+# end
+
+# Base.:+(BCs::MultiDimDirectionalBC...) = compose(BCs...)
+
+# """
+# Qx, Qy,... = decompose(Q::ComposedMultiDimBC{T,N,M})
+# -------------------------------------------------------------------------------------
+# Decomposes a ComposedMultiDimBC in to components that extend along each dimension individually
+# """
+# decompose(Q::ComposedMultiDimBC) = Tuple([MultiDimBC(Q.BCs[i], i) for i in 1:ndims(Q)])
+
+# getaxis(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = D
+# getboundarytype(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = B
+
+# Base.ndims(Q::MultiDimensionalBC{T,N}) where {T,N} = N
+
+# Base.:*(BC::AtomicBC, u::AbstractArray) = MultiDimBC{1}(BC, size(u)) * u
+
+# function Base.:*(Q::MultiDimDirectionalBC{T, B, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
+#     @assert perpsize(u, D) == size(Q.BCs) "Size of the BCs array in the MultiDimBC is incorrect, needs to be $(perpsize(u,D)) to extend dimension $D, got $(size(Q.BCs))"
+#     lower, upper = slice_rmul(Q.BCs, u, D)
+#     return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
+# end
+
+# function Base.:*(Q::MultiDimDirectionalBC{T, PeriodicBC{T}, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
+#     lower = selectdim(u, D, 1)
+#     upper = selectdim(u, D, size(u,D))
+#     return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
+# end
+
+
+# function Base.:*(Q::ComposedMultiDimBC{T, B, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
+#     for dim in 1:N
+#         @assert perpsize(u, dim) == size(Q.BCs[dim]) "Size of the BCs array for dimension $dim in the MultiDimBC is incorrect, needs to be $(perpsize(u,dim)), got $(size(Q.BCs[dim]))"
+#     end
+#     out = slice_rmul.(Q.BCs, fill(u, N), 1:N)
+#     lower = [A[1] for A in out]
+#     upper = [A[2] for A in out]
+#     return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
+# end
+
+# function Base.:*(Q::ComposedMultiDimBC{T, PeriodicBC{T}, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
+#     lower = [selectdim(u, d, 1) for d in 1:N]
+#     upper = [selectdim(u, d, size(u, d)) for d in 1:N]
+#     return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
diff --git a/test/DerivativeOperators/3D_laplacian.jl b/test/DerivativeOperators/3D_laplacian.jl
index e5a0e0b8b..90f9a541d 100644
--- a/test/DerivativeOperators/3D_laplacian.jl
+++ b/test/DerivativeOperators/3D_laplacian.jl
@@ -12,7 +12,10 @@ Dyy = CenteredDifference{2}(2, 4, dy, length(y))
 Dzz = CenteredDifference{3}(2, 4, dz, length(z))
 
 A = Dxx+Dyy+Dzz
-Q = compose(Dirichlet0BC(Float64, length.(s))...)
+q1 = MultiDimBC{1}(3,6,dx)
+q2 = MultiDimBC{2}(3,6,dy)
+q3 = MultiDimBC{3}(3,6,dz)
+Q = compose(q1,q2,q3)
 
 dt = dx/(sqrt(3)*3e8)
 t = 0.0:dt:10/3e8
diff --git a/test/DerivativeOperators/bc_coeff_compositions.jl b/test/DerivativeOperators/bc_coeff_compositions.jl
index ddd5a52c7..fa60d1120 100644
--- a/test/DerivativeOperators/bc_coeff_compositions.jl
+++ b/test/DerivativeOperators/bc_coeff_compositions.jl
@@ -140,13 +140,17 @@ end
     # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
 
     M = rand(N,10)
-    Qx = MultiDimBC{1}(Q, size(M))
-    Am = L*Qx
+    s = size(M)
+
+    ghost_LQM = zeros(N,10)
+    for i in 1:size(M,2)
+        mul!(view(ghost_LQM,:,i), L*Q, M[:,i])
+    end
+
     LQM = zeros(N,10)
     for i in 1:10
         mul!(view(LQM,:,i), L, Q*M[:,i])
     end
-    ghost_LQM = Am*M
     @test ghost_LQM ≈ LQM
 
     u = rand(N + 2)
@@ -365,16 +369,15 @@ end
     @test f.(x) ≈ ghost_f ≈ analytic_f
 
     # Check that left division with matrices works
-    M = [f2.(x) f2.(x)]
-    Qx = MultiDimBC{1}(Q, size(M))
-
-    Am = L*Qx
-    ghost_fM = Am \ M
-    s = size(M)
-    analytic_fM = analytic_Am \ reshape(M, prod(s))
-    @test ghost_fM ≈ reshape(analytic_fM, s)
+    # M = [f2.(x) f2.(x)]
+    # Qx = MultiDimBC{1}(Q, size(M))
+    # Am = L*Qx
+    # ghost_fM = Am \ M
+    # s = size(M)
+    # analytic_fM = analytic_Am \ reshape(M, prod(s))
+    # @test ghost_fM ≈ reshape(analytic_fM, s)
 
-    fM_temp = zeros(N,2)
+    # fM_temp = zeros(N,2)
 
     # Test \ Inhomogenous BC
     # f(x) = -x^2 + x + 4.0
@@ -408,24 +411,24 @@ end
     @test f.(x) ≈ ghost_f ≈ analytic_f
 
     # Check \ for Matrix
-    M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]
+    # M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]
 
-    s = size(M2)
-    Qx = MultiDimBC{1}(Q, size(M2))
+    # s = size(M2)
+    # Qx = MultiDimBC{1}(Q, size(M2))
 
-    analytic_QLm, analytic_Qbm = Array(Qx, s)
+    # analytic_QLm, analytic_Qbm = Array(Qx, s)
 
-    analytic_ALm = analytic_Lm*analytic_QLm
-    analytic_Abm = analytic_Lm*analytic_Qbm
+    # analytic_ALm = analytic_Lm*analytic_QLm
+    # analytic_Abm = analytic_Lm*analytic_Qbm
 
-    Am = L*Qx
-    analytic_M = analytic_ALm \ (reshape(M2 , prod(s)).- analytic_Abm)
-    ghost_M = Am \ M2
-    @test reshape(analytic_M, s) ≈ ghost_M
+    # Am = L*Qx
+    # analytic_M = analytic_ALm \ (reshape(M2 , prod(s)).- analytic_Abm)
+    # ghost_M = Am \ M2
+    # @test reshape(analytic_M, s) ≈ ghost_M
 
-    # Additionally test that A\M2 ≈ [f, 2.0(f-4.0)+4.0, 10.0(f-4.0)+4.0]
-    M = [f.(x) 2.0*(f.(x) .- 4.0).+4.0 10.0*(f.(x) .- 4.0).+4.0]
-    @test M ≈ reshape(analytic_M, s) ≈ ghost_M
+    # # Additionally test that A\M2 ≈ [f, 2.0(f-4.0)+4.0, 10.0(f-4.0)+4.0]
+    # M = [f.(x) 2.0*(f.(x) .- 4.0).+4.0 10.0*(f.(x) .- 4.0).+4.0]
+    # @test M ≈ reshape(analytic_M, s) ≈ ghost_M
 end
 
 @testset "Test Left Division L4 (fourth order)" begin
@@ -478,14 +481,14 @@ end
     _cond(A::GhostDerivativeOperator, u) = cond(sparse(A,length(u))[1] |> Array)
     @test norm(analytic_u - ghost_u) / norm(ghost_u) < _cond(A,u) * eps() #
 
-    M2 = [u 2.0*u 10.0*u]
-    s = size(M2)
-    Qx = MultiDimBC{1}(Q, size(M2))
-    Am = L*Qx
-    #Somehow the operator is singular
-    @test_broken analytic_M = analytic_Am \ (reshape(M2, prod(s)) .- repeat(analytic_Ab, 3))
-    @test_broken ghost_M = Am \ M2
-    @test_broken reshape(analytic_M, s) ≈ ghost_M
+    # M2 = [u 2.0*u 10.0*u]
+    # s = size(M2)
+    # Qx = MultiDimBC{1}(Q, size(M2))
+    # Am = L*Qx
+    # #Somehow the operator is singular
+    # @test_broken analytic_M = analytic_Am \ (reshape(M2, prod(s)) .- repeat(analytic_Ab, 3))
+    # @test_broken ghost_M = Am \ M2
+    # @test_broken reshape(analytic_M, s) ≈ ghost_M
 
 end
 
diff --git a/test/DerivativeOperators/multi_dim_bc_test.jl b/test/DerivativeOperators/multi_dim_bc_test.jl
index d67fab92e..67e6323fc 100644
--- a/test/DerivativeOperators/multi_dim_bc_test.jl
+++ b/test/DerivativeOperators/multi_dim_bc_test.jl
@@ -1,104 +1,108 @@
 using LinearAlgebra, SparseArrays, DiffEqOperators, Random, Test
-################################################################################
-# Test 2d extension
-################################################################################
 
-Random.seed!(7373)
-
-#Create Array
-n = 8
-m = 15
-A = rand(n,m)
-
-#Create atomic BC
-q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
-q2 = PeriodicBC(Float64)
-
-BCx = vcat(fill(q1, div(m,2)), fill(q2, m-div(m,2)))  #The size of BCx has to be all size components *except* for x
-BCy = vcat(fill(q1, div(n,2)), fill(q2, n-div(n,2)))
-
-
-Qx = MultiDimBC{1}(BCx)
-Qy = MultiDimBC{2}(BCy)
-
-Ax = Qx*A
-Ay = Qy*A
-
-@test size(Ax)[1] == size(A)[1]+2
-@test size(Ay)[2] == size(A)[2]+2
-
-for j in 1:m
-    @test Ax[:, j]  == Array(BCx[j]*A[:, j])
-end
-for i in 1:n
-    @test Ay[i,:] == Array(BCy[i]*A[i,:])
-end
-
-
-################################################################################
-# Test 3d extension
-################################################################################
-
-#Create Array
-n = 8
-m = 11
-o = 12
-A = rand(n,m, o)
-
-#Create atomic BC
-q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
-q2 = PeriodicBC(Float64)
-
-BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (m-div(m,2), o)))  #The size of BCx has to be all size components *except* for x
-BCy = vcat(fill(q1, (div(n,2), o)), fill(q2, (n-div(n,2), o)))
-BCz = fill(Dirichlet0BC(Float64), (n,m))
-
-Qx = MultiDimBC{1}(BCx)
-Qy = MultiDimBC{2}(BCy)
-Qz = MultiDimBC{3}(Dirichlet0BC(Float64), size(A)) #Test the other constructor
-
-Ax = Qx*A
-Ay = Qy*A
-Az = Qz*A
-
-Q = compose(Qx,Qy,Qz)
-QL, Qb = Array(Q, size(A))
-QLs, Qbs = sparse(Q, size(A))
-
-A_conc = QL*reshape(A, prod(size(A))) .+ Qb
-A_conc_sp = QLs*reshape(A,prod(size(A))) .+ Qbs
-
-#test BC concretization
-A_arr = Array(Q*A)
-@test reshape(A_arr, prod(size(A_arr))) ≈ A_conc_sp ≈ A_conc
-
-@test size(Ax)[1] == size(A)[1]+2
-@test size(Ay)[2] == size(A)[2]+2
-@test size(Az)[3] == size(A)[3]+2
-for j in 1:m, k in 1:o
-    @test Ax[:, j, k] == Array(BCx[j, k]*A[:, j, k])
-end
-for i in 1:n, k in 1:o
-    @test Ay[i, :, k] == Array(BCy[i, k]*A[i, :, k])
-end
-for i in 1:n, j in 1:m
-    @test Az[i, j, :] == Array(BCz[i, j]*A[i, j, :])
-end
-
-#test compositions to higher dimension
-for N in 2:6
-    sizes = rand(4:7, N)
-    local A = rand(sizes...)
-
-    Q1_N = RobinBC(Tuple(rand(3)), Tuple(rand(3)), fill(0.1, N), 4.0, size(A))
-
-    local Q = compose(Q1_N...)
-
-    A1_N = Q1_N.*fill(A, N)
-
-    local A_arr = Array(Q*A)
-    Q_l, Q_b = sparse(Q, size(A))
-
-    @test A_arr ≈ Array(compose(A1_N...))
-    @test A_arr ≈ reshape(Q_l*reshape(A, length(A)) .+ Q_b, size(A_arr)) #Test concretization
-end
+@testset "Test 2D & 3D MultiDimBC" begin
+	################################################################################
+	# Test 2d extension
+	################################################################################
+
+	s = x, y = (-1.9:0.1:1.9, -1.9:0.1:1.9)                 # domain of unpadded function
+	dx = dy = x[2] - x[1]
+	paraboloid(x::T, y::T) where T = 2*(x^2+y^2) - 4        # declare an elliptic paraboloid function
+	u0 = [paraboloid(X, Y) for X in x, Y in y]
+
+	q1 = MultiDimBC{1}(2,6,dx)                              # Padding along x
+	q2 = MultiDimBC{2}(2,6,dy)                              # Padding along y
+	Q = compose(q1,q2)                                      # composition ofMultiDimBC operators
+	u1 = Q*u0
+
+	s = x, y = (-2.0:0.1:2.0, -2.0:0.1:2.0)
+	u_pad = [paraboloid(X, Y) for X in x, Y in y]           # padded along x & y of 
+
+	for I in CartesianIndices(u_pad)
+		@test u1[I] ≈ u_pad[I] atol = 0.05
+	end
+
+	dx = dy = 0.1*ones(40)                                  # non-uniform grid
+	q1 = MultiDimBC{1}(2,6,dx)
+	q2 = MultiDimBC{2}(2,6,dy)
+	Q = compose(q1,q2)
+	u1 = Q*u0
+
+	for I in CartesianIndices(u_pad)
+		@test u1[I] ≈ u_pad[I] atol = 0.05
+	end
+
+
+	################################################################################
+	# Test 3d extension
+	################################################################################
+
+	#Create hyperboloid Array
+	s = x, y, z = (-1.9:0.1:1.9, -1.9:0.1:1.9, -1.9:0.1:1.9)
+	dx = dy = dz = x[2] - x[1]
+	hyperboloid(x::T, y::T, z::T) where T = 4*x^2+ 9*y^2 - z^2
+	u0 = [hyperboloid(X, Y, Z) for X in x, Y in y, Z in z]
+
+	# Create MultiDimBCs & compose
+	q1 = MultiDimBC{1}(3,6,dx)
+	q2 = MultiDimBC{2}(3,6,dy)
+	q3 = MultiDimBC{3}(3,6,dz)
+	Q = compose(q1,q2,q3)
+	u1 = Q*u0
+
+	# Padded analytical solution
+	s = x, y, z = (-2.0:0.1:2.0, -2.0:0.1:2.0, -2.0:0.1:2.0)
+	u_pad = [hyperboloid(X, Y, Z) for X in x, Y in y, Z in z]
+
+	# Testing
+	for I in CartesianIndices(u_pad)
+		@test u1[I] ≈ u_pad[I] atol = 0.2
+	end
+
+	# Test for non-uniform grids
+	dx = dy = dz = 0.1*ones(40)
+	q1 = MultiDimBC{1}(3,6,dx)
+	q2 = MultiDimBC{2}(3,6,dy)
+	q3 = MultiDimBC{3}(3,6,dz)
+	Q = compose(q1,q2,q3)
+	u1 = Q*u0
+
+	for I in CartesianIndices(u_pad)
+		@test u1[I] ≈ u_pad[I] atol = 0.2
+	end
+
+	# test BC concretization
+	# A_arr = Array(Q*A)
+	# @test reshape(A_arr, prod(size(A_arr))) ≈ A_conc_sp ≈ A_conc
+
+	# @test size(Ax)[1] == size(A)[1]+2
+	# @test size(Ay)[2] == size(A)[2]+2
+	# @test size(Az)[3] == size(A)[3]+2
+	# for j in 1:m, k in 1:o
+	#     @test Ax[:, j, k] == Array(BCx[j, k]*A[:, j, k])
+	# end
+	# for i in 1:n, k in 1:o
+	#     @test Ay[i, :, k] == Array(BCy[i, k]*A[i, :, k])
+	# end
+	# for i in 1:n, j in 1:m
+	#     @test Az[i, j, :] == Array(BCz[i, j]*A[i, j, :])
+	# end
+
+	# #test compositions to higher dimension
+	# for N in 2:6
+	#     sizes = rand(4:7, N)
+	#     local A = rand(sizes...)
+
+	#     Q1_N = RobinBC(Tuple(rand(3)), Tuple(rand(3)), fill(0.1, N), 4.0, size(A))
+
+	#     local Q = compose(Q1_N...)
+
+	#     A1_N = Q1_N.*fill(A, N)
+
+	#     local A_arr = Array(Q*A)
+	#     Q_l, Q_b = sparse(Q, size(A))
+
+	#     @test A_arr ≈ Array(compose(A1_N...))
+	#     @test A_arr ≈ reshape(Q_l*reshape(A, length(A)) .+ Q_b, size(A_arr)) #Test concretization
+	# end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 0a4aa0937..817f30aa0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -19,7 +19,7 @@ if GROUP == "All" || GROUP == "OperatorInterface"
     @time @safetestset "Validate Regular Derivative Operators" begin include("DerivativeOperators/regular_operator_validation.jl") end
     @time @safetestset "Validate and Compare Generic Operators" begin include("DerivativeOperators/generic_operator_validation.jl") end
     @time @safetestset "Validate Boundary Padded Array Concretization" begin include("DerivativeOperators/boundary_padded_array.jl") end
-    #@time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("DerivativeOperators/multi_dim_bc_test.jl") end
+    @time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("DerivativeOperators/multi_dim_bc_test.jl") end
     @time @safetestset "2nd order check" begin include("DerivativeOperators/2nd_order_check.jl") end
     @time @safetestset "Non-linear Diffusion" begin include("DerivativeOperators/Fast_Diffusion.jl") end
     @time @safetestset "KdV" begin include("DerivativeOperators/KdV.jl") end # 2-Soliton case needs implementation