diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl
index fe9918484..117de2b0a 100644
--- a/src/DiffEqOperators.jl
+++ b/src/DiffEqOperators.jl
@@ -30,12 +30,12 @@ include("derivative_operators/multi_dim_bc_operators.jl")
 
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
+include("derivative_operators/coefficients.jl")
 include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/convolutions.jl")
 include("derivative_operators/ghost_derivative_operator.jl")
 include("derivative_operators/derivative_operator_functions.jl")
-include("derivative_operators/coefficient_functions.jl")
 
 ### Composite Operators
 include("composite_operators.jl")
diff --git a/src/derivative_operators/coefficient_functions.jl b/src/derivative_operators/coefficient_functions.jl
deleted file mode 100644
index a3e81cce3..000000000
--- a/src/derivative_operators/coefficient_functions.jl
+++ /dev/null
@@ -1,67 +0,0 @@
-"""
-```
-compute_coeffs(coeff_func, current_coeffs)
-```
-Calculates the coefficients for the stencil of UpwindDifference operators.
-"""
-function compute_coeffs!(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number}
-    return current_coeffs .+= coeff_func
-end
-
-function compute_coeffs!(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number}
-     return current_coeffs[:] += coeff_func
-end
-
-# Coefficient functions when coeff_func is a Function and current_coeffs exists
-function compute_coeffs!(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number}
-    if hasmethod(coeff_func, (Vector{T},))
-        current_coeffs[:] = coeff_func(current_coeffs)
-    else
-        map!(coeff_func, current_coeffs, current_coeffs)
-    end
-    return current_coeffs
-    # if hasmethod(coeff_func, (Int,)) # assume we want to provide the index of coefficients to coeff_func
-    #     for i = 1:length(current_coeffs)
-    #         current_coeffs[i] += coeff_func(i)
-    #     end
-    #     return current_coeffs
-    # elseif hasmethod(coeff_func, (Vector{Int},))  # assume we want to provide the index of coefficients to coeff_func
-    #     current_coeffs[:] += coeff_func(collect(1:length(current_coeffs)))
-    #     return current_coeffs
-    # elseif hasmethod(coeff_func, (UnitRange{Int},)) # assume we want to provide the index of coefficients to coeff_func
-    #     current_coeffs[:] += coeff_func(1:length(current_coeffs))
-    #     return current_coeffs
-    # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients
-    #     map!(coeff_func, current_coeffs, current_coeffs)
-    #     return current_coeffs
-    # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients
-    #     current_coeffs[:] = coeff_func(current_coeffs)
-    #     return current_coeffs
-    # else
-    #     error("Coefficient functions with the arguments of $coeff_func have not been implemented.")
-    # end
-end
-
-compute_coeffs(coeff_func::Number, current_coeffs::AbstractVector{T}) where {T<:Number} = current_coeffs .+ coeff_func
-compute_coeffs(coeff_func::AbstractVector{T}, current_coeffs::AbstractVector{T}) where {T<:Number} = coeff_func + current_coeffs
-
-function compute_coeffs(coeff_func::Function, current_coeffs::AbstractVector{T}) where {T<:Number}
-    if hasmethod(coeff_func, (Vector{T},))
-        return coeff_func(current_coeffs)
-    else
-        return map(coeff_func, current_coeffs)
-    end
-    # if hasmethod(coeff_func, (Int,))
-    #     return current_coeffs + map(coeff_func, collect(1:length(current_coeffs)))
-    # elseif hasmethod(coeff_func, (Vector{Int},))
-    #     return current_coeffs + coeff_func(collect(1:length(current_coeffs)))
-    # elseif hasmethod(coeff_func, (UnitRange{Int},))
-    #     return current_coeffs + coeff_func(1:length(current_coeffs))
-    # elseif hasmethod(coeff_func, (Float64,)) # assume we want coeff_func to operate on existing coefficients
-    #     return map(coeff_func, current_coeffs)
-    # elseif hasmethod(coeff_func, (Vector{Float64},)) # assume we want to coeff_func to operate on existing coefficients
-    #     return coeff_func(current_coeffs)
-    # else
-    #     error("Coefficient functions with the arguments of coeff_func $coeff_func have not been implemented.")
-    # end
-end
diff --git a/src/derivative_operators/coefficients.jl b/src/derivative_operators/coefficients.jl
new file mode 100644
index 000000000..4b4b9cfa4
--- /dev/null
+++ b/src/derivative_operators/coefficients.jl
@@ -0,0 +1,32 @@
+"""
+```
+init_coefficients(coeff_func, len::Int)
+```
+
+Return the initial value of an operator's `coefficients` field based on the type
+of `coeff_func`.
+"""
+init_coefficients(coeff_func::Nothing, len::Int) = nothing
+
+init_coefficients(coeff_func::Number, len::Int) = coeff_func * ones(typeof(coeff_func), len)
+
+function init_coefficients(coeff_func::AbstractVector{T}, len::Int) where T <: Number
+    coeff_func
+end
+
+init_coefficients(coeff_func::Function, len::Int) = ones(Float64, len)
+
+
+
+"""
+```
+get_coefficient(coefficients, index)
+```
+"""
+get_coefficient(coefficients::AbstractVector, index::Int) = coefficients[index]
+
+# FIXME: I don't think this case is used anymore
+get_coefficient(coefficients::Number, index::Int) = coefficients
+
+# FIXME: Why use "true" here for the value 1?
+get_coefficient(coefficients::Nothing, index::Int) = true
diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl
index 0567eb63b..3b0a0380e 100644
--- a/src/derivative_operators/concretization.jl
+++ b/src/derivative_operators/concretization.jl
@@ -12,12 +12,12 @@ function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) wh
     bstl = A.boundary_stencil_length
 
     coeff   = A.coefficients
