diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl
index 66d1f6ef7..7087ad166 100644
--- a/src/DiffEqOperators.jl
+++ b/src/DiffEqOperators.jl
@@ -7,6 +7,7 @@ import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, isconstan
 using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices
 using LazyBandedMatrices, ModelingToolkit
 
+abstract type AbstractDiffEqAffineOperator{T} end
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} end
@@ -22,7 +23,7 @@ include("utils.jl")
 include("boundary_padded_arrays.jl")
 
 ### Boundary Operators
-include("derivative_operators/BC_operators.jl")
+include("derivative_operators/bc_operators.jl")
 include("derivative_operators/multi_dim_bc_operators.jl")
 
 ### Derivative Operators
@@ -30,7 +31,6 @@ include("derivative_operators/fornberg.jl")
 include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/convolutions.jl")
-include("derivative_operators/concretization.jl")
 include("derivative_operators/ghost_derivative_operator.jl")
 include("derivative_operators/derivative_operator_functions.jl")
 include("derivative_operators/coefficient_functions.jl")
@@ -42,8 +42,11 @@ include("MOL_discretization.jl")
 
 include("docstrings.jl")
 
+### Concretizations
+include("derivative_operators/concretization.jl")
+
 # The (u,p,t) and (du,u,p,t) interface
-for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition]
+for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition, GhostDerivativeOperator]
   (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
   (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
 end
@@ -53,8 +56,9 @@ export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
        CenteredDifference, UpwindDifference
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
-       MultiDimDirectionalBC, ComposedMultiDimBC,
-       compose, decompose, perpsize
+       MultiDimDirectionalBC, ComposedMultiDimBC
+export compose, decompose, perpsize
+
 export GhostDerivativeOperator
 export MOLFiniteDifference
 end # module
diff --git a/src/boundary_padded_arrays.jl b/src/boundary_padded_arrays.jl
index 47bf5f2bc..7bf9e8877 100644
--- a/src/boundary_padded_arrays.jl
+++ b/src/boundary_padded_arrays.jl
@@ -18,6 +18,8 @@ end
 
 Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
 Base.size(Q::BoundaryPaddedVector) = (length(Q.u) + 2,)
+unpadded_size(Q::BoundaryPaddedVector) = size(Q.u)
+unpadded_size(Q::AbstractArray) = size(Q)
 Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q)
 
 function Base.getindex(Q::BoundaryPaddedVector,i)
@@ -47,6 +49,7 @@ function Base.size(Q::BoundaryPaddedArray)
     S[getaxis(Q)] += 2
     return Tuple(S)
 end
+unpadded_size(Q::BoundaryPaddedArray) = size(Q.u)
 
 """
 A = compose(padded_arrays::BoundaryPaddedArray...)
@@ -95,6 +98,7 @@ 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
+unpadded_size(Q::ComposedBoundaryPaddedArray) = size(Q.u)
 
 """
 Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray)
diff --git a/src/composite_operators.jl b/src/composite_operators.jl
index 10a78ca48..67b84a474 100644
--- a/src/composite_operators.jl
+++ b/src/composite_operators.jl
@@ -17,6 +17,7 @@ struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: Abst
   op::OpType
 end
 *(α::DiffEqScalar{T,F}, L::AbstractDiffEqLinearOperator{T}) where {T,F} = DiffEqScaledOperator(α, L)
+*(α::Number, L::AbstractDiffEqLinearOperator{T}) where T = DiffEqScaledOperator(DiffEqScalar(convert(T,α)), L)
 -(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
 getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
 Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
@@ -27,15 +28,19 @@ opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
 getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
 getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} =
   L.coeff * L.op[I...]
-*(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op * x)
-*(x::AbstractVecOrMat, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
-/(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op / x)
-/(x::AbstractVecOrMat, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
-\(L::DiffEqScaledOperator, x::AbstractVecOrMat) = 1/L.coeff * (L.op \ x)
-\(x::AbstractVecOrMat, L::DiffEqScaledOperator) = L.coeff * (x \ L)
-mul!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
-  lmul!(L.coeff, mul!(Y, L.op, B))
-ldiv!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
+*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
+*(x::AbstractArray, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
+/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
+/(x::AbstractArray, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
+\(L::DiffEqScaledOperator, x::AbstractArray) = 1/L.coeff * (L.op \ x)
+\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
+for N in (2,3)
+  @eval begin
+    mul!(Y::AbstractArray{T,$N}, L::DiffEqScaledOperator{T}, B::AbstractArray{T,$N}) where {T} =
+        lmul!(Y, L.coeff, mul!(Y, L.op, B))
+  end
+end
+ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray) =
   lmul!(1/L.coeff, ldiv!(Y, L.op, B))
 factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
 for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
@@ -46,18 +51,19 @@ end
 
 # Linear Combination
 struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{T}}},
-  C<:AbstractVector{T}} <: AbstractDiffEqCompositeOperator{T}
-  ops::O
-  cache::C
-  function DiffEqOperatorCombination(ops; cache=nothing)
-    T = eltype(ops[1])
-    if cache == nothing
-      cache = Vector{T}(undef, size(ops[1], 1))
-      fill!(cache,0)
+    C<:AbstractVector{T}} <: AbstractDiffEqCompositeOperator{T}
+    ops::O
+    cache::C
+    function DiffEqOperatorCombination(ops; cache=nothing)
+        T = eltype(ops[1])
+        for i in 2:length(ops)
+            @assert size(ops[i]) == size(ops[1]) "Operators must be of the same size to be combined! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively"
+        end
+        if cache == nothing
+            cache = zeros(T, size(ops[1], 1))
+        end
+        new{T,typeof(ops),typeof(cache)}(ops, cache)
     end
