diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl
index 199f051ce..472962bd1 100644
--- a/src/DiffEqOperators.jl
+++ b/src/DiffEqOperators.jl
@@ -20,8 +20,12 @@ include("basic_operators.jl")
 include("matrixfree_operators.jl")
 include("jacvec_operators.jl")
 
+### Boundary Padded Arrays
+include("boundary_padded_arrays.jl")
+
 ### Boundary Operators
 include("derivative_operators/BC_operators.jl")
+include("derivative_operators/multi_dim_bc_operators.jl")
 
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
@@ -47,6 +51,9 @@ export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
        CenteredDifference, UpwindDifference
-export RobinBC, GeneralBC
+export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
+       MultiDimDirectionalBC, ComposedMultiDimBC,
+       compose, decompose, perpsize
+
 export GhostDerivativeOperator
 end # module
diff --git a/src/boundary_padded_arrays.jl b/src/boundary_padded_arrays.jl
new file mode 100644
index 000000000..19983ab82
--- /dev/null
+++ b/src/boundary_padded_arrays.jl
@@ -0,0 +1,157 @@
+# Boundary Padded Arrays
+abstract type AbstractBoundaryPaddedArray{T, N} <: AbstractArray{T, N} end
+
+"""
+A vector type that extends a vector u with ghost points at either end
+"""
+struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}} <: AbstractBoundaryPaddedArray{T, 1}
+    l::T
+    r::T
+    u::T2
+end
+Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
+Base.size(Q::BoundaryPaddedVector) = (length(Q.u) + 2,)
+Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q)
+
+function Base.getindex(Q::BoundaryPaddedVector,i)
+    if i == 1
+        return Q.l
+    elseif i == length(Q)
+        return Q.r
+    else
+        return Q.u[i-1]
+    end
+end
+
+"""
+Higher dimensional generalization of BoundaryPaddedVector, pads an array of dimension N along the dimension D with 2 Arrays of dimension N-1, stored in lower and upper
+
+"""
+struct BoundaryPaddedArray{T, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T,N}
+    lower::B #an array of dimension M = N-1, used to extend the lower index boundary
+    upper::B #Ditto for the upper index boundary
+    u::V
+end
+
+getaxis(Q::BoundaryPaddedArray{T,D,N,M,V,B}) where {T,D,N,M,V,B} = D
+
+function Base.size(Q::BoundaryPaddedArray)
+    S = [size(Q.u)...]
+    S[getaxis(Q)] += 2
+    return Tuple(S)
+end
+
+"""
+A = compose(padded_arrays::BoundaryPaddedArray...)
+
+-------------------------------------------------------------------------------------
+
+Example:
+A = compose(Ax, Ay, Az) # 3D domain
+A = compose(Ax, Ay) # 2D Domain
+
+Composes BoundaryPaddedArrays that extend the same u for each different dimension that u has in to a ComposedBoundaryPaddedArray
+
+Ax Ay and Az can be passed in any order, as long as there is exactly one BoundaryPaddedArray that extends each dimension.
+"""
+function compose(padded_arrays::BoundaryPaddedArray...)
+    N = ndims(padded_arrays[1])
+    Ds = getaxis.(padded_arrays)
+    (length(padded_arrays) == N) || throw(ArgumentError("The padded_arrays must cover every dimension - make sure that the number of padded_arrays is equal to ndims(u)."))
+    for D in Ds
+        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple Arrays that extend along dimension $D - make sure every dimension has a unique extension"))
+    end
+    reduce((|), fill(padded_arrays[1].u, (length(padded_arrays),)) .== getfield.(padded_arrays, :u)) || throw(ArgumentError("The padded_arrays do not all extend the same u!"))
+    padded_arrays = padded_arrays[sortperm([Ds...])]
+    lower = [padded_array.lower for padded_array in padded_arrays]
+    upper = [padded_array.upper for padded_array in padded_arrays]
+
+    ComposedBoundaryPaddedArray{gettype(padded_arrays[1]),N,N-1,typeof(padded_arrays[1].u),typeof(lower[1])}(lower, upper, padded_arrays[1].u)
+end
+
+# Composed BoundaryPaddedArray
+
+struct ComposedBoundaryPaddedArray{T, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T, N}
+    lower::Vector{B}
+    upper::Vector{B}
+    u::V
+end
+
+# Aliases
+AbstractBoundaryPaddedMatrix{T} = AbstractBoundaryPaddedArray{T,2}
+AbstractBoundaryPadded3Tensor{T} = AbstractBoundaryPaddedArray{T,3}
+
+BoundaryPaddedMatrix{T, D, V, B} = BoundaryPaddedArray{T, D, 2, 1, V, B}
+BoundaryPadded3Tensor{T, D, V, B} = BoundaryPaddedArray{T, D, 3, 2, V, B}
+
+ComposedBoundaryPaddedMatrix{T,V,B} = ComposedBoundaryPaddedArray{T,2,1,V,B}
+ComposedBoundaryPadded3Tensor{T,V,B} = ComposedBoundaryPaddedArray{T,3,2,V,B}
+
+Base.size(Q::ComposedBoundaryPaddedArray) = size(Q.u).+2
+
+"""
+Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray)
+
+-------------------------------------------------------------------------------------
+
+Decomposes a ComposedBoundaryPaddedArray in to components that extend along each dimension individually
+"""
+decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A), ndims(A), ndims(A)-1, typeof(lower[1])}(A.lower[i], A.upper[i], A.u) for i in 1:ndims(A)])
+
+
+Base.length(Q::AbstractBoundaryPaddedArray) = reduce((*), size(Q))
+Base.firstindex(Q::AbstractBoundaryPaddedArray, d::Int) = 1
+Base.lastindex(Q::AbstractBoundaryPaddedArray) = length(Q)
+Base.lastindex(Q::AbstractBoundaryPaddedArray, d::Int) = size(Q)[d]
+gettype(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = T
+Base.ndims(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = N
+
+add_dim(A::AbstractArray, i) = reshape(A, size(A)...,i)
+add_dim(i) = i
+
+function Base.getindex(Q::BoundaryPaddedArray{T,D,N,M,V,B}, _inds::Vararg{Int,N}) where {T,D,N,M,V,B} #supports range and colon indexing!
+    inds = [_inds...]
+    S = size(Q)
+    dim = D
+    otherinds = inds[setdiff(1:N, dim)]
+    @assert length(S) == N
+    if inds[dim] == 1
+        return Q.lower[otherinds...]
+    elseif inds[dim] == S[dim]
+        return Q.upper[otherinds...]
+    elseif typeof(inds[dim]) <: Integer
+        inds[dim] = inds[dim] - 1
+        return Q.u[inds...]
+    elseif typeof(inds[dim]) == Colon
+        if mapreduce(x -> typeof(x) != Colon, (|), otherinds)
+            return vcat(Q.lower[otherinds...],  Q.u[inds...], Q.upper[otherinds...])
+        else
+            throw("A colon on the extended dim is as yet incompatible with additional colons")
+        end
+    elseif typeof(inds[dim]) <: AbstractArray
+        throw("Range indexing not yet supported!")
+    end
+end
+
+function Base.getindex(Q::ComposedBoundaryPaddedArray{T, N, M, V, B} , inds::Vararg{Int, N}) where {T, N, M, V, B} #as yet no support for range indexing or colon indexing
+    S = size(Q)
+    @assert reduce((&), inds .<= S)
+    for (dim, index) in enumerate(inds)
+        if index == 1
+            _inds = inds[setdiff(1:N, dim)]
+            if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds)
+                return zero(T)
+            else
+                return Q.lower[dim][(_inds.-1)...]
+            end
+        elseif index == S[dim]
+            _inds = inds[setdiff(1:N, dim)]
+            if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds)
+                return zero(T)
+            else
+                return Q.upper[dim][(_inds.-1)...]
+            end
+        end
+     end
+    return Q.u[(inds.-1)...]
+end
diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index 1923957b6..97b0949e6 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -1,152 +1,173 @@
 abstract type AbstractBC{T} <: AbstractDiffEqLinearOperator{T} end
 