-    get_coeff = if coeff isa AbstractVector
-        i -> coeff[i]
+    get_coeff(i) = if coeff isa AbstractVector
+        get_coefficient(coeff, i)
     elseif coeff isa Number
-        i -> coeff
+        coeff
     else
-        i -> true
+        true
     end
 
     for i in 1:bl
@@ -451,7 +451,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
     stencils = A.stencil_coefs
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             cur_stencil = stencils
             L[i,i+1:i+stl] = cur_coeff*cur_stencil
@@ -462,7 +462,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         cur_stencil = stencils
         cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
         if cur_coeff >= 0
@@ -473,7 +473,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             cur_stencil = stencils
             cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
@@ -495,7 +495,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
     coeff   = A.coefficients
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
         else
@@ -504,7 +504,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
         else
@@ -513,7 +513,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
         else
@@ -540,7 +540,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
     stencils = A.stencil_coefs
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             cur_stencil = stencils
             L[i,i+1:i+stl] = cur_coeff*cur_stencil
@@ -551,7 +551,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         cur_stencil = stencils
         cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
         if cur_coeff >= 0
@@ -562,7 +562,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             cur_stencil = stencils
             cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
@@ -584,7 +584,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
     coeff   = A.coefficients
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
         else
@@ -593,7 +593,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
         else
@@ -602,7 +602,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
         else
@@ -629,7 +629,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
     stencils = A.stencil_coefs
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             cur_stencil = stencils
             L[i,i+1:i+stl] = cur_coeff*cur_stencil
@@ -640,7 +640,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         cur_stencil = stencils
         cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
         if cur_coeff >= 0
@@ -651,7 +651,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             cur_stencil = stencils
             cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
@@ -674,7 +674,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
     L = BandedMatrix{T}(Zeros(len, len+2), (stl-2, stl))
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
         else
@@ -683,7 +683,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
     end
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
         else
@@ -692,7 +692,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
     end
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
         else
diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl
index 2e2475a51..bd51804e0 100644
--- a/src/derivative_operators/convolutions.jl
+++ b/src/derivative_operators/convolutions.jl
@@ -33,7 +33,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A
         for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
             xtempi = zero(T)
             cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
-            cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+            cur_coeff = get_coefficient(coeff, i)
             for idx in 1:A.stencil_length
                 xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
             end
@@ -43,7 +43,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A
         for i in [(1+A.boundary_point_count):(A.boundary_point_count+offset); (length(x_temp)-A.boundary_point_count-offset+1):(length(x_temp)-A.boundary_point_count)]
             xtempi = zero(T)
             cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
-            cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+            cur_coeff = get_coefficient(coeff, i)
             for idx in 1:A.stencil_length
                 xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
             end
@@ -58,7 +58,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A:
     coeff   = A.coefficients
     for i in 1 : A.boundary_point_count
         cur_stencil = stencil[i]
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i)
         xtempi = cur_coeff*stencil[i][1]*x[1]
         for idx in 2:A.boundary_stencil_length
             xtempi += cur_coeff * cur_stencil[idx] * x[idx]
@@ -73,7 +73,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, x::AbstractVector{T2}, A
     coeff   = A.coefficients
     for i in 1 : A.boundary_point_count
         cur_stencil = stencil[i]
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i)
         xtempi = cur_coeff*stencil[i][end]*x[end]
         for idx in (A.boundary_stencil_length-1):-1:1
             xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
@@ -93,7 +93,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
     for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
         xtempi = zero(T)
         cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i)
         cur_stencil = cur_coeff >= 0 ? cur_stencil : A.derivative_order % 2 == 0 ? reverse(cur_stencil) : -1*reverse(cur_stencil)
         for idx in 1:A.stencil_length
             x_idx = cur_coeff < 0 ? x[i - A.stencil_length + 1 + idx] : x[i + idx]
@@ -163,7 +163,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
     coeff   = A.coefficients
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             xtempi = zero(T)
             cur_stencil = A.stencil_coefs[1,i-bpc]
@@ -190,7 +190,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::D
     coeff   = A.coefficients
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             xtempi = zero(T)
             cur_stencil = A.low_boundary_coefs[1,i]
@@ -218,7 +218,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
     coeff   = A.coefficients
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             xtempi = zero(T)
             cur_stencil = A.high_boundary_coefs[2,i-len+bpc]
@@ -250,7 +250,7 @@ function convolve_interior!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector
     for i in (2+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)-1
         xtempi = zero(T)
         cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i - A.boundary_point_count)
         @inbounds for idx in 1:A.stencil_length
             xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1]
         end
@@ -264,7 +264,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector,
     coeff   = A.coefficients
     for i in 1 : A.boundary_point_count
         cur_stencil = stencil[i]
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i)
         xtempi = cur_coeff*cur_stencil[1]*_x.l
         @inbounds for idx in 2:A.boundary_stencil_length
             xtempi += cur_coeff * cur_stencil[idx] * _x.u[idx-1]
@@ -277,7 +277,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector,
     i = 1 + A.boundary_point_count
     xtempi = zero(T)
     cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs
-    cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
+    cur_coeff = get_coefficient(coeff, i - A.boundary_point_count)
     xtempi = cur_coeff*cur_stencil[1]*_x.l
     @inbounds for idx in 2:A.stencil_length
         xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1]
@@ -296,7 +296,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector
     i = length(x_temp)-A.boundary_point_count
     xtempi = zero(T)
     cur_stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-A.boundary_point_count] : A.stencil_coefs
-    cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i-A.boundary_point_count] : coeff isa Number ? coeff : true
+    cur_coeff = get_coefficient(coeff, i - A.boundary_point_count)
     xtempi = cur_coeff*cur_stencil[end]*_x.r
     @inbounds for idx in 1:A.stencil_length-1
         xtempi += cur_coeff * cur_stencil[idx] * x[(i-1) - (mid-idx) + 1]
