diff --git a/Project.toml b/Project.toml
index e32bf89..10125cd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ContinuumArrays"
 uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
-version = "0.18.7"
+version = "0.19"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -31,14 +31,14 @@ ArrayLayouts = "1.0"
 BandedMatrices = "1"
 BlockArrays = "1"
 DomainSets = "0.6, 0.7"
-FastTransforms = "0.15, 0.16"
+FastTransforms = "0.15, 0.16, 0.17"
 FillArrays = "1.0"
-InfiniteArrays = "0.14, 0.15"
+InfiniteArrays = "0.15"
 Infinities = "0.1"
 IntervalSets = "0.7"
 LazyArrays = "2"
 Makie = "0.20, 0.21, 0.22"
-QuasiArrays = "0.11.8"
+QuasiArrays = "0.12"
 RecipesBase = "1.0"
 StaticArrays = "1.0"
 julia = "1.10"
diff --git a/docs/src/index.md b/docs/src/index.md
index 082a3ad..ffa9e94 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -20,12 +20,35 @@ product of the specified sizes, similar to `plan_fft`.
 5. `plotgrid(::MyBasis, n...)`: return `n` grid points suitable for plotting the basis. The default value for `n` is 10,000. 
 
 
-## Routines
+## Differential Operators
 
 
 ```@docs
 Derivative
 ```
+
+```@docs
+Laplacian
+```
+
+```@docs
+AbsLaplacian
+```
+
+```@docs
+laplacian
+```
+
+```@docs
+abslaplacian
+```
+
+```@docs
+weaklaplacian
+```
+
+## Routines
+
 ```@docs
 transform
 ```
diff --git a/src/ContinuumArrays.jl b/src/ContinuumArrays.jl
index 0c3d19f..d440d13 100644
--- a/src/ContinuumArrays.jl
+++ b/src/ContinuumArrays.jl
@@ -19,11 +19,12 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize, _cutdim,
                     AbstractQuasiFill, UnionDomain, sum_size, sum_layout, _cumsum, cumsum_layout, applylayout, equals_layout, layout_broadcasted, PolynomialLayout, dot_size,
-                    diff_layout, diff_size
+                    diff_layout, diff_size, AbstractQuasiVecOrMat
 import InfiniteArrays: Infinity, InfAxes
 import AbstractFFTs: Plan
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients 
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, .., transform, expand, plan_transform, basis, coefficients,
+        weaklaplacian, laplacian, Laplacian, AbsLaplacian, abslaplacian
 
 
 
@@ -107,7 +108,7 @@ include("plotting.jl")
 
 sum_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = _sum(expand(a), dims)
 dot_size(::InfiniteCardinal{1}, a, b) = dot(expand(a), expand(b))
-diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, dims) = diff(expand(a); dims=dims)
+diff_size(::Tuple{InfiniteCardinal{1}, Vararg{Integer}}, a, order...; dims...) = diff(expand(a), order...; dims...)
 function copy(d::Dot{<:ExpansionLayout,<:ExpansionLayout,<:AbstractQuasiArray,<:AbstractQuasiArray})
     a,b = d.A,d.B
     P,c = basis(a),coefficients(a)
diff --git a/src/bases/bases.jl b/src/bases/bases.jl
index 91709cf..6079bc6 100644
--- a/src/bases/bases.jl
+++ b/src/bases/bases.jl
@@ -421,6 +421,8 @@ basis_axes(ax, v) = error("Overload for $ax")
 coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any}}) where N = v.args[2]
 coefficients(v::ApplyQuasiArray{<:Any,N,typeof(*),<:Tuple{Any,Any,Vararg{Any}}}) where N = ApplyArray(*, tail(v.args)...)
 
+coefficients(B::Basis{T}) where T = SquareEye{T}((axes(B,2),))
+
 
 function unweighted(lay::ExpansionLayout, a)
     wP,c = arguments(lay, a)