+
+abstract type AtomicBC{T} <: AbstractBC{T} end
+
 """
-Robin, General, and in general Neumann and Dirichlet BCs are all affine opeartors, meaning that they take the form Qx = Qax + Qb.
+Robin, General, and in general Neumann, Dirichlet and Bridge BCs are all affine opeartors, meaning that they take the form Q*x = Qa*x + Qb.
+"""
+abstract type AffineBC{T} <: AtomicBC{T} end
+
 """
-abstract type AffineBC{T,V} <: AbstractBC{T} end
+q = PeriodicBC{T}()
 
-struct PeriodicBC{T} <: AbstractBC{T}
+Qx, Qy, ... = PeriodicBC{T}(size(u)) #When all dimensions are to be extended with a periodic boundary condition.
 
+-------------------------------------------------------------------------------------
+Creates a periodic boundary condition, where the lower index end of some u is extended with the upper index end and vice versa.
+It is not reccomended to concretize this BC type in to a BandedMatrix, since the vast majority of bands will be all 0s. SpatseMatrix concretization is reccomended.
+"""
+struct PeriodicBC{T} <: AtomicBC{T}
 end
 
 """
-  The variables in l are [αl, βl, γl], and correspond to a BC of the form al*u(0) + bl*u'(0) = cl
+  q = RobinBC(left_coefficients, right_coefficients, dx::T, approximation_order) where T # When this BC extends a dimension with a uniform step size
+
+  q = RobinBC(left_coefficients, right_coefficients, dx::Vector{T}, approximation_order) where T # When this BC extends a dimension with a non uniform step size. dx should be the vector of step sizes for the whole dimension
+
+-------------------------------------------------------------------------------------
 
+  The variables in l are [αl, βl, γl], and correspond to a BC of the form αl*u(0) + βl*u'(0) = γl imposed on the lower index boundary.
+  The variables in r are [αl, βl, γl], and correspond to an analagous boundary on the higher index end.
   Implements a robin boundary condition operator Q that acts on a vector to give an extended vector as a result
   Referring to (https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf)
-
   Write vector b̄₁ as a vertical concatanation with b0 and the rest of the elements of b̄ ₁, denoted b̄`₁, the same with ū into u0 and ū`. b̄`₁ = b̄`_2 = fill(β/Δx, length(stencil)-1)
   Pull out the product of u0 and b0 from the dot product. The stencil used to approximate u` is denoted s. b0 = α+(β/Δx)*s[1]
   Rearrange terms to find a general formula for u0:= -b̄`₁̇⋅ū`/b0 + γ/b0, which is dependent on ū` the robin coefficients and Δx.
   The non identity part of Qa is qa:= -b`₁/b0 = -β.*s[2:end]/(α+β*s[1]/Δx). The constant part is Qb = γ/(α+β*s[1]/Δx)
   do the same at the other boundary (amounts to a flip of s[2:end], with the other set of boundary coeffs)
 """
-struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T,V}
-    a_l::V
+struct RobinBC{T} <: AffineBC{T}
+    a_l::Vector{T}
     b_l::T