@@ -304,7 +304,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T1}, _x::BoundaryPaddedVector
     x_temp[i] = xtempi + !overwrite*x_temp[i]
     for i in 1 : A.boundary_point_count
         cur_stencil = stencil[i]
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[bc_start + i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, bc_start + i)
         xtempi = cur_coeff*cur_stencil[end]*_x.r
         @inbounds for idx in (A.boundary_stencil_length-1):-1:1
             xtempi += cur_coeff * cur_stencil[end-idx] * _x.u[end-idx+1]
@@ -331,7 +331,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
     for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
         xtempi = zero(T)
         cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_point_count] : stencil
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff = get_coefficient(coeff, i)
         cur_stencil = cur_coeff >= 0 ? cur_stencil : A.derivative_order % 2 == 0 ? reverse(cur_stencil) : -1*reverse(cur_stencil)
         for idx in 1:A.stencil_length
             index = cur_coeff < 0 ? i - A.stencil_length + 1 + idx : i + idx
@@ -355,7 +355,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
     x = _x.u
 
     for i in 1:A.boundary_point_count
-        cur_coeff = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         xtempi = 0.0
         if cur_coeff >= 0 && i+A.stencil_length <= length(_x)
             cur_stencil = eltype(upwind_stencils) <: AbstractVector ? upwind_stencils[i] : upwind_stencils
@@ -383,7 +383,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
     _x_len = length(_x)
 
     for i in 1:A.boundary_point_count
-        cur_coeff = coeff[x_temp_len-A.boundary_point_count+i]
+        cur_coeff = get_coefficient(coeff, x_temp_len-A.boundary_point_count+i)
         xtempi = 0.0
         if cur_coeff < 0 && _x_len-A.stencil_length - A.boundary_point_count + i >= 1
             cur_stencil = eltype(downwind_stencils) <: AbstractVector ? downwind_stencils[i] : downwind_stencils
@@ -419,7 +419,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
 
 
     for i in bpc+1:len-bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             xtempi = zero(T)
             cur_stencil = A.stencil_coefs[1,i-bpc]
@@ -450,7 +450,7 @@ function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
     _x_len = length(_x)
 
     for i in 1:bpc
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff >= 0
             xtempi = zero(T)
             cur_stencil = A.low_boundary_coefs[1,i]
@@ -482,7 +482,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector,
     _x_len = length(_x)
 
     for i in len-bpc+1:len
-        cur_coeff   = coeff[i]
+        cur_coeff = get_coefficient(coeff, i)
         if cur_coeff < 0
             xtempi = zero(T)
             cur_stencil = A.high_boundary_coefs[2,i-len+bpc]
diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl
index e87535959..528a5a9fc 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 : fill!(Vector{T}(undef,len),0)
+    coefficients            = init_coefficients(coeff_func, len)
 
     DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
         typeof(low_boundary_coefs),typeof(coefficients),
@@ -103,7 +103,7 @@ function CenteredDifference{N}(derivative_order::Int,
                                                                   calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1]
     high_boundary_coefs      = convert(SVector{boundary_point_count},_high_boundary_coefs)
 
-    coefficients            = coeff_func isa Nothing ? nothing : zeros(T,len)
+    coefficients            = init_coefficients(coeff_func, len)
 
     DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs),
         typeof(low_boundary_coefs),typeof(coefficients),
@@ -128,12 +128,12 @@ constructs a DerivativeOperator that automatically implements upwinding.
 
 ### Inputs
 * `dx::T` or `dx::Vector{T}`: grid spacing
-* `coeff_func`: function mapping index in the grid to coefficient at that grid location
+* `coeff_func`: Same as for DerivativeOperator.
 
 ### Examples
 julia> drift = [1., 1., -1.]
