diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index 33242f58..3d741474 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -2,7 +2,7 @@ module Interpolations
 
 using Base.Cartesian
 
-import Base: size, eltype, getindex
+import Base: size, eltype, getindex, ndims
 
 export 
     Interpolation,
@@ -50,8 +50,10 @@ end
 # However, all prefilters copy the array, so do that here as well
 prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A)
 
-size(itp::Interpolation, d::Integer) = size(itp.coefs, d)
-size(itp::Interpolation) = size(itp.coefs)
+size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) =
+    size(itp.coefs, d) - 2*padding(IT())
+size(itp::Interpolation) = tuple(collect([size(itp,i) for i in 1:ndims(itp)]))
+ndims(itp::Interpolation) = ndims(itp.coefs)
 eltype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
@@ -68,6 +70,14 @@ include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
 
+# If nothing else is specified, don't pad at all
+padding(::InterpolationType) = 0
+padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT())
+function similar_with_padding(A, it::InterpolationType)
+    pad = padding(it)
+    coefs = Array(eltype(A), [size(A,i)+2pad for i in 1:ndims(A)]...)
+    coefs, pad
+end
 
 # This creates getindex methods for all supported combinations
 for IT in (
@@ -130,22 +140,28 @@ for IT in (
         Quadratic{LinearBC,OnCell}
     )
     @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
-        ret = similar(A)
+        ret, pad = similar_with_padding(A,it)
         szs = collect(size(A))
         strds = collect(strides(A))
+        strdsR = collect(strides(ret))
 
         for dim in 1:N
             n = szs[dim]
             szs[dim] = 1
 
-            M = prefiltering_system_matrix(eltype(A), n, it)
+            M, b = prefiltering_system(eltype(A), n+2pad, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i
                 strt = 1 + sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
+                strtR = 1 + sum([(cc[i]-1)*strdsR[i] for i in 1:length(cc)])
                 rng = range(strt, strds[dim], n)
-                ret[rng] = M \ vec(A[rng])
-            end            
+                rngR = range(strtR, strdsR[dim], size(ret,dim))
+
+                b[1+pad:end-pad] += A[rng]
+                ret[rngR] = M \ b
+                b[1+pad:end-pad], A[rng]
+            end
             szs[dim] = n
         end
         ret
diff --git a/src/quadratic.jl b/src/quadratic.jl
index fd2b80c8..48a7c845 100644
--- a/src/quadratic.jl
+++ b/src/quadratic.jl
@@ -12,106 +12,12 @@ end
 
 function indices(q::Quadratic, N)
     quote
+        pad = padding($q)
         @nexprs $N d->begin
             ixp_d = ix_d + 1
             ixm_d = ix_d - 1
 
-            # fx_d is the distance from the middle point to the interpolation point
-            fx_d = x_d - convert(typeof(x_d), ix_d)
-        end
-    end
-end
-
-function bc_gen(::Quadratic{Flat,OnGrid}, N)
-    quote
-        # After extrapolation has been handled, 1 <= x_d <= size(itp,d)
-        # The index is simply the closest integer.
-         @nexprs $N d->(ix_d = iround(x_d))
-        # end
-    end
-end
-
-function bc_gen(::Quadratic{Flat,OnCell}, N)
-    quote
-        # After extrapolation has been handled, 0.5 <= x_d <= size(itp,d)+.5
-        @nexprs $N d->begin
-            # The index is the closest integer...
-            ix_d = iround(x_d)
-
-            #...except in the case where x_d is actually at the upper edge;
-            # since size(itp,d)+.5 is rounded toward size(itp,d)+1,
-            # it needs special treatment*.
-            if x_d == size(itp,d)+.5
-                ix_d -= 1
-            end
-        end
-    end
-end
-function bc_gen(::Quadratic{LinearBC,OnCell}, N)
-    quote
-         @nexprs $N d->begin
-            if x_d < 1
-                # extrapolate towards -∞
-                fx_d = x_d - convert(typeof(x_d), 1)
-                k = itp[2] - itp[1]
-
-                return itp[1] + k * fx_d
-            end
-            if x_d > size(itp, d)
-                # extrapolate towards ∞
-                s_d = size(itp, d)
-                fx_d = x_d - convert(typeof(x_d), s_d)
-                k = itp[s_d] - itp[s_d - 1]
-
-                return itp[s_d] + k * fx_d
-            end
-
-            ix_d = iround(x_d)
-         end
-    end
-end
-
-function indices(::Quadratic{Flat,OnGrid}, N)
-    quote
-        @nexprs $N d->begin
-            ixp_d = ix_d + 1
-            ixm_d = ix_d - 1
-
-            fx_d = x_d - convert(typeof(x_d), ix_d)
-
-            if ix_d == size(itp,d)
-                ixp_d = ixm_d
-            end
-            if ix_d == 1
-                ixm_d = ixp_d
-            end
-        end
-    end
-end
-function indices(::Quadratic{Flat,OnCell}, N)
-    quote
-        @nexprs $N d->begin
-            ixp_d = ix_d + 1
-            ixm_d = ix_d - 1
-
-            fx_d = x_d - convert(typeof(x_d), ix_d)
-
-            if ix_d == size(itp,d)
-                ixp_d = ix_d
-            end
-            if ix_d == 1
-                ixm_d = ix_d
-            end
-        end
-    end
-end
-function indices(::Quadratic{LinearBC,OnCell}, N)
-    quote
-        @nexprs $N d->begin
-            ixp_d = ix_d + 1
-            ixm_d = ix_d - 1
-
-            fx_d = x_d - convert(typeof(x_d), ix_d)
+            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
         end
     end
 end
@@ -126,7 +32,7 @@ function coefficients(::Quadratic, N)
     end
 end
 
-# This assumes integral values ixm_d, ix_d, and ixp_d (typically ixm_d = ix_d-1, ixp_d = ix_d+1, except at boundaries),
+# This assumes integral values ixm_d, ix_d, and ixp_d,
 # coefficients cm_d, c_d, and cp_d, and an array itp.coefs
 function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
     if length(offsets) < N
@@ -140,47 +46,58 @@ function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
     end
 end
 
-function unmodified_system_matrix{T}(::Type{T}, n::Int, ::Quadratic)
+# Quadratic interpolation has 1 extra coefficient at each boundary
+# Therefore, make the coefficient array 1 larger in each direction,
+# in each dimension.
+padding(::Quadratic) = 1
+
+function inner_system_diags{T}(::Type{T}, n::Int, ::Quadratic)
     du = fill(convert(T,1/8), n-1)
     d = fill(convert(T,3/4),n)
     dl = copy(du)
+    d[1] = d[end] = du[1] = dl[end] = convert(T,0)
     (dl,d,du)
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{ExtendInner})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    d[1] = d[end] = 9/8
-    du[1] = dl[end] = -1/4
-    MT = lufact!(Tridiagonal(dl, d, du))
-    U = zeros(T,n,2)
-    V = zeros(T,2,n)
-    C = zeros(T,2,2)
-
-    C[1,1] = C[2,2] = 1/8
-    U[1,1] = U[n,2] = 1.
-    V[1,3] = V[2,n-2] = 1.
-
-    Woodbury(MT, U, C, V)
-end
-
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    d[1] += 1/8
-    d[end] += 1/8
-    lufact!(Tridiagonal(dl, d, du))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = -1
+    du[1] = dl[end] = 1
+    lufact!(Tridiagonal(dl, d, du)), zeros(T, n)
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-    du[1] += 1/8
-    dl[end] += 1/8
-    lufact!(Tridiagonal(dl, d, du))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = convert(T, -1)
+    du[1] = dl[end] = convert(T, 0)
+
+    rowspec = zeros(T,n,2)
+    # first row     last row
+    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    colspec = zeros(T,2,n)
+    # third col     third-to-last col
+    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    valspec = zeros(T,2,2)
+    # val for [1,3], val for [n,n-2]
+    valspec[1,1] = valspec[2,2] = convert(T, 1)
+
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{LinearBC,OnCell})
-    dl,d,du = unmodified_system_matrix(T,n,q)
-
-    d[1] = d[end] = 1
-    du[1] = dl[end] = 0
-    lufact!(Tridiagonal(dl, d, du))
+function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{LinearBC})
+    dl,d,du = inner_system_diags(T,n,q)
+    d[1] = d[end] = convert(T, 1)
+    du[1] = dl[end] = convert(T, -2)
+
+    rowspec = zeros(T,n,2)
+    # first row     last row
+    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    colspec = zeros(T,2,n)
+    # third col     third-to-last col
+    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    valspec = zeros(T,2,2)
+    # val for [1,3], val for [n,n-2]
+    valspec[1,1] = valspec[2,2] = convert(T, 1)
+
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
 end