-    a_r::V
+    a_r::Vector{T}
     b_r::T
-    function RobinBC(l::AbstractArray{T}, r::AbstractArray{T}, dx::AbstractArray{T}, order = one(T)) where {T}
+    function RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dx::T, order = 1) where {T}
         αl, βl, γl = l
         αr, βr, γr = r
-        dx_l, dx_r = dx
 
         s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order
 
-        a_l = -s[2:end]./(αl*dx_l/βl + s[1])
-        a_r = s[end:-1:2]./(αr*dx_r/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
+        a_l = -s[2:end]./(αl*dx/βl + s[1])
+        a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
 
-        b_l = γl/(αl+βl*s[1]/dx_l)
-        b_r = γr/(αr-βr*s[1]/dx_r)
+        b_l = γl/(αl+βl*s[1]/dx)
+        b_r = γr/(αr-βr*s[1]/dx)
 
-        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
+        return new{T}(a_l, b_l, a_r, b_r)
+    end
+    function RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
+        αl, βl, γl = l
+        αr, βr, γr = r
+
+        s_index = Array(one(T):convert(T,order+1))
+        dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]
+
+        s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
+        denom_l = αl+βl*s[1]/dx_l[1]
+        denom_r = αr-βr*s[1]/dx_r[end]
+
+        a_l = -βl.*s[2:end]./(denom_l*dx_l[2:end])
+        a_r = βr.*s[end:-1:2]./(denom_r*dx_r[1:(end-1)]) # for other boundary stencil is flippedlr with *opposite sign*
+
+        b_l = γl/denom_l
+        b_r = γr/denom_r
+
+        return new{T}(a_l, b_l, a_r, b_r)
     end
 end
 
+
+
 """
+q = GeneralBC(α_leftboundary, α_rightboundary, dx::T, approximation_order)
+
+-------------------------------------------------------------------------------------
+
 Implements a generalization of the Robin boundary condition, where α is a vector of coefficients.
 Represents a condition of the form α[1] + α[2]u[0] + α[3]u'[0] + α[4]u''[0]+... = 0
 Implemented in a similar way to the RobinBC (see above).
-
 This time there are multiple stencils for multiple derivative orders - these can be written as a matrix S.
 All components that multiply u(0) are factored out, turns out to only involve the first colum of S, s̄0. The rest of S is denoted S`. the coeff of u(0) is s̄0⋅ᾱ[3:end] + α[2].
 the remaining components turn out to be ᾱ[3:end]⋅(S`ū`) or equivalantly (transpose(ᾱ[3:end])*S`)⋅ū`. Rearranging, a stencil q_a to be dotted with ū` upon extension can readily be found, along with a constant component q_b
 """
-struct GeneralBC{T, V<:AbstractVector{T}} <:AffineBC{T,V}
-    a_l::V
+struct GeneralBC{T} <:AffineBC{T}
+    a_l::Vector{T}
     b_l::T
-    a_r::V
+    a_r::Vector{T}
     b_r::T
-    function GeneralBC(αl::AbstractArray{T}, αr::AbstractArray{T}, dx::AbstractArray{T}, order = 1) where {T}
-        dx_l, dx_r = dx
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
         nl = length(αl)
         nr = length(αr)
         S_l = zeros(T, (nl-2, order+nl-2))
         S_r = zeros(T, (nr-2, order+nr-2))
 
         for i in 1:(nl-2)
-            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, nl-2-i-order))] #am unsure if the length of the dummy_x is correct here
+            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
         end
+
         for i in 1:(nr-2)
-            S_r[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, nr-2-i-order))]
+            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
         end
-        s0_l = S_l[:,1] ; Sl = S_l[2:end,:]
-        s0_r = -S_r[:,1] ; Sr = -S_r[2:end,:]
+        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
+        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
         denoml = αl[2] .+ αl[3:end] ⋅ s0_l
         denomr = αr[2] .+ αr[3:end] ⋅ s0_r
 
-        a_l = -transpose(αl) * Sl ./denoml
-        a_r = -transpose(αr) * Sr ./denomr
+        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
+        a_r = -transpose(transpose(αr[3:end]) * Sr) ./denomr
 
         b_l = -αl[1]/denoml
         b_r = -αr[1]/denomr
-        new{T, typeof(a_l)}(a_l,b_l,reverse!(a_r),b_r)
+        new{T}(a_l,b_l,reverse!(a_r),b_r)
     end
-end
 
-#implement Neumann and Dirichlet as special cases of RobinBC
-NeumannBC(α::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dx, order)
-DirichletBC(α::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where T = RobinBC([one(T), zero(T), α[1]], [one(T), zero(T), α[2]], dx, order)
-# other acceptable argument signatures
-RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], [dx_l, dx_r], order)
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
 
-# this  is 'boundary padded vector' as opposed to 'boundary padded array' to distinguish it from the n dimensional implementation that will eventually be neeeded
-struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}} <: AbstractVector{T}
-    l::T
-    r::T
-    u::T2
-end
+        nl = length(αl)
+        nr = length(αr)
+        dx_l, dx_r = (dx[1:(order+nl-2)], reverse(dx[(end-order-nr+3):end]))
+        S_l = zeros(T, (nl-2, order+nl-2))
+        S_r = zeros(T, (nr-2, order+nr-2))
 
-Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+        for i in 1:(nl-2)
+            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
+        end
 