-julia> L1 = UpwindDifference(1, 1, 1., 3, i -> drift[i])
-julia> L2 = UpwindDifference(1, 1, 1., 3, i -> 1.)
+julia> L1 = UpwindDifference(1, 1, 1., 3, drift)
+julia> L2 = UpwindDifference(1, 1, 1., 3, 1.)
 julia> Q = Neumann0BC(1, 1.)
 julia> Array(L1 * Q)[1]
 3×3 Array{Float64,2}:
@@ -161,21 +161,18 @@ function UpwindDifference{N}(derivative_order::Int,
 
     low_boundary_x         = 0.0:(boundary_stencil_length-1)
     L_boundary_deriv_spots = 1.0:boundary_stencil_length - 2.0
-    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots]
-    low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
+    _low_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots]
+    low_boundary_coefs     = convert(SVector{boundary_point_count},_low_boundary_coefs)
 
-    high_boundary_x         = 0.0:-1.0:-(boundary_stencil_length-1.0)
+    high_boundary_x        = 0.0:-1.0:-(boundary_stencil_length-1.0)
     R_boundary_deriv_spots = -1.0:-1.0:-(boundary_stencil_length-2.0)
-    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots]
-    high_boundary_coefs = convert(SVector{boundary_point_count},_high_boundary_coefs)
+    _high_boundary_coefs   = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots]
+    high_boundary_coefs    = convert(SVector{boundary_point_count},_high_boundary_coefs)
 
-    coefficients = zeros(T,len)
-    if coeff_func != nothing
-        compute_coeffs!(coeff_func, coefficients)
-    end
+    coefficients           = init_coefficients(coeff_func, len)
 
     DerivativeOperator{T,N,true,T,typeof(stencil_coefs),
-        typeof(low_boundary_coefs),Vector{T},
+        typeof(low_boundary_coefs),typeof(coefficients),
         typeof(coeff_func)}(
         derivative_order, approximation_order, dx, len, stencil_length,
         stencil_coefs,
@@ -214,14 +211,10 @@ function UpwindDifference{N}(derivative_order::Int,
     _downwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2:i+1])) for i in len-boundary_point_count+1:len])
     high_boundary_coefs = [_upwind_coefs ; _downwind_coefs]
 