-    # TODO: safecheck dimensions
-    new{T,typeof(ops),typeof(cache)}(ops, cache)
-  end
 end
 +(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops)
 +(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2))
@@ -72,10 +78,10 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) =
 size(L::DiffEqOperatorCombination, args...) = size(L.ops[1], args...)
 getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops)
 getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N} = sum(op -> op[I...], L.ops)
-*(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op * x, L.ops)
-*(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
-/(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op / x, L.ops)
-\(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
+*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
+*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
+/(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops)
+\(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
 function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector)
   mul!(y, L.ops[1], b)
   for op in L.ops[2:end]
@@ -92,7 +98,10 @@ struct DiffEqOperatorComposition{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{
   caches::C
   function DiffEqOperatorComposition(ops; caches=nothing)
     T = eltype(ops[1])
-    # TODO: safecheck dimensions
+    for i in 2:length(ops)
+      @assert size(ops[i-1], 1) == size(ops[i], 2) "Operations do not have compatable sizes! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively."
+    end
+
     if caches == nothing
       # Construct a list of caches to be used by mul! and ldiv!
       caches = []
@@ -122,12 +131,12 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorComposition) =
 size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2))
 size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m]
 opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops)
-*(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op*acc, L.ops; init=x)
-*(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
-/(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op/acc, L.ops; init=x)
-/(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
-\(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)
-\(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc\op, reverse(L.ops); init=x)
+*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x)
+*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
+/(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op/acc, L.ops; init=x)
+/(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
+\(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)
+\(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc\op, reverse(L.ops); init=x)
 function mul!(y::AbstractVector, L::DiffEqOperatorComposition, b::AbstractVector)
   mul!(L.caches[1], L.ops[1], b)
   for i in 2:length(L.ops) - 1
diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/bc_operators.jl
similarity index 71%
rename from src/derivative_operators/BC_operators.jl
rename to src/derivative_operators/bc_operators.jl
index 0f05711e0..c6e09cef8 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/bc_operators.jl
@@ -1,4 +1,4 @@
-abstract type AbstractBC{T} <: AbstractDiffEqLinearOperator{T} end
+abstract type AbstractBC{T} <: AbstractDiffEqAffineOperator{T} end
 
 
 abstract type AtomicBC{T} <: AbstractBC{T} end
@@ -14,9 +14,32 @@ struct NeumannBC{N} end
 struct Neumann0BC{N} end
 struct DirichletBC{N} end
 struct Dirichlet0BC{N} end
+
+"""
+q = PeriodicBC{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}
     PeriodicBC(T::Type) = new{T}()
 end
+
+"""
+  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}
     a_l::V
     b_l::T
@@ -57,8 +80,11 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
     end
 end
 
+stencil(q::AffineBC{T}, N::Int) where T = ([transpose(q.a_l) transpose(zeros(T, N-length(q.a_l)))], [transpose(zeros(T, N-length(q.a_r))) transpose(q.a_r)])
+affine(q::AffineBC) = (q.b_l, q.b_r)
 
-
+stencil(q::PeriodicBC{T}, N::Int) where T= ([transpose(zeros(T, N-1)) one(T)], [one(T) transpose(zeros(T, N-1))])
+affine(q::PeriodicBC{T}) where T = (zero(T), zero(T))
 """
 q = GeneralBC(α_leftboundary, α_rightboundary, dx::T, approximation_order)
 
diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl
index 30adbbcd3..ee538c6bd 100644
--- a/src/derivative_operators/concretization.jl
+++ b/src/derivative_operators/concretization.jl
@@ -163,57 +163,194 @@ function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
 end
 
 # Multi dimensional BC operators
+_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
 
 """
-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 similarly sized array of the affine parts.
+This is confusing, but it does work
 """
-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]
+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
 
-    return (Q_L, Q_b)
+    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))) #setup 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 anmd 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
 
 """
-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 similarly sized array of the affine parts.
+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, 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]
+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 (Q_L, Q_b)
+    return (QL, Qb)
 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]
+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
 
-    return (Q_L, Q_b)
+    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))) #setup 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 anmd 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
 