-Base.size(Q::AbstractBC) = (Inf, Inf) #Is this nessecary?
-Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
-Base.size(Q::BoundaryPaddedVector) = (length(Q),)
-Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q)
-
-function Base.getindex(Q::BoundaryPaddedVector,i)
-    if i == 1
-        return Q.l
-    elseif i == length(Q)
-        return Q.r
-    else
-        return Q.u[i-1]
-    end
-end
+        for i in 1:(nr-2)
+            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
+        end
+        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
+        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-function LinearAlgebra.Array(Q::AffineBC{T,V}, N::Int) where {T,V}
-    Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
-    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
-    return (Array(Q_L), Q_b)
-end
+        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
+        denomr = αr[2] .+ αr[3:end] ⋅ s0_r
 
-function SparseArrays.SparseMatrixCSC(Q::AffineBC{T,V}, N::Int) where {T,V}
-    Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
-    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
-    return (Q_L, Q_b)
-end
+        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
+        a_r = -transpose(transpose(αr[3:end]) * Sr) ./denomr
 
-function SparseArrays.sparse(Q::AffineBC{T,V}, N::Int) where {T,V}
-    SparseMatrixCSC(Q,N)
+        b_l = -αl[1]/denoml
+        b_r = -αr[1]/denomr
+        new{T}(a_l,b_l,reverse!(a_r),b_r)
+    end
 end
 
-LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T = Array([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))])
-SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T = [transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))]
-SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
 
-function LinearAlgebra.Array(Q::BoundaryPaddedVector)
-    return [Q.l; Q.u; Q.r]
-end
 