-    # Compute coefficients
-    coefficients = zeros(T,len)
-    if coeff_func != nothing
-        compute_coeffs!(coeff_func, coefficients)
-    end
+    coefficients = init_coefficients(coeff_func, len)
 
     DerivativeOperator{T,N,true,typeof(dx),typeof(stencil_coefs),
-        typeof(low_boundary_coefs),Vector{T},
+        typeof(low_boundary_coefs),typeof(coefficients),
         typeof(coeff_func)}(
         derivative_order, approximation_order, dx, len, stencil_length,
         stencil_coefs,
diff --git a/test/coefficient_functions.jl b/test/coefficient_functions.jl
deleted file mode 100644
index 5e7a5cf23..000000000
--- a/test/coefficient_functions.jl
+++ /dev/null
@@ -1,88 +0,0 @@
-using Test, DiffEqOperators
-
-# Set up coefficient functions to test with
-vec_fcn = Vector{Function}(undef,0)
-f1(x::Float64) = sin(x)
-f2(x::Vector{Float64}) = sin.(x)
-push!(vec_fcn, f1)
-push!(vec_fcn, f2)
-
-@testset "Test coefficient functions when current_coeffs exists" begin
-
-    vec_fcn_ans = Vector{Vector{Float64}}(undef,0)
-    for i = 1:2
-        push!(vec_fcn_ans, [sin(1.1), sin(2.3), sin(4.5)])
-    end
-
-    current_coeffs = [1.1, 2.3, 4.5]
-
-    # Test inplace versions
-    current_coeffs1 = Vector{Float64}(undef,3)
-    @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ 2.54 .* ones(3)
-    current_coeffs1 = Vector{Float64}(undef,3)
-    @test DiffEqOperators.compute_coeffs!([1.1, 2.3, 4.5], current_coeffs1) ≈ current_coeffs
-    @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ [3.64, 4.84, 7.04]
-    @test all(current_coeffs .!= current_coeffs1)
-    current_coeffs1 = copy(current_coeffs)
-    @test DiffEqOperators.compute_coeffs!([2.54, 1.03, .18], current_coeffs1) ≈ [3.64, 3.33, 4.68]
-    @test all(current_coeffs .!= current_coeffs1)
-    for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans)
-        current_coeffs1 = copy(current_coeffs)
-        @test DiffEqOperators.compute_coeffs!(fcn, current_coeffs1) ≈ fcn_ans
-        @test all(current_coeffs1 .!= current_coeffs)
-    end
-
-    # Test not in place versions
-    current_coeffs1 = Vector{Float64}(undef,3)
-    @test DiffEqOperators.compute_coeffs(2.54, current_coeffs) ≈ [3.64, 4.84, 7.04]
-    @test DiffEqOperators.compute_coeffs([2.54, 1.03, 0.18], current_coeffs) ≈ [3.64, 3.33, 4.68]
-    for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans)
-        @test DiffEqOperators.compute_coeffs(fcn, current_coeffs) ≈ fcn_ans
-    end
-end
-
-@testset "Check coefficients of DerivativeOperators are properly computed" begin
-
-    # Set up operators to test construction with different coeff_func
-    vec_op = Vector{DerivativeOperator}(undef,0)
-    dxv = [1.4, 1.1, 2.3, 3.6, 1.5]
-    push!(vec_op, UpwindDifference(1, 1, 1., 3, 1.))
-    push!(vec_op, UpwindDifference(1, 1, 1., 3, [1., 2.25, 3.5]))
-    push!(vec_op, UpwindDifference(1, 1, 1., 3, cos))
-    push!(vec_op, UpwindDifference(1, 1, dxv, 3, 1.))
-    push!(vec_op, UpwindDifference(1, 1, dxv, 3, [1., 2.25, 3.5]))
-    push!(vec_op, UpwindDifference(1, 1, dxv, 3, cos))
-
-
-    vec_op_ans = Vector{Any}(undef,0)
-    push!(vec_op_ans, ones(3))
-    push!(vec_op_ans, [1., 2.25, 3.5])
-    push!(vec_op_ans, ones(3))
-    push!(vec_op_ans, ones(3))
-    push!(vec_op_ans, [1., 2.25, 3.5])
-    push!(vec_op_ans, ones(3))
-
-    # Check constructors
-    @test UpwindDifference(1, 1, 1., 3).coefficients + ones(3) - ones(3) ≈ zeros(3)
-    @test CenteredDifference(2, 2, 1., 3, 0.).coefficients + ones(3) - ones(3) ≈ zeros(3)
-    for (L, op_ans) in zip(vec_op, vec_op_ans)
-        @test L.coefficients ≈ op_ans
-    end
-    push!(vec_op, CenteredDifference(2, 2, 1., 3)) == nothing
-
-    # Compute answers to coeff_func * operator
-    func_mul_op_ans1 = Vector{Any}(undef,0)
-    for i = 1:2
-        push!(func_mul_op_ans1, zeros(3))
-    end
-    func_mul_op_ans2 = Vector{Any}(undef,0)
-    for i = 1:2
-        push!(func_mul_op_ans2, sin(1) .* ones(3))
-    end
-    func_mul_op_ans3 = Vector{Any}(undef,0)
-    for i = 1:2
-        push!(func_mul_op_ans3, [0., sin(1.5), 0.])
-    end
-end
-
-nothing
diff --git a/test/runtests.jl b/test/runtests.jl
index 77afc7b78..ca48bc139 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -29,7 +29,6 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Convolutions" begin include("convolutions.jl") end
     @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end
     @time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end
-    @time @safetestset "Coefficient Functions" begin include("coefficient_functions.jl") end
     @time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
     @time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end
     @time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("MOL_1D_Linear_Convection.jl") end