-"""
-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))
+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
+################################################################################
 # HIgher Dimensional Concretizations. The following concretizations return two dimensional arrays
 # which operate on flattened vectors. Mshape is the size of the unflattened array on which A is operating on.
 
@@ -564,3 +701,52 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
     end
     return L
 end
+
+# GhostDerivativeOperator Concretizations
+################################################################################
+function LinearAlgebra.Array(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
+    Q_l, Q_b = Array(A.Q,N)
+    return (Array(A.L,N)*Q_l, Array(A.L,N)*Q_b)
+end
+
+function LinearAlgebra.Array(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F,N,I<:Int}
+    Q_l, Q_b = Array(A.Q,s)
+    return (Array(A.L,s)*Q_l, Array(A.L,s)*Q_b)
+end
+
+function BandedMatrices.BandedMatrix(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
+    Q_l, Q_b = BandedMatrix(A.Q,N)
+    return (BandedMatrix(A.L,N)*Q_l, BandedMatrix(A.L,N)*Q_b)
+end
+
+function BandedMatrices.BandedMatrix(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F, N, I<:Int}
+    Q_l,Q_b = BandedMatrix(A.Q,s)
+    return (BandedMatrix(A.L,s)*Q_l, BandedMatrix(A.L,s)*Q_b)
+end
+
+function SparseArrays.SparseMatrixCSC(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
+    Q_l, Q_b = SparseMatrixCSC(A.Q,N)
+    return (SparseMatrixCSC(A.L,N)*Q_l, SparseMatrixCSC(A.L,N)*Q_b)
+end
+
+function SparseArrays.SparseMatrixCSC(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F,N,I<:Int}
+    Q_l, Q_b = SparseMatrixCSC(A.Q,s)
+    return (SparseMatrixCSC(A.L,s)*Q_l, SparseMatrixCSC(A.L,s)*Q_b)
+end
+
+
+function SparseArrays.sparse(A::GhostDerivativeOperator{T, E, F}, N::Int=A.L.len) where {T,E,F}
+    return SparseMatrixCSC(A,N)
+end
+
+
+function SparseArrays.sparse(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F, N,I}
+    return SparseMatrixCSC(A,s)
+end
+
+################################################################################
+# Composite Opeartor Concretizations
+################################################################################
+Array(L::DiffEqScaledOperator, s) = L.coeff * Array(L.op, s)
+Array(L::DiffEqOperatorCombination, s) = sum(Array.(L.ops, fill(s, length(L.ops))))
+Array(L::DiffEqOperatorComposition, s) = prod(Array.(reverse(L.ops), fill(s, length(L.ops))))
diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl
index fe8f17808..55495ec80 100644
--- a/src/derivative_operators/convolutions.jl
+++ b/src/derivative_operators/convolutions.jl
@@ -306,7 +306,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector
         cur_stencil = stencil[i]
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[bc_start + i] : coeff isa Number ? coeff : true
         xtempi = cur_coeff*cur_stencil[end]*_x.r
-        @inbounds for idx in A.stencil_length:-1:1
+        @inbounds for idx in (A.boundary_stencil_length-1):-1:1
             xtempi += cur_coeff * cur_stencil[end-idx] * _x.u[end-idx+1]
         end
         x_temp[bc_start + i] = xtempi + !overwrite*x_temp[bc_start + i]
diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl
index 4d804ed4f..6f86a41ad 100644
--- a/src/derivative_operators/derivative_operator.jl
+++ b/src/derivative_operators/derivative_operator.jl
@@ -53,7 +53,7 @@ function CenteredDifference{N}(derivative_order::Int,
 
     high_boundary_coefs      = convert(SVector{boundary_point_count},reverse(map(reverse, _low_boundary_coefs*(-1)^derivative_order)))
 
-    coefficients            = coeff_func isa Nothing ? nothing : Vector{T}(undef,len)
+    coefficients            = coeff_func isa Nothing ? nothing : fill!(Vector{T}(undef,len),0)
 
     DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
         typeof(low_boundary_coefs),typeof(coefficients),
@@ -239,3 +239,6 @@ diff_axis(A::DerivativeOperator{T,N}) where {T,N} = N
 function ==(A1::DerivativeOperator, A2::DerivativeOperator)
     return all([eval(:($A1.$name == $A2.$name)) for name in fieldnames(DerivativeOperator)])
 end
+function Laplacian(aor::Int, dxyz::Union{NTuple{N, T}, NTuple{N,AbstractVector{T}}}, s::NTuple{N,I}, coeff_func=nothing) where {T,N,I<:Int}
+    return sum(CenteredDifference{i}(2, aor, dxyz[i], s[i], coeff_func) for i in 1:N)
+end
diff --git a/src/derivative_operators/derivative_operator_functions.jl b/src/derivative_operators/derivative_operator_functions.jl
index a50f767fc..bf88db786 100644
--- a/src/derivative_operators/derivative_operator_functions.jl
+++ b/src/derivative_operators/derivative_operator_functions.jl
@@ -17,28 +17,43 @@
 # Fallback mul! implementation for a single DerivativeOperator operating on an AbstractArray
 function LinearAlgebra.mul!(x_temp::AbstractArray{T}, A::DerivativeOperator{T,N}, M::AbstractArray{T}; overwrite = true) where {T,N}
 
-    # Check that x_temp has correct dimensions
+    # Check that x_temp has valid dimensions, allowing unnecessary padding in M
     v = zeros(ndims(x_temp))
-    v[N] = 2
-    @assert [size(x_temp)...]+v == [size(M)...]
+    v .= 2
+    @assert all(([size(x_temp)...] .== [size(M)...]) 
+        .| (([size(x_temp)...] .+ v) .== [size(M)...])
+        )
 
     # Check that axis of differentiation is in the dimensions of M and x_temp
-    ndimsM = ndims(M)
-    @assert N <= ndimsM
+    ndims_M = ndims(M)
+    @assert N <= ndims_M
+    @assert size(x_temp, N) + 2 == size(M, N) # differentiated dimension must be padded
 
-    dimsM = [axes(M)...]
     alldims = [1:ndims(M);]
     otherdims = setdiff(alldims, N)
 
     idx = Any[first(ind) for ind in axes(M)]
-    itershape = tuple(dimsM[otherdims]...)
     nidx = length(otherdims)
+
+    dims_M = [axes(M)...]
+    dims_x_temp = [axes(x_temp)...]
+    minimal_padding_indices = map(enumerate(dims_x_temp)) do (dim, val)
+        if dim == N || length(dims_x_temp[dim]) == length(dims_M[dim])
+            Colon()
+        else
+            dims_M[dim][begin+1:end-1]
+        end
+    end
+    minimally_padded_M = view(M, minimal_padding_indices...)
+
+    itershape = tuple(dims_x_temp[otherdims]...)
     indices = Iterators.drop(CartesianIndices(itershape), 0)
 
     setindex!(idx, :, N)
     for I in indices
+        # replace all elements of idx with corresponding elt of I, except at index N
         Base.replace_tuples!(nidx, idx, idx, otherdims, I)
-        mul!(view(x_temp, idx...), A, view(M, idx...), overwrite = overwrite)
+        mul!(view(x_temp, idx...), A, view(minimally_padded_M, idx...), overwrite = overwrite)
     end
 end
 
@@ -48,25 +63,38 @@ for MT in [2,3]
     @eval begin
         function LinearAlgebra.mul!(x_temp::AbstractArray{T,$MT}, A::DerivativeOperator{T,N,false,T2,S1,S2,T3}, M::AbstractArray{T,$MT}) where
                                                                             {T,N,T2,SL,S1<:SArray{Tuple{SL},T,1,SL},S2,T3<:Union{Nothing,Number}}
-            # Check that x_temp has correct dimensions
+            # Check that x_temp has valid dimensions, allowing unnecessary padding in M
             v = zeros(ndims(x_temp))
-            v[N] = 2
-            @assert [size(x_temp)...]+v == [size(M)...]
+            v .= 2
+            @assert all(([size(x_temp)...] .== [size(M)...]) 
+                .| (([size(x_temp)...] .+ v) .== [size(M)...])
+                )
 
             # Check that axis of differentiation is in the dimensions of M and x_temp
-            ndimsM = ndims(M)
-            @assert N <= ndimsM
+            ndims_x_temp = ndims(x_temp)
+            @assert N <= ndims_x_temp
+            @assert size(x_temp, N) + 2 == size(M, N) # differentiated dimension must be padded
 
             # Determine padding for NNlib.conv!
             bpc = A.boundary_point_count
-            pad = zeros(Int64,ndimsM)
+            pad = zeros(Int64, ndims_x_temp)
             pad[N] = bpc
 
             # Reshape x_temp for NNlib.conv!
             _x_temp = reshape(x_temp, (size(x_temp)...,1,1))
 
             # Reshape M for NNlib.conv!
-            _M = reshape(M, (size(M)...,1,1))
+            dims_M = [axes(M)...]
+            dims_x_temp = [axes(x_temp)...]
+            minimal_padding_indices = map(enumerate(dims_x_temp)) do (dim, val)
+                if dim == N || length(dims_x_temp[dim]) == length(dims_M[dim])
+                    Colon()
+                else
+                    dims_M[dim][begin+1:end-1]
+                end
+            end
+            minimally_padded_M = view(M, minimal_padding_indices...)
+            _M = reshape(minimally_padded_M, (size(minimally_padded_M)...,1,1))
 
             # Setup W, the kernel for NNlib.conv!
             s = A.stencil_coefs
@@ -84,20 +112,20 @@ for MT in [2,3]
 
             # Now deal with boundaries
             if bpc > 0
-                dimsM = [axes(M)...]
-                alldims = [1:ndims(M);]
+                dims_minimal_M = [axes(minimally_padded_M)...]
+                alldims = [1:ndims(minimally_padded_M);]
                 otherdims = setdiff(alldims, N)
 
-                idx = Any[first(ind) for ind in axes(M)]
-                itershape = tuple(dimsM[otherdims]...)
+                idx = Any[first(ind) for ind in axes(minimally_padded_M)]
+                itershape = tuple(dims_minimal_M[otherdims]...)
                 nidx = length(otherdims)
                 indices = Iterators.drop(CartesianIndices(itershape), 0)
 
                 setindex!(idx, :, N)
                 for I in indices
                     Base.replace_tuples!(nidx, idx, idx, otherdims, I)
-                    convolve_BC_left!(view(x_temp, idx...), view(M, idx...), A)
-                    convolve_BC_right!(view(x_temp, idx...), view(M, idx...), A)
+                    convolve_BC_left!(view(x_temp, idx...), view(minimally_padded_M, idx...), A)
+                    convolve_BC_right!(view(x_temp, idx...), view(minimally_padded_M, idx...), A)
                 end
             end
         end
diff --git a/src/derivative_operators/ghost_derivative_operator.jl b/src/derivative_operators/ghost_derivative_operator.jl
index 38fb060e0..0cf72b16a 100644
--- a/src/derivative_operators/ghost_derivative_operator.jl
+++ b/src/derivative_operators/ghost_derivative_operator.jl
@@ -1,9 +1,9 @@
-struct GhostDerivativeOperator{T<:Real, E<:AbstractDerivativeOperator{T}, F<:AbstractBC{T}} <: AbstractDiffEqLinearOperator{T}
+struct GhostDerivativeOperator{T, E<:AbstractDiffEqLinearOperator{T}, F<:AbstractBC{T}} <: AbstractDiffEqLinearOperator{T}
     L :: E
     Q :: F
 end
 
-function *(L::AbstractDerivativeOperator{T}, Q::AbstractBC{T}) where{T}
+function *(L::AbstractDiffEqLinearOperator{T}, Q::AbstractBC{T}) where{T}
     return GhostDerivativeOperator{T, typeof(L), typeof(Q)}(L,Q)
 end
 
@@ -11,64 +11,29 @@ function *(L::AbstractDiffEqCompositeOperator{T}, Q::AbstractBC{T}) where{T}
     return sum(map(op -> op * Q, L.ops))
 end
 
-function LinearAlgebra.mul!(x::AbstractVector, A::GhostDerivativeOperator, u::AbstractVector)
-    @assert length(u) == A.L.len == length(x)
+function *(A::GhostDerivativeOperator{T1}, u::AbstractArray{T2}) where {T1,T2}
+    #TODO Implement a function domaincheck(L::AbstractDiffEqLinearOperator, u) to see if components of L along each dimension match the size of u
+    x = zeros(promote_type(T1,T2), unpadded_size(u))
     LinearAlgebra.mul!(x, A.L, A.Q*u)
-end
-
-function *(A::GhostDerivativeOperator{T1}, u::AbstractVector{T2}) where {T1,T2}
-    @assert length(u) == A.L.len
-    x = zeros(promote_type(T1,T2), A.L.len)
-    LinearAlgebra.mul!(x, A, u)
     return x
 end
 
-function LinearAlgebra.mul!(M_temp::AbstractMatrix, A::GhostDerivativeOperator, M::AbstractMatrix)
-    @assert size(M,1) == size(M_temp,1) == A.L.len
-    for i in 1:size(M,2)
-        mul!(view(M_temp,:,i), A, view(M,:,i))
-    end
-end
-
-function *(A::GhostDerivativeOperator{T1}, M::AbstractMatrix{T2}) where {T1,T2}
-    @assert size(M,1) == A.L.len
-    M_temp = zeros(promote_type(T1,T2), A.L.len, size(M,2))
-    LinearAlgebra.mul!(M_temp, A, M)
-    return M_temp
-end
-
-
-function LinearAlgebra.ldiv!(x::AbstractVector, A::GhostDerivativeOperator, u::AbstractVector)
-    @assert length(x) == A.L.len
-    (AL,Ab) = Array(A)
-    LinearAlgebra.ldiv!(x, lu!(AL), u-Ab)
-end
-
-
-function \(A::GhostDerivativeOperator{T1}, u::AbstractVector{T2}) where {T1,T2}
-    @assert length(u) == A.L.len
-    x = zeros(promote_type(T1,T2),size(A,2))
-    LinearAlgebra.ldiv!(x, A, u)
-    return x
-end
 
-function LinearAlgebra.ldiv!(M_temp::AbstractMatrix, A::GhostDerivativeOperator, M::AbstractMatrix)
-    @assert size(M_temp) == size(M)
-    @assert A.L.len == size(M,1)
-    (AL,Ab) = Array(A)
-    LinearAlgebra.ldiv!(M_temp, lu!(AL), M .- Ab)
+function \(A::GhostDerivativeOperator, u::AbstractArray) # FIXME should have T1,T2 and promote result
+    #TODO implement check that A has compatible size with u
+    s = size(u)
+    (A_l,A_b) = sparse(A, s)
+    x = A_l \ Vector(reshape(u, length(u)) .- A_b) #Has to be converted to vector to work, A_b being sparse was causing a conversion to sparse.
+    return reshape(x, s)
 end
 
-function \(A::GhostDerivativeOperator{T1}, M::AbstractMatrix{T2}) where {T1,T2}
-    @assert A.L.len == size(M,1)
-    M_temp = zeros(promote_type(T1,T2), A.L.len, size(M,2))
-    LinearAlgebra.ldiv!(M_temp, A, M)
-    return M_temp
+function \(A::GhostDerivativeOperator, u::AbstractVector) # FIXME as above, should promote_type(T1,T2)
+    # TODO: is this specialization to u::AbstractVector really any faster?
+    @assert length(u) == size(A.L, 1)
+    (A_l,A_b) = sparse(A,length(u))
+    A_l \ Vector(u .- A_b)
 end
 
-# # Interface with other GhostDerivativeOperator and with
-# # AbstractDiffEqCompositeOperator
-
 # update coefficients
 function DiffEqBase.update_coefficients!(A::GhostDerivativeOperator,u,p,t)
     DiffEqBase.update_coefficients!(A.L,u,p,t)
@@ -93,22 +58,5 @@ Base.size(A::GhostDerivativeOperator) = (A.L.len, A.L.len)
 Base.size(A::GhostDerivativeOperator,i::Integer) = size(A)[i]
 Base.length(A::GhostDerivativeOperator) = reduce(*, size(A))
 
-
-# Concretizations, will be moved to concretizations.jl later
-function LinearAlgebra.Array(A::GhostDerivativeOperator,N::Int=A.L.len)
-    return (Array(A.L,N)*Array(A.Q,A.L.len)[1], Array(A.L,N)*Array(A.Q,A.L.len)[2])
-end
-
-function BandedMatrices.BandedMatrix(A::GhostDerivativeOperator,N::Int=A.L.len)
-    return (BandedMatrix(A.L,N)*Array(A.Q,A.L.len)[1], BandedMatrix(A.L,N)*Array(A.Q,A.L.len)[2])
-end
-
-function SparseArrays.SparseMatrixCSC(A::GhostDerivativeOperator,N::Int=A.L.len)
-    return (SparseMatrixCSC(A.L,N)*SparseMatrixCSC(A.Q,A.L.len)[1], SparseMatrixCSC(A.L,N)*SparseMatrixCSC(A.Q,A.L.len)[2])
-end
-
-function SparseArrays.sparse(A::GhostDerivativeOperator,N::Int=A.L.len)
-    return SparseMatrixCSC(A,N)
-end
-
 @inline ==(A1::GhostDerivativeOperator, A2::GhostDerivativeOperator) = A1.L == A2.L && A1.Q == A2.Q
+
diff --git a/src/derivative_operators/multi_dim_bc_operators.jl b/src/derivative_operators/multi_dim_bc_operators.jl
index 8cce46cd6..833faef39 100644
--- a/src/derivative_operators/multi_dim_bc_operators.jl
+++ b/src/derivative_operators/multi_dim_bc_operators.jl
@@ -51,7 +51,7 @@ struct MultiDimDirectionalBC{T<:Number, B<:AtomicBC{T}, D, N, M} <: MultiDimensi
 end
 
 struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
-    BCs::Vector{Array{B, M}}
+    BCs::Vector{Array{B, M}} # Why isn't this a vector of MultiDimBCs?
 end
 
 
@@ -102,7 +102,7 @@ NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where {T,dim} = RobinBC{dim}((zero
 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), 1.0, 2.0, 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)
+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)
@@ -117,7 +117,6 @@ GeneralBC{dim}(αl::AbstractVector{T}, αr::AbstractVector{T}, dx, order, s) whe
 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...)
 
@@ -153,23 +152,40 @@ 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)])
+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), typeof(lower)}(lower, upper, u)
+    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)
-    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), typeof(out[1][1])}([A[1] for A in out], [A[2] for A in out], u)
+    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
diff --git a/src/utils.jl b/src/utils.jl
index b76f31be9..54436a87a 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -13,6 +13,14 @@ function unit_indices(N::Int) #create unit CartesianIndex for each dimension
     Tuple(out)
 end
 
+function cartesian_to_linear(I::CartesianIndex, s)  #Not sure if there is a builtin that does this - convert cartesian index to linear index of an array of size s
+    out = I[1]
+    for i in 1:length(s)-1
+        out += (I[i+1]-1)*prod(s[1:i])
+    end
+    return out
+end
+
 add_dims(A::AbstractArray, n::Int; dims::Int = 1) = cat(ndims(A) + n, A, dims = dims)
 
 ""
diff --git a/test/3D_laplacian.jl b/test/3D_laplacian.jl
new file mode 100644
index 000000000..e5a0e0b8b
--- /dev/null
+++ b/test/3D_laplacian.jl
@@ -0,0 +1,37 @@
+using DiffEqOperators
+
+s = x, y, z = (-5:0.2:5, -5:0.2:5, -5:0.2:5)
+dx = dy = dz = x[2] - x[1]
+
+ricker(x::T, y::T, z::T) where T = (4*(x^2+y^2+z^2) - 6)*exp(-(x^2+y^2+z^2))
+
+u0 = [ricker(X, Y, Z) for Z in z, Y in y, X in x]
+
+Dxx = CenteredDifference{1}(2, 4, dx, length(x))
+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))...)
+
+dt = dx/(sqrt(3)*3e8)
+t = 0.0:dt:10/3e8
+
+f(u,p,t) = (3e8)^2 .*(A*Q*u)
+
+
+
+function steptime(u,uold,uolder)
+    return ((dt^2 .*f(u,0,0) .+ 5u .- 4uold .+ uolder)./2, u, uold)
+end
+let uolder = deepcopy(u0), uold = deepcopy(u0), u = deepcopy(u0)
+    u,uold,uolder = steptime(u0,uold,uolder)
+    # @gif for ti in t #4th order time stepper
+    #         u, uold, uolder = steptime(u,uold,uolder)
+    #         heatmap(u)
+    # end
+    for ti in t #4th order time stepper
+        u, uold, uolder = steptime(u,uold,uolder)
+    end
+end
+
diff --git a/test/bc_coeff_compositions.jl b/test/bc_coeff_compositions.jl
index 3524f5cdc..4d22b960b 100644
--- a/test/bc_coeff_compositions.jl
+++ b/test/bc_coeff_compositions.jl
@@ -137,12 +137,15 @@ end
     @test (A1 + B1) * u == A1 * u + B1 * u
 
     # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
+
     M = rand(N,10)
+    Qx = MultiDimBC{1}(Q, size(M))
+    Am = L*Qx
     LQM = zeros(N,10)
     for i in 1:10
         mul!(view(LQM,:,i), L, Q*M[:,i])
     end
-    ghost_LQM = A*M
+    ghost_LQM = Am*M
     @test ghost_LQM ≈ LQM
 
     u = rand(N + 2)
@@ -340,6 +343,16 @@ end
     analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
     analytic_AL = analytic_L*analytic_QL
 
+    analytic_QM = zeros((length(x)+2)*2, length(x)*2)
+    interior = CartesianIndices(Tuple([2:length(x)+1, 1:2]))
+    I1 = CartesianIndex(1,0)
+    for I in interior
+        i = DiffEqOperators.cartesian_to_linear(I, (length(x)+2, 2)) #helper function, see utils.jl
+        j = DiffEqOperators.cartesian_to_linear(I-I1, (length(x), 2))
+        analytic_QM[i,j] = 1.0
+    end
+    analytic_Am = kron(Diagonal(ones(2)), analytic_L)*analytic_QM
+
     # No affine component to the this system
     analytic_f = analytic_AL \ f2.(x)
     ghost_f = A \ f2.(x)
@@ -350,19 +363,17 @@ end
     # Additionally test that A\f2.(x) ≈ f.(x)
     @test f.(x) ≈ ghost_f ≈ analytic_f
 
-    # Check ldiv!
-    f_temp = zeros(N)
-    ldiv!(f_temp, A, f2.(x))
-    @test f_temp ≈ ghost_f ≈ analytic_f
-
     # Check that left division with matrices works
-    ghost_fM = A \ [f2.(x) f2.(x)]
-    analytic_fM = analytic_AL \ [f2.(x) f2.(x)]
-    @test ghost_fM ≈ analytic_fM
+    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)
-    ldiv!(fM_temp, A, [f2.(x) f2.(x)])
-    @test fM_temp ≈ analytic_fM
 
     # Test \ Inhomogenous BC
     # f(x) = -x^2 + x + 4.0
@@ -384,6 +395,7 @@ end
     analytic_Qb = [4.0; zeros(N); 4.0]
     analytic_AL = analytic_L*analytic_QL
     analytic_Ab = analytic_L*analytic_Qb
+    analytic_Lm = kron(Diagonal(ones(3)), analytic_L)
 
     analytic_f = analytic_AL \ (f2.(x) - analytic_Ab)
     ghost_f = A \ f2.(x)
@@ -394,25 +406,25 @@ end
     # Additionally test that A\f2.(x) ≈ f.(x)
     @test f.(x) ≈ ghost_f ≈ analytic_f
 
-    # Check ldiv!
-    f_temp = zeros(N)
-    ldiv!(f_temp, A, f2.(x))
-    @test f_temp ≈ ghost_f ≈ analytic_f
-
     # Check \ for Matrix
     M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]
-    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
-    ghost_M = A \ M2
-    @test analytic_M ≈ ghost_M
 
-    # Check ldiv! for Matrix
-    M_temp = zeros(N,3)
-    ldiv!(M_temp, A, M2)
-    @test M_temp ≈ analytic_M ≈ ghost_M
+    s = size(M2)
+    Qx = MultiDimBC{1}(Q, size(M2))
+
+    analytic_QLm, analytic_Qbm = Array(Qx, s)
+
+    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
 
     # 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 ≈ M_temp ≈ analytic_M ≈ ghost_M
+    @test M ≈ reshape(analytic_M, s) ≈ ghost_M
 end
 
 @testset "Test Left Division L4 (fourth order)" begin
@@ -424,7 +436,7 @@ end
     u = sin.(x)
 
     L = CenteredDifference(4, 4, dx, N)
-    Q = RobinBC((1.0, 0.0, sin(0.0)), (1.0, 0.0, sin(0.2+dx)), dx)
+    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, sin(0.2+dx)), dx)
     A = L*Q
 
     analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
@@ -433,25 +445,63 @@ end
     analytic_Qb = [zeros(N+1); sin(0.2+dx)]
     analytic_Ab = analytic_L*analytic_Qb
 
-    analytic_u = analytic_AL \ (u - analytic_Ab)
+
+    analytic_QM = zeros((length(x)+2)*3, length(x)*3)
+    interior = CartesianIndices(Tuple([2:length(x)+1, 1:3]))
+    I1 = CartesianIndex(1,0)
+    for I in interior
+        i = DiffEqOperators.cartesian_to_linear(I, (length(x)+2, 3)) #helper function, see utils.jl
+        j = DiffEqOperators.cartesian_to_linear(I-I1, (length(x), 3))
+        analytic_QM[i,j] = 1.0
+    end
+    analytic_Am = kron(Diagonal(ones(3)), analytic_L)*analytic_QM
+
+    analytic_u = analytic_AL \ Vector(u .- analytic_Ab)
     ghost_u = A \ u
 
-    # Check that A\u.(x) is consistent with analytic_AL \ u.(x)
-    @test analytic_u ≈ ghost_u
+    As_l, As_b = sparse(A, length(u))
 
-    # Check ldiv!
-    u_temp = zeros(N)
-    ldiv!(u_temp, A, u)
-    @test u_temp ≈ ghost_u ≈ analytic_u
+    @test As_l ≈ analytic_AL #This is true
+    @test As_b ≈  analytic_Ab #This is true
+    @test Vector(u .- analytic_Ab) ≈ Vector(u .- As_b) #this is true
+    u_translated = u .- analytic_Ab |> Vector
+    # TODO validate analytic_AL as having reasonable condition (what is reasonable condition?)
+    @test norm(analytic_AL \ u_translated - As_l \ u_translated) / norm(analytic_AL \ u_translated) < cond(analytic_AL)*eps()
+
+    sp_u = As_l\Vector(u .- As_b)
+    # Check that A\u.(x) is consistent with analytic_AL \ u.(x)
+    @test ghost_u ≈ sp_u
+    # TODO make cond function in ghost_derivative_operator?
+    # however, note: ArgumentError: 2-norm condition number is not implemented for sparse matrices, try cond(Array(A), 2) instead
+    # so the actual function would probably just propagate that error
+    _cond(A::GhostDerivativeOperator, u) = cond(sparse(A,length(u))[1] |> Array)
+    @test norm(analytic_u - ghost_u) / norm(ghost_u) < _cond(A,u) * eps() #
 
-    # Check \ for Matrix
     M2 = [u 2.0*u 10.0*u]
-    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
-    ghost_M = A \ M2
-    @test analytic_M ≈ ghost_M
-
-    # Check ldiv! for Matrix
-    M_temp = zeros(N,3)
-    ldiv!(M_temp, A, M2)
-    @test M_temp ≈ ghost_M ≈ analytic_M
+    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
+
+@testset "Test Operator and BC combinations" begin
+    N = 40
+    x = range(-pi, stop = pi, length=N)
+    Δx = x[2]-x[1]
+    u₀=1.0
+    Γ=1.0
+    Dx=u₀*CenteredDifference{1}(1,4,Δx,N)
+    Dxx=Γ*CenteredDifference{1}(2,4,Δx,N)
+    Q=NeumannBC((2x[1]+2, 2x[end]+2), Δx, 6)
+
+    A = Dx*Q + Dxx*Q
+
+    y = A*(x.^2)
+    analytic_y = 2x.+2
+
+    @test_broken y ≈ analytic_y
 end
diff --git a/test/MultiDimBC_test.jl b/test/multi_dim_bc_test.jl
similarity index 69%
rename from test/MultiDimBC_test.jl
rename to test/multi_dim_bc_test.jl
index d2092cce6..0a7dc374f 100644
--- a/test/MultiDimBC_test.jl
+++ b/test/multi_dim_bc_test.jl
@@ -1,8 +1,10 @@
-using LinearAlgebra, DiffEqOperators, Random, Test
+using LinearAlgebra, SparseArrays, DiffEqOperators, Random, Test
 ################################################################################
 # Test 2d extension
 ################################################################################
 
+Random.seed!(7373)
+
 #Create Array
 n = 8
 m = 15
@@ -53,18 +55,23 @@ BCz = fill(Dirichlet0BC(Float64), (n,m))
 
 Qx = MultiDimBC{1}(BCx)
 Qy = MultiDimBC{2}(BCy)
-Qz = Dirichlet0BC{3}(Float64, size(A)) #Test the other constructor
-
-Q1 = (Qx+Qy+Qz)
-Q2 = compose(Qx,Qy,Qz) #test addition combinations
-@test_broken Q1 == Q2 #This fails
-@test all(Q1.BCs .== Q2.BCs) # but this passes so it does actually work
-@test_broken Qtmp = Qx + Qz
-@test_skip Qz+Qx+Qy == Qy+Qtmp
+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
@@ -79,16 +86,19 @@ for i in 1:n, j in 1:m
 end
 
 #test compositions to higher dimension
-for N in 2:7
+for N in 2:6
     sizes = rand(4:7, N)
-    A = rand(sizes...)
+    local A = rand(sizes...)
 
-    Q1_N = Neumann0BC(Float64, Tuple(ones(N)), 3.0, size(A))
+    Q1_N = RobinBC(Tuple(rand(3)), Tuple(rand(3)), fill(0.1, N), 4.0, size(A))
 
-    Q = compose(Q1_N...)
+    local Q = compose(Q1_N...)
 
     A1_N = Q1_N.*fill(A, N)
 
-    A_extended = Q*A
-    @test Array(A_extended) == Array(compose(A1_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
diff --git a/test/runtests.jl b/test/runtests.jl
index d0ffa4dfb..261f07c9c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,7 +20,7 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Validate Regular Derivative Operators" begin include("regular_operator_validation.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 "Validate Higher Dimensional Boundary Extension" begin include("multi_dim_bc_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
@@ -33,6 +33,7 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
     @time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end
     @time @safetestset "Basic SDO Examples" begin include("BasicSDOExamples.jl") end
+    @time @safetestset "3D laplacian Test" begin include("3D_laplacian.jl") end
     # @time @safetestset "Linear Complementarity Problem Examples" begin include("lcp.jl"); include("lcp_split.jl") end
 end