-function Base.convert(::Type{Array},A::AbstractBC{T}) where T
-    Array(A)
-end
+#implement Neumann and Dirichlet as special cases of RobinBC
+NeumannBC(α::AbstractVector{T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dx, order)
+DirichletBC(αl::T, αr::T) where T = RobinBC([one(T), zero(T), αl], [one(T), zero(T), αr], 1.0, 2.0 )
+#specialized constructors for Neumann0 and Dirichlet0
+Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
+Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC([zero(T), zero(T)], dx, order)
 
-function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
-    SparseMatrixCSC(A)
-end
+# other acceptable argument signatures
+#RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
 
-function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
-    SparseMatrixCSC(A)
-end
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
+
+Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?
+
+gettype(Q::AbstractBC{T}) where T = T
diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl
index ea913f517..e1edad4e6 100644
--- a/src/derivative_operators/concretization.jl
+++ b/src/derivative_operators/concretization.jl
@@ -95,3 +95,155 @@ end
 function Base.convert(::Type{AbstractMatrix},A::DerivativeOperator{T}) where T
     BandedMatrix(A)
 end
+
+
+################################################################################
+# Boundary Padded Array concretizations
+################################################################################
+
+function LinearAlgebra.Array(Q::BoundaryPaddedArray{T,D,N,M,V,B}) where {T,D,N,M,V,B}
+    S = size(Q)
+    out = zeros(T, S...)
+    dim = D
+    ulowview = selectdim(out, dim, 1)
+    uhighview = selectdim(out, dim, S[dim])
+    uview = selectdim(out, dim, 2:(S[dim]-1))
+    ulowview .= Q.lower
+    uhighview .= Q.upper
+    uview .= Q.u
+    return out
+end
+
+function LinearAlgebra.Array(Q::ComposedBoundaryPaddedArray{T,N,M,V,B}) where {T,N,M,V,B}
+    S = size(Q)
+    out = zeros(T, S...)
+    dimset = 1:N
+    uview = out
+    for dim in dimset
+        ulowview = selectdim(out, dim, 1)
+        uhighview = selectdim(out, dim, S[dim])
+        uview = selectdim(uview, dim, 2:(S[dim]-1))
+        for (index, otherdim) in enumerate(setdiff(dimset, dim))
+            ulowview = selectdim(ulowview, index, 2:(S[otherdim]-1))
+            uhighview = selectdim(uhighview, index, 2:(S[otherdim]-1))
+        end
+        ulowview .= Q.lower[dim]
+        uhighview .= Q.upper[dim]
+    end
+    uview .= Q.u
+    return out
+end
+
+function Base.convert(::Type{AbstractArray}, A::AbstractBoundaryPaddedArray)
+    Array(A)
+end
+
+################################################################################
+# Boundary Condition Operator concretizations
+################################################################################
+
+#Atomic BCs
+function LinearAlgebra.Array(Q::AffineBC{T}, N::Int) where {T}
+    Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
+    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
+    return (Array(Q_L), Q_b)
+end
+
+function SparseArrays.SparseMatrixCSC(Q::AffineBC{T}, N::Int) where {T}
+    Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
+    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
+    return (Q_L, Q_b)
+end
+
+function BandedMatrices.BandedMatrix(Q::AffineBC{T}, N::Int) where {T}
+    Q_l = BandedMatrix{T}(Eye(N), (length(Q.a_r)-1, length(Q.a_l)-1))
+    inbands_setindex!(Q_L, Q.a_l, 1, 1:length(Q.a_l))
+    inbands_setindex!(Q_L, Q.a_r, N, (N-length(Q.a_r)+1):N)
+    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
+    return (Q_L, Q_b)
+end
+
+function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
+    SparseMatrixCSC(Q,N)
+end
+
+LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T = (Array([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))]), zeros(T, N))
+SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T = ([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))], zeros(T, N))
+SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
+function BandedMatrices.BandedMatrix(Q::PeriodicBC{T}, N::Int) where T #Not reccomended!
+    Q_array = BandedMatrix{T}(Eye(N), (N-1, N-1))
+    Q_array[1, end] = one(T)
+    Q_array[1, 1] = zero(T)
+    Q_array[end, 1] = one(T)
+    Q_array[end, end] = zero(T)
+
+    return (Q_array, zeros(T, N))
+end
+
+function LinearAlgebra.Array(Q::BoundaryPaddedVector)
+    return [Q.l; Q.u; Q.r]
+end
+
+function Base.convert(::Type{Array},A::AbstractBC{T}) where T
+    Array(A)
+end
+
+function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
+    SparseMatrixCSC(A)
+end
+
+function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
+    SparseMatrixCSC(A)
+end
+
+# Multi dimensional BC operators
+
+"""
+Returns a tuple, the first element of which is an array of the shape of the boundary,
+filled with the linear operator parts of the respective Atomic BCs.
+the second element is a simularly sized array of the affine parts.
+"""
+function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, K}, M) where {T, B, D,N,K}
+    bc_tuples = Array.(Q.BCs, fill(M, size(Q.BCs)))
+    Q_L = [bc_tuple[1] for bc_tuple in bc_tuples]
+    inds = Array(1:N)
+    inds[1], inds[D] = inds[D], inds[1]
+    Q_b = [permutedims(add_dims(bc_tuple[2], N-1), inds) for bc_tuple in bc_tuples]
+
+    return (Q_L, Q_b)
+end
+
+"""
+Returns a tuple, the first element of which is a sparse array of the shape of the boundary,
+filled with the linear operator parts of the respective Atomic BCs.
+the second element is a simularly sized array of the affine parts.
+"""
+function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, K}, M) where {T, B, D,N,K}
+    bc_tuples = sparse.(Q.BCs, fill(M, size(Q.BCs)))
+    Q_L = [bc_tuple[1] for bc_tuple in bc_tuples]
+    inds = Array(1:N)
+    inds[1], inds[D] = inds[D], inds[1]
+    Q_b = [permutedims(add_dims(bc_tuple[2], N-1), inds) for bc_tuple in bc_tuples]
+
+    return (Q_L, Q_b)
+end
+
+SparseArrays.sparse(Q::MultiDimDirectionalBC, N) = SparseMatrixCSC(Q, N)
+
+function BandedMatrices.BandedMatrix(Q::MultiDimDirectionalBC{T, B, D, N, K}, M) where {T, B, D,N,K}
+    bc_tuples = BandedMatrix.(Q.BCs, fill(M, size(Q.BCs)))
+    Q_L = [bc_tuple[1] for bc_tuple in bc_tuples]
+    inds = Array(1:N)
+    inds[1], inds[D] = inds[D], inds[1]
+    Q_b = [permutedims(add_dims(bc_tuple[2], N-1),inds) for bc_tuple in bc_tuples]
+
+    return (Q_L, Q_b)
+end
+
+"""
+Returns a Tuple of MultiDimDirectionalBC Array concretizations, one for each dimension
+"""
+LinearAlgebra.Array(Q::ComposedMultiDimBC, Ns) = Tuple(Array.(Q.BCs, Ns))
+SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC, Ns...) = Tuple(sparse.(Q.BCs, Ns))
+SparseArrays.sparse(Q::ComposedMultiDimBC, Ns) = SparseMatrixCSC(Q, Ns)
+BandedMatrices.BandedMatrix(Q::ComposedMultiDimBC, Ns) = Tuple(BandedMatrix.(Q.BCs, Ns))
diff --git a/src/derivative_operators/multi_dim_bc_operators.jl b/src/derivative_operators/multi_dim_bc_operators.jl
new file mode 100644
index 000000000..9af52c2b0
--- /dev/null
+++ b/src/derivative_operators/multi_dim_bc_operators.jl
@@ -0,0 +1,161 @@
+
+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 slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
+    @assert N != 1
+    u_temp = similar(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
+
+"""
+slicemul is the only limitation on the BCs here being used up to arbitrary dimension, an N dimensional implementation is needed.
+"""
+
+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}}
+end
+
+"""
+A multiple dimensional 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 unifrom 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
+"""
+MultiDimBC(BC::Array{B,N}, dim::Integer) where {N, B<:AtomicBC} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
+#s should be size of the domain
+MultiDimBC(BC::B, s, dim::Integer) where  {B<:AtomicBC} = 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 boundarties
+
+PeriodicBC{T}(s) where T = MultiDimBC(PeriodicBC{T}(), s)
+
+NeumannBC(α::AbstractVector{T}, dxyz, order, s) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dxyz, order, s)
+DirichletBC(αl::T, αr::T, s) where T = RobinBC([one(T), zero(T), αl], [one(T), zero(T), αr], [ones(T, si) for si in s], 2.0, s)
+
+Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
+Neumann0BC(T::Type, dxyz, order, s) = NeumannBC([zero(T), zero(T)], dxyz, order, s)
+
+RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), s[setdiff(1:length(s), dim)])) for dim in 1:length(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), s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
+
+
+perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)] #the size of A perpendicular to dim
+
+"""
+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 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...)
+    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
+
+"""
+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.BC[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
+
+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), typeof(lower)}(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)
+    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), typeof(out[1][1])}([A[1] for A in out], [A[2] for A in out], u)
+end
diff --git a/test/MultiDimBC_test.jl b/test/MultiDimBC_test.jl
new file mode 100644
index 000000000..af680be81
--- /dev/null
+++ b/test/MultiDimBC_test.jl
@@ -0,0 +1,88 @@
+using LinearAlgebra, DiffEqOperators, Random, Test
+################################################################################
+# Test 2d extension
+################################################################################
+
+#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(BCx, 1)
+Qy = MultiDimBC(BCy, 2)
+
+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(BCx, 1)
+Qy = MultiDimBC(BCy, 2)
+Qz = MultiDimBC(Dirichlet0BC(Float64), size(A), 3) #Test the other constructor
+
+Ax = Qx*A
+Ay = Qy*A
+Az = Qz*A
+
+@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:7
+    sizes = rand(4:7, N)
+    A = rand(sizes...)
+
+    Q1_N = Neumann0BC(Float64, Tuple(ones(N)), 3.0, size(A))
+
+    Q = compose(Q1_N...)
+
+    A1_N = Q1_N.*fill(A, N)
+
+    A_extended = Q*A
+    @test Array(A_extended) == Array(compose(A1_N...))
+end
diff --git a/test/bc_coeff_compositions.jl b/test/bc_coeff_compositions.jl
index 6dfbf2a57..212b3ecae 100644
--- a/test/bc_coeff_compositions.jl
+++ b/test/bc_coeff_compositions.jl
@@ -36,13 +36,13 @@ end
     al = rand()
     bl = rand()
     cl = rand()