@@ -660,29 +662,45 @@ cumsum_layout(::ExpansionLayout, A, dims) = cumsum_layout(ApplyLayout{typeof(*)}
 ###
 # diff
 ###
+diff_layout(::AbstractBasisLayout, Vm, order...; dims...) = error("Overload diff(::$(typeof(Vm)))")
+function diff_layout(::AbstractBasisLayout, a, order::Int; dims...)
+    order < 0 && throw(ArgumentError("order must be non-negative"))
+    order == 0 && return a
+    isone(order) ? diff(a) : diff(diff(a), order-1)
+end
 
-diff_layout(::AbstractBasisLayout, Vm, dims...) = error("Overload diff(::$(typeof(Vm)))")
+function diff_layout(::SubBasisLayout, Vm, order::Int; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    diff(parent(Vm), order)[:,parentindices(Vm)[2]]
+end
 
-function diff_layout(::SubBasisLayout, Vm, dims::Integer=1)
+function diff_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
-    diff(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
+    diff(parent(Vm), order...)[:,parentindices(Vm)[2]]
 end
 
-function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, dims::Integer=1)
+
+function diff_layout(::WeightedBasisLayout{SubBasisLayout}, Vm, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
     w = weight(Vm)
     V = unweighted(Vm)
-    view(diff(w .* parent(V)), parentindices(V)...)
+    view(diff(w .* parent(V), order...), parentindices(V)...)
 end
 
-function diff_layout(::MappedBasisLayouts, V, dims::Integer=1)
-    kr = basismap(V)
-    @assert kr isa AbstractAffineQuasiVector
-    D = diff(demap(V); dims=dims)
+diff_layout(::MappedBasisLayouts, V, order::Int; dims...) = diff_mapped(basismap(V), V, order::Int; dims...)
+diff_layout(::MappedBasisLayouts, V, order...; dims...) = diff_mapped(basismap(V), V, order...; dims...)
+
+function diff_mapped(kr::AbstractAffineQuasiVector, V; dims...)
+    D = diff(demap(V); dims...)
     view(basis(D), kr, :) * (kr.A*coefficients(D))
 end
 
-diff_layout(::ExpansionLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)
+function diff_mapped(kr::AbstractAffineQuasiVector, V, order::Int; dims...)
+    D = diff(demap(V), order; dims...)
+    view(basis(D), kr, :) * (kr.A^order*coefficients(D))
+end
+
+diff_layout(::ExpansionLayout, A, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
 
 
 ####
@@ -708,6 +726,57 @@ function grammatrix_layout(::MappedBasisLayouts, P)
     grammatrix(Q)/kr.A
 end
 
+"""
+    abslaplacian(A, α=1)
+
+computes ``|Δ|^α * A``. 
+"""
+abslaplacian(A, order...; dims...) = abslaplacian_layout(MemoryLayout(A), A, order...; dims...)
+abslaplacian_layout(layout, A, order...; dims...) = abslaplacian_axis(axes(A,1), A, order...; dims...)
+abslaplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = -diff(A, 2order; dims...)
+
+"""
+    abslaplacian(A, k=1)
+
+computes ``Δ^k * A``. 
+"""
+laplacian(A, order...; dims...) = laplacian_layout(MemoryLayout(A), A, order...; dims...)
+laplacian_layout(layout, A, order...; dims...) = laplacian_axis(axes(A,1), A, order...; dims...)
+laplacian_axis(::Inclusion{<:Number}, A, order=1; dims...) = diff(A, 2order; dims...)
+
+
+laplacian_layout(::ExpansionLayout, A, order...; dims...) = laplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
+abslaplacian_layout(::ExpansionLayout, A, order...; dims...) = abslaplacian_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
+
+function abslaplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    abslaplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
+end
+
+function laplacian_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    laplacian(parent(Vm), order...)[:,parentindices(Vm)[2]]
+end
+
+function laplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
+    a = arguments(LAY, V)
+    dims == 1 || throw(ArgumentError("cannot take laplacian a vector along dimension $dims"))
+    *(laplacian(a[1], order...), tail(a)...)
+end
+
+function abslaplacian_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
+    a = arguments(LAY, V)
+    dims == 1 || throw(ArgumentError("cannot take abslaplacian a vector along dimension $dims"))
+    *(abslaplacian(a[1], order...), tail(a)...)
+end
+
+
+
+"""
+    weaklaplacian(A)
+
+represents the weak Laplacian.
+"""
 weaklaplacian(A) = weaklaplacian_layout(MemoryLayout(A), A)
 weaklaplacian_layout(_, A) = weaklaplacian_axis(axes(A,1), A)
 weaklaplacian_axis(::Inclusion{<:Number}, A) = -(diff(A)'diff(A))
diff --git a/src/bases/basisconcat.jl b/src/bases/basisconcat.jl
index 448c02f..5ecf01e 100644
--- a/src/bases/basisconcat.jl
+++ b/src/bases/basisconcat.jl
@@ -8,9 +8,9 @@ abstract type AbstractConcatBasis{T} <: Basis{T} end
 
 copy(S::AbstractConcatBasis) = S
 
-function diff(S::AbstractConcatBasis; dims::Integer)
+function diff(S::AbstractConcatBasis, order...; dims::Integer=1)
     dims == 1 || error("not implemented")
-    args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args; dims=dims))
+    args = arguments.(Ref(ApplyLayout{typeof(*)}()), diff.(S.args, order...; dims=dims))
     all(length.(args) .== 2) || error("Not implemented")
     concatbasis(S, map(first, args)...) * mortar(Diagonal([map(last, args)...]))
 end
@@ -112,4 +112,4 @@ function QuasiArrays._getindex(::Type{IND}, A::HvcatBasis{T}, (x,j)::IND) where
 end
 
 
-diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}; dims::Integer=1) = hcat((diff.(H.args; dims=dims))...)
\ No newline at end of file
+diff(H::ApplyQuasiMatrix{<:Any,typeof(hcat)}, order...; dims::Integer=1) = hcat((diff.(H.args, order...; dims=dims))...)
\ No newline at end of file
diff --git a/src/operators.jl b/src/operators.jl
index b017227..bd29f67 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -108,51 +108,97 @@ show(io::IO, ::MIME"text/plain", δ::DiracDelta) = show(io, δ)
 # Differentiation
 #########
 
+abstract type AbstractDifferentialQuasiMatrix{T} <: LazyQuasiMatrix{T} end
+
 """
 Derivative(axis)
 
 represents the differentiation (or finite-differences) operator on the
 specified axis.
 """
-struct Derivative{T,D} <: LazyQuasiMatrix{T}
-    axis::Inclusion{T,D}
+struct Derivative{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
 end
 
-Derivative{T}(axis::Inclusion{<:Any,D}) where {T,D} = Derivative{T,D}(axis)
-Derivative{T}(domain) where T = Derivative{T}(Inclusion(domain))
+"""
+Laplacian(axis)
+
+represents the laplacian operator `Δ` on the
+specified axis.
+"""
+struct Laplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
+end
 
-Derivative(L::AbstractQuasiMatrix) = Derivative(axes(L,1))
+"""
+AbsLaplacian(axis)
 
-show(io::IO, a::Derivative) = summary(io, a)
-function summary(io::IO, D::Derivative)
-    print(io, "Derivative(")
-    summary(io,D.axis)
-    print(io,")")
+represents the positive-definite/negative/absolute-value laplacian operator `|Δ| ≡ -Δ` on the
+specified axis.
+"""
+struct AbsLaplacian{T,D<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
+    axis::D
+    order::Order
 end
 
-axes(D::Derivative) = (D.axis, D.axis)
-==(a::Derivative, b::Derivative) = a.axis == b.axis
-copy(D::Derivative) = Derivative(copy(D.axis))
+operatororder(D) = something(D.order, 1)
+
+show(io::IO, a::AbstractDifferentialQuasiMatrix) = summary(io, a)
+axes(D::AbstractDifferentialQuasiMatrix) = (D.axis, D.axis)
+==(a::AbstractDifferentialQuasiMatrix, b::AbstractDifferentialQuasiMatrix) = a.axis == b.axis && operatororder(a) == operatororder(b)
+copy(D::AbstractDifferentialQuasiMatrix) = D
+
+
 
-@simplify function *(D::Derivative, B::AbstractQuasiMatrix)
+@simplify function *(D::AbstractDifferentialQuasiMatrix, B::AbstractQuasiVecOrMat)
     T = typeof(zero(eltype(D)) * zero(eltype(B)))
-    diff(convert(AbstractQuasiMatrix{T}, B); dims=1)
+    operatorcall(D, convert(AbstractQuasiArray{T}, B), D.order)
 end
 
-@simplify function *(D::Derivative, B::AbstractQuasiVector)
-    T = typeof(zero(eltype(D)) * zero(eltype(B)))
-    diff(convert(AbstractQuasiVector{T}, B))
+^(D::AbstractDifferentialQuasiMatrix{T}, k::Integer) where T = similaroperator(D, D.axis, operatororder(D) .* k)
+
+function view(D::AbstractDifferentialQuasiMatrix, kr::Inclusion, jr::Inclusion)
+    @boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
+    D
 end
 
+operatorcall(D::AbstractDifferentialQuasiMatrix, B, order) = operatorcall(D)(B, order)
+operatorcall(D::AbstractDifferentialQuasiMatrix, B, ::Nothing) = operatorcall(D)(B)
 
 
+operatorcall(::Derivative) = diff
+operatorcall(::Laplacian) = laplacian
+operatorcall(::AbsLaplacian) = abslaplacian
 
-^(D::Derivative, k::Integer) = ApplyQuasiArray(^, D, k)
 
+for Op in (:Derivative, :Laplacian, :AbsLaplacian)
+    @eval begin
+        $Op{T, D}(axis::D, order) where {T,D<:Inclusion} = $Op{T,D,typeof(order)}(axis, order)
+        $Op{T, D}(axis::D) where {T,D<:Inclusion} = $Op{T,D,Nothing}(axis, nothing)
+        $Op{T}(axis::D, order...) where {T,D<:Inclusion} = $Op{T,D}(axis, order...)
+        $Op{T}(domain, order...) where T = $Op{T}(Inclusion(domain), order...)
 
-function view(D::Derivative, kr::Inclusion, jr::Inclusion)
-    @boundscheck axes(D,1) == kr == jr || throw(BoundsError(D,(kr,jr)))
-    D
+        $Op(ax::AbstractQuasiVector{T}, order...) where T = $Op{eltype(eltype(ax))}(ax, order...)
+        $Op(L::AbstractQuasiMatrix, order...) = $Op(axes(L,1), order...)
+
+        similaroperator(D::$Op, ax, ord) = $Op{eltype(D)}(ax, ord)
+
+        simplifiable(::typeof(*), A::$Op, B::$Op) = Val(true)
+        *(D1::$Op, D2::$Op) = similaroperator(convert(AbstractQuasiMatrix{promote_type(eltype(D1),eltype(D2))}, D1), D1.axis, operatororder(D1)+operatororder(D2))
+
+
+        function summary(io::IO, D::$Op)
+            print(io, "$($Op)(")
+            summary(io, D.axis)
+            if !isnothing(D.order)
+                print(io, ", ")
+                print(io, D.order)
+            end
+            print(io,")")
+        end
+    end
 end
 
 # struct Multiplication{T,F,A} <: AbstractQuasiMatrix{T}
@@ -161,11 +207,17 @@ end
 # end
 
 
-const Identity{T,D} = QuasiDiagonal{T,Inclusion{T,D}}
-
-Identity(d::Inclusion) = QuasiDiagonal(d)
-
 struct OperatorLayout <: AbstractLazyLayout end
-MemoryLayout(::Type{<:Derivative}) = OperatorLayout()
+MemoryLayout(::Type{<:AbstractDifferentialQuasiMatrix}) = OperatorLayout()
 # copy(M::Mul{OperatorLayout, <:ExpansionLayout}) = simplify(M)
 # copy(M::Mul{OperatorLayout, <:AbstractLazyLayout}) = M.A * expand(M.B)
+
+
+# Laplacian
+
+abs(Δ::Laplacian{T}) where T = AbsLaplacian{T}(axes(Δ,1), Δ.order)
+-(Δ::Laplacian{<:Any,<:Any,Nothing}) = abs(Δ)
+-(Δ::AbsLaplacian{T,<:Any,Nothing}) where T = Laplacian{T}(Δ.axis)
+
+^(Δ::AbsLaplacian{T}, k::Real) where T = AbsLaplacian{T}(Δ.axis, operatororder(Δ)*k)
+^(Δ::AbsLaplacian{T}, k::Integer) where T = AbsLaplacian{T}(Δ.axis, operatororder(Δ)*k)
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index a2b8c63..6da07b3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -39,7 +39,7 @@ import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, App
         D = Derivative(x)
         @test D == Derivative{Float64}(x) == Derivative{Float64}(-1..1)
         @test D*x ≡ QuasiOnes(x)
-        @test D^2 * x ≡ QuasiZeros(x)
+        @test D^2 * x == zeros(x)
         @test D*[x D*x] == [D*x D^2*x]
         @test stringmime("text/plain", D) == "Derivative(Inclusion($(-1..1)))"
         @test_throws DimensionMismatch Derivative(Inclusion(0..1)) * x
diff --git a/test/test_splines.jl b/test/test_splines.jl
index fb9771d..68466c4 100644
--- a/test/test_splines.jl
+++ b/test/test_splines.jl
@@ -1,6 +1,6 @@
 using ContinuumArrays, LinearAlgebra, Base64, FillArrays, QuasiArrays, BandedMatrices, BlockArrays, Test
 using QuasiArrays: ApplyQuasiArray, ApplyStyle, MemoryLayout, mul, MulQuasiMatrix, Vec
-import LazyArrays: MulStyle, LdivStyle, arguments, applied, apply
+import LazyArrays: MulStyle, LdivStyle, arguments, applied, apply, simplifiable
 import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, SubBasisLayout, AdjointMappedBasisLayouts, MappedBasisLayout, plan_grid_transform, weaklaplacian
 
 @testset "Splines" begin
@@ -45,6 +45,8 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
             @test sum(H/H \ exp.(x)) ≈ ℯ-1 atol=1E-5
             @test last(cumsum(H/H \ exp.(x))) ≈ sum(H/H\exp.(x))
         end
+
+        @test coefficients(H) ≡ Eye(size(H,2))
     end
 
     @testset "LinearSpline" begin
@@ -163,6 +165,28 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
         @test typeof(diff(L)) == typeof(diff(L; dims=1)) == typeof(D*L)
         @test_throws ErrorException diff(L; dims=2)
 
+        @test diff(L[:,1:2])[1.1,:] == diff(L)[1.1,1:2]
+
+        @test diff(L,0) ≡ L
+        @test diff(f,0) ≡ f
+        @test diff(L,2)[1.1,:] == laplacian(L)[1.1,:] == -abslaplacian(L)[1.1,:] == laplacian(L,1)[1.1,:] == -abslaplacian(L,1)[1.1,:]
+
+        @test diff(L[:,1:2],2)[1.1,:] == diff(L,2)[1.1,1:2]
+        @test diff(f,2)[1.1] == laplacian(f)[1.1] == laplacian(f,1)[1.1] == -abslaplacian(f)[1.1] == -abslaplacian(f,1)[1.1]
+
+        @test laplacian(L[:,1:2])[1.1,:] == laplacian(L)[1.1,1:2] == -abslaplacian(L[:,1:2])[1.1,:] == -abslaplacian(L)[1.1,1:2]
+
+        Δ = Laplacian(L)
+        @test (Δ * L)[1.1,:] == -(abs(Δ) * L)[1.1,:] == laplacian(L)[1.1,:]
+        @test simplifiable(*, Δ, L) == simplifiable(*, abs(Δ), L) == Val(true)
+        @test -abs(Δ) == Δ
+        @test -Δ == abs(Δ)
+        @test Δ^2 == Δ*Δ
+        @test abs(Δ)^2 == abs(Δ^2) == abs(Δ)^2.0
+        @test simplifiable(*, Δ, Δ) == simplifiable(*, abs(Δ), abs(Δ)) == Val(true)
+        @test summary(Δ) == "Laplacian(Inclusion(1 .. 3))"
+        @test summary(Δ^2) == "Laplacian(Inclusion(1 .. 3), 2)"
+
         M = applied(*, (D*L).args..., [1,2,4])
         @test eltype(materialize(M)) == Float64
 
@@ -404,10 +428,12 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout,
         @test @inferred(grid(L[y,:])) ≈ (grid(L) .+ 1) ./ 2
 
         D = Derivative(x)
-        @test (D*L[y,:])[0.1,1] ≈ -9
+        @test (D*L[y,:])[0.1,1] ≈ diff(L[y,:])[0.1,1] ≈ -9
         @test H[y,:] \ (D*L[y,:]) isa BandedMatrix
         @test H[y,:] \ (D*L[y,:]) ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]
 
+        @test all(iszero,diff(L[y,:],2)[0.1,:])
+
         B = L[y,2:end-1]
         @test MemoryLayout(typeof(B)) isa MappedBasisLayout
         @test B[0.1,1] == L[2*0.1-1,2]