-    dx_l = rand()
+
     ar = rand()
     br = rand()
     cr = rand()
-    dx_r = rand()
+    dx = rand()
 
-    Q = RobinBC(al, bl, cl, dx_l, ar, br, cr, dx_r)
+    Q = RobinBC([al, bl, cl], [ar, br, cr], dx)
     N = 20
     L = CenteredDifference(4,4, 1.0, N)
     L2 = CenteredDifference(2,4, 1.0, N)
@@ -112,7 +112,7 @@ end
     N = length(x)
 
     L = CenteredDifference(2, 2, dx, N)
-    Q = RobinBC(1.0, 0.0, 0.0, dx, 1.0, 0.0, 0.0, dx)
+    Q = RobinBC([1.0, 0.0, 0.0], [1.0, 0.0, 0.0], dx)
     A = L*Q
 
     analytic_L = second_derivative_stencil(N) ./ dx^2
@@ -155,7 +155,7 @@ end
     N = length(x)
 
     L = CenteredDifference(2, 2, dx, N)
-    Q = RobinBC(1.0, 0.0, 4.0, dx, 1.0, 0.0, 4.0, dx)
+    Q = RobinBC([1.0, 0.0, 4.0], [1.0, 0.0, 4.0], dx)
     A = L*Q
 
     analytic_L = second_derivative_stencil(N) ./ dx^2
@@ -203,7 +203,7 @@ end
     u = sin.(x)
 
     L = CenteredDifference(4, 4, dx, N)
-    Q = RobinBC(1.0, 0.0, sin(0.0), dx, 1.0, 0.0, sin(0.2+dx), dx)
+    Q = RobinBC([1.0, 0.0, sin(0.0)], [1.0, 0.0, sin(0.2+dx)], dx)
     A = L*Q
 
     analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
diff --git a/test/boundary_padded_array.jl b/test/boundary_padded_array.jl
new file mode 100644
index 000000000..f5c4a4940
--- /dev/null
+++ b/test/boundary_padded_array.jl
@@ -0,0 +1,64 @@
+using LinearAlgebra, DiffEqOperators, Random, Test
+################################################################################
+# Test BoundaryPaddedArray up to 5 dimensions
+################################################################################
+
+for dimensionality in 2:5
+    for dim in 1:dimensionality
+        sizes = rand(4:10, dimensionality)
+        A = rand(sizes...)
+        lower = Array(selectdim(A, dim, 1))
+        upper = Array(selectdim(A, dim, size(A)[dim]))
+
+        Apad = DiffEqOperators.BoundaryPaddedArray{Float64, dim, dimensionality, dimensionality-1, typeof(A), typeof(lower)}(lower, upper, selectdim(A, dim, 2:(size(A)[dim]-1)))
+
+        @test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix
+
+        for I in CartesianIndices(A)  #test getindex for all indicies of Apad
+            @test A[I] == Apad[I]
+        end
+    end
+end
+
+################################################################################
+# Test ComposedBoundaryPaddedMatrix
+################################################################################
+
+n = 5
+m = 7
+A = rand(n,m)
+A[1,1] = A[end,1] = A[1,end] = A[end,end] = 0.0
+
+lower = Vector[A[1,2:(end-1)], A[2:(end-1),1]]
+upper = Vector[A[end,2:(end-1)], A[2:(end-1),end]]
+
+Apad = DiffEqOperators.ComposedBoundaryPaddedMatrix{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1)])
+
+@test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix
+
+for i in 1:n, j in 1:m #test getindex for all indicies of Apad
+    @test A[i,j] == Apad[i,j]
+end
+
+################################################################################
+# Test ComposedBoundaryPaddedTensor{3}
+################################################################################
+
+n = 4
+m = 5
+o = 7
+A = rand(n,m,o)
+A[1,1,:] = A[end,1, :] = A[1,end, :] = A[end,end, :] = zeros(o)
+A[1,:,1] = A[end, :, 1] = A[1,:,end] = A[end,:,end] = zeros(m)
+A[:,1,1] = A[:,end,1] = A[:,1,end] = A[:,end,end] = zeros(n)
+
+lower = Matrix[A[1,2:(end-1), 2:(end-1)], A[2:(end-1),1, 2:(end-1)], A[2:(end-1), 2:(end-1), 1]]
+upper = Matrix[A[end, 2:(end-1), 2:(end-1)], A[2:(end-1), end, 2:(end-1)], A[2:(end-1), 2:(end-1), end]]
+
+Apad = DiffEqOperators.ComposedBoundaryPadded3Tensor{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1), 2:(end-1)])
+
+@test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix
+
+for I in CartesianIndices(A) #test getindex for all indicies of Apad
+    @test A[I] == Apad[I]
+end
diff --git a/test/generic_operator_validation.jl b/test/generic_operator_validation.jl
index 93d5a988b..23f1b48ea 100644
--- a/test/generic_operator_validation.jl
+++ b/test/generic_operator_validation.jl
@@ -1,7 +1,7 @@
 using DiffEqOperators
 
 n = 100
-x=0.0:0.01:2π
+x=0.0:0.005:2π
 xprime = x[2:(end-1)]
 dx=diff(x)
 y = exp.(π*x)
@@ -9,34 +9,22 @@ y_im = exp.(π*im*x)
 yim_ = y_im[2:(end-1)]
 y_ = y[2:(end-1)]
 
-@test_broken for dor in 1:6, aor in 1:8
+@test_broken for dor in 1:6, aor in 1:6
 
-    D1 = CenteredDifference(dor,aor,dx[1],length(x))
-    D2 = DiffEqOperators.FiniteDifference{Float64}(dor,aor,dx,length(x))
-    D = (D1,D2)
-    @test_broken convert(Array, D1) ≈ convert(Array, D2)
+    D1 = CenteredDifference(dor,aor,dx,length(x))
 
-    #take derivatives
+    #take derivative
     yprime1 = D1*y
-    yprime2 = D2*y
 
     #test result
     @test_broken yprime1 ≈ (π^dor)*y_ # test operator with known derivative of exp(kx)
-    @test_broken yprime2 ≈ (π^dor)*y_ # test operator with known derivative of exp(kx)
-
-    #test equivalance
-    @test_broken yprime1 ≈ yprime2
 
     #take derivatives
     y_imprime1 = D1*y_im
-    y_imprime2 = D2*y_im
 
     #test result
     @test_broken y_imprime1 ≈ ((pi*im)^dor)*yim_ # test operator with known derivative of exp(jkx)
-    @test_broken y_imprime2 ≈ ((pi*im)^dor)*yim_ # test operator with known derivative of exp(jkx)
 
-    #test equivalance
-    @test_broken y_imprime1 ≈ y_imprime2
 
     #TODO: implement specific tests for the left and right boundary regions, waiting until after update
 end
diff --git a/test/jacvec_operators.jl b/test/jacvec_operators.jl
index 12ac8b74c..5490cdbf6 100644
--- a/test/jacvec_operators.jl
+++ b/test/jacvec_operators.jl
@@ -1,64 +1,66 @@
-using DiffEqBase, DiffEqOperators, ForwardDiff, LinearAlgebra, Test
-const A = rand(300,300)
-f(du,u) = mul!(du,A,u)
-f(u) = A*u
-x = rand(300)
-v = rand(300)
-du = similar(x)
-
-cache1 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
-cache2 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
-@test DiffEqOperators.num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
-@test DiffEqOperators.num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
-@test DiffEqOperators.num_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
-
-@test DiffEqOperators.auto_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
-@test DiffEqOperators.auto_jacvec!(du, f, x, v, cache1, cache2) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
-@test DiffEqOperators.auto_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
-
-f(du,u,p,t) = mul!(du,A,u)
-f(u,p,t) = A*u
-L = JacVecOperator(f,x)
-@test L*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
-@test L*v ≈ DiffEqOperators.auto_jacvec(f, x, v)
-@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, x, v)
-DiffEqBase.update_coefficients!(L,v,nothing,nothing)
-@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, v, v)
-
-L = JacVecOperator(f,x,autodiff=false)
-DiffEqBase.update_coefficients!(L,x,nothing,nothing)
-@test L*x ≈ DiffEqOperators.num_jacvec(f, x, x)
-@test L*v ≈ DiffEqOperators.num_jacvec(f, x, v)
-@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, x, v) rtol=1e-6
-DiffEqBase.update_coefficients!(L,v,nothing,nothing)
-@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6
-
-L2 = JacVecOperator{Float64}(f)
-DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
-@test L2*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
-@test L2*v ≈ DiffEqOperators.auto_jacvec(f, v, v)
-@test mul!(du,L2,x) ≈ DiffEqOperators.auto_jacvec(f, x, x)
-
-L2 = JacVecOperator{Float64}(f,autodiff=false)
-DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
-@test L2*x ≈ DiffEqOperators.num_jacvec(f, x, x)
-@test L2*v ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6
-
-using OrdinaryDiffEq
-function lorenz(du,u,p,t)
- du[1] = 10.0(u[2]-u[1])
- du[2] = u[1]*(28.0-u[3]) - u[2]
- du[3] = u[1]*u[2] - (8/3)*u[3]
-end
-u0 = [1.0;0.0;0.0]
-tspan = (0.0,100.0)
-ff1 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0))
-ff2 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0,autodiff=false))
-for ff in [ff1, ff2]
-  prob = ODEProblem(ff,u0,tspan)
-  @test solve(prob,TRBDF2()).retcode == :Success
-  @test solve(prob,TRBDF2(linsolve=LinSolveGMRES())).retcode == :Success
-  @test solve(prob,Exprb32()).retcode == :Success
-  @test_broken sol = solve(prob,Rosenbrock23())
-  @test_broken sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))
-end
+using DiffEqBase, DiffEqOperators, ForwardDiff, LinearAlgebra, Test
+const A = rand(300,300)
+f(du,u) = mul!(du,A,u)
+f(u) = A*u
+x = rand(300)
+v = rand(300)
+du = similar(x)
+
+cache1 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
+cache2 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
+@test DiffEqOperators.num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
+@test DiffEqOperators.num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
+@test DiffEqOperators.num_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
+
+@test DiffEqOperators.auto_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
+@test DiffEqOperators.auto_jacvec!(du, f, x, v, cache1, cache2) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
+@test DiffEqOperators.auto_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
+
+f(du,u,p,t) = mul!(du,A,u)
+f(u,p,t) = A*u
+L = JacVecOperator(f,x)
+@test L*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
+@test L*v ≈ DiffEqOperators.auto_jacvec(f, x, v)
+@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, x, v)
+DiffEqBase.update_coefficients!(L,v,nothing,nothing)
+@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, v, v)
+
+L = JacVecOperator(f,x,autodiff=false)
+DiffEqBase.update_coefficients!(L,x,nothing,nothing)
+@test L*x ≈ DiffEqOperators.num_jacvec(f, x, x)
+@test L*v ≈ DiffEqOperators.num_jacvec(f, x, v)
+@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, x, v) rtol=1e-6
+DiffEqBase.update_coefficients!(L,v,nothing,nothing)
+@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6
+
+L2 = JacVecOperator{Float64}(f)
+DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
+@test L2*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
+@test L2*v ≈ DiffEqOperators.auto_jacvec(f, v, v)
+@test mul!(du,L2,x) ≈ DiffEqOperators.auto_jacvec(f, x, x)
+
+L2 = JacVecOperator{Float64}(f,autodiff=false)
+DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
+@test L2*x ≈ DiffEqOperators.num_jacvec(f, x, x)
+@test L2*v ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6
+
+using OrdinaryDiffEq
+function lorenz(du,u,p,t)
+ du[1] = 10.0(u[2]-u[1])
+ du[2] = u[1]*(28.0-u[3]) - u[2]
+ du[3] = u[1]*u[2] - (8/3)*u[3]
+end
+u0 = [1.0;0.0;0.0]
+tspan = (0.0,100.0)
+ff1 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0))
+ff2 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0,autodiff=false))
+
+
+for ff in [ff1, ff2]
+  prob = ODEProblem(ff,u0,tspan)
+  @test solve(prob,TRBDF2()).retcode == :Success
+  @test solve(prob,TRBDF2(linsolve=LinSolveGMRES())).retcode == :Success
+  @test solve(prob,Exprb32()).retcode == :Success
+  @test_broken sol = solve(prob,Rosenbrock23())
+  @test_broken sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))
+end
diff --git a/test/robin.jl b/test/robin.jl
index 8ad527e22..e9a62a88b 100644
--- a/test/robin.jl
+++ b/test/robin.jl
@@ -4,32 +4,32 @@ using LinearAlgebra, DiffEqOperators, Random, Test
 al = rand(5)
 bl = rand(5)
 cl = rand(5)
-dx_l = rand(5)
+dx = rand(5)
 ar = rand(5)
 br = rand(5)
 cr = rand(5)
-dx_r = rand(5)
+
 
 # Construct 5 arbitrary RobinBC operators
 for i in 1:5
-    Q = RobinBC(al[i], bl[i], cl[i], dx_l[i], ar[i], br[i], cr[i], dx_r[i])
+    Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]], dx[i])
 
     Q_L, Q_b = Array(Q,5i)
 
     #Check that Q_L is is correctly computed
     @test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
-    @test Q_L[1,:] ≈ [1 / (1-al[i]*dx_l[i]/bl[i]); zeros(5i-1)]
-    @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx_r[i]/br[i])]
+    @test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+    @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
 
     #Check that Q_b is computed correctly
-    @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx_l[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx_r[i])]
+    @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
 
     # Construct the extended operator and check that it correctly extends u to a (5i+2)
     # vector, along with encoding boundary condition information.
     u = rand(5i)
 
     Qextended = Q*u
-    CorrectQextended = [(cl[i]-(bl[i]/dx_l[i])*u[1])/(al[i]-bl[i]/dx_l[i]); u; (cr[i]+ (br[i]/dx_r[i])*u[5i])/(ar[i]+br[i]/dx_r[i])]
+    CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
     @test length(Qextended) ≈ 5i+2
 
     # Check concretization
@@ -42,4 +42,19 @@ for i in 1:5
 
 end
 
+#3rd order RobinBC, calculated with left stencil [-11/6 3 -3/2 1/3], right stencil [-1/3 3/2 -3 11/6] and [α,β,γ] = [1 6 10]
+u0 = -4/10
+uend = 125/12
+u = Vector(1.0:10.0)
+Q = RobinBC([1.0, 6.0, 10.0], [1.0, 6.0, 10.0], 1.0, 3)
+urobinextended = Q*u
+@test urobinextended.l ≈ u0
+@test urobinextended.r ≈ uend
+# General BC should be equivalent
+G = GeneralBC([-10.0, 1.0, 6.0], [-10.0, 1.0, 6.0], 1.0, 3)
+ugeneralextended = G*u
+@test ugeneralextended.l ≈ u0
+@test ugeneralextended.r ≈ uend
+
+
 #TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary
diff --git a/test/runtests.jl b/test/runtests.jl
index b2068c9ce..516e438d7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,9 +8,12 @@ import Base: isapprox
 @time @safetestset "BC and Coefficient Compositions" begin include("bc_coeff_compositions.jl") end
 @time @safetestset "Derivative Operators Interface" begin include("derivative_operators_interface.jl") end
 @time @safetestset "Validate and Compare Generic Operators" begin include("generic_operator_validation.jl") end
+@time @safetestset "Validate Boundary Padded Array Concretization" begin include("boundary_padded_array.jl") end
+@time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("MultiDimBC_test.jl") end
 #@time @safetestset "2nd order check" begin include("2nd_order_check.jl") end
 #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails
 #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end
 @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end
 @time @safetestset "Convolutions" begin include("convolutions.jl") end
 @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end
+