diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl
index df3af2f..ccc7560 100644
--- a/src/MultivariateOrthogonalPolynomials.jl
+++ b/src/MultivariateOrthogonalPolynomials.jl
@@ -32,7 +32,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
        Zernike, ZernikeWeight, zerniker, zernikez,
        PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
-       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
+       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial, kernel,
        grammatrix, oneto
 
 include("ModalInterlace.jl")
@@ -41,5 +41,7 @@ include("disk.jl")
 include("rectdisk.jl")
 include("triangle.jl")
 
+include("kernels.jl")
+
 
 end # module
diff --git a/src/kernels.jl b/src/kernels.jl
new file mode 100644
index 0000000..414e63c
--- /dev/null
+++ b/src/kernels.jl
@@ -0,0 +1,19 @@
+const Kernel{T,F,D1,V,D2} = BroadcastQuasiMatrix{T,F,Tuple{D1,QuasiAdjoint{V,Inclusion{V,D2}}}}
+
+"""
+    kernel(f)
+
+returns `K::AbstractQuasiMatrix` such that `K[x,y] == f[SVector(x,y)]`.
+"""
+function kernel(f::AbstractQuasiVector{<:Real})
+    P², c = f.args
+    P,Q = P².args
+    P * c.array * Q'
+end
+
+@simplify function *(K::Kernel, Q::OrthogonalPolynomial)
+    x,y = axes(K)
+    Px,Py = legendre(x),legendre(y)
+    kernel(expand(RectPolynomial(Px, Py), splat(K.f))) * Q
+end
+
diff --git a/src/rect.jl b/src/rect.jl
index 475bf1a..06f5672 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -49,6 +49,7 @@ end
     RectPolynomial(A,U) * KronTrav(M, Eye{eltype(M)}(∞))
 end
 
+
 function weaklaplacian(P::RectPolynomial)
     A,B = P.args
     Δ_A,Δ_B = weaklaplacian(A), weaklaplacian(B)
@@ -59,7 +60,7 @@ end
 function \(P::RectPolynomial, Q::RectPolynomial)
     PA,PB = P.args
     QA,QB = Q.args
-    KronTrav(PA\QA, PB\QB)
+    KronTrav(PB\QB, PA\QA)
 end
 
 @simplify function *(Ac::QuasiAdjoint{<:Any,<:RectPolynomial}, B::RectPolynomial)
@@ -68,7 +69,11 @@ end
     KronTrav(PA'QA, PB'QB)
 end
 
+"""
+   ApplyPlan(f, plan)
 
+applies a plan and then the function f. 
+"""
 struct ApplyPlan{T, F, Pl}
     f::F
     plan::Pl
@@ -78,7 +83,6 @@ ApplyPlan(f, P) = ApplyPlan{eltype(P), typeof(f), typeof(P)}(f, P)
 
 *(A::ApplyPlan, B::AbstractArray) = A.f(A.plan*B)
 
-basis_axes(d::Inclusion{<:Any,<:ProductDomain}, v) = KronPolynomial(map(d -> basis(Inclusion(d)),components(d.domain))...)
 
 struct TensorPlan{T, Plans}
     plans::Plans
@@ -98,6 +102,8 @@ function checkpoints(P::RectPolynomial)
     SVector.(x, y')
 end
 
+basis_axes(d::Inclusion{<:Any,<:ProductDomain}, v) = KronPolynomial(map(d -> basis(Inclusion(d)),components(d.domain))...)
+
 function plan_grid_transform(P::KronPolynomial{d,<:Any,<:Fill}, B::Tuple{Block{1}}, dims=1:1) where d
     @assert dims == 1
 
@@ -149,4 +155,13 @@ function transform_ldiv(K::KronPolynomial{d,V,<:Fill{<:Legendre}}, f::Union{Abst
     T = KronPolynomial{d}(Fill(ChebyshevT{V}(), size(K.args)...))
     dat = (T \ f).array
     DiagTrav(pad(FastTransforms.th_cheb2leg(paddeddata(dat)), axes(dat)...))
-end
\ No newline at end of file
+end
+
+
+function Base.summary(io::IO, P::RectPolynomial)
+    A,B = P.args
+    summary(io, A)
+    print(io, " ⊗ ")
+    summary(io, B)
+end
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 873aa52..4e0baa8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,4 +6,5 @@ include("test_modalinterlace.jl")
 include("test_disk.jl")
 include("test_rectdisk.jl")
 include("test_triangle.jl")
+include("test_kernels.jl")
 # include("test_dirichlettriangle.jl")
diff --git a/test/test_kernels.jl b/test/test_kernels.jl
new file mode 100644
index 0000000..926f822
--- /dev/null
+++ b/test/test_kernels.jl
@@ -0,0 +1,37 @@
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, FillArrays, LazyBandedMatrices, BlockArrays, StaticArrays, Test
+using BlockBandedMatrices: _BandedBlockBandedMatrix, _BlockBandedMatrix
+using Base: oneto
+
+@testset "Kernels" begin
+    @testset "Basic" begin
+        P = Legendre()
+        T = Chebyshev()
+        P² = RectPolynomial(Fill(P, 2))
+        T² = RectPolynomial(Fill(T, 2))
+
+        k = (x,y) -> exp(x*cos(y))
+        K = kernel(expand(P², splat(k)))
+        @test K[0.1,0.2] ≈ k(0.1,0.2)
+        @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+        K = kernel(expand(T², splat(k)))
+        @test (K * expand(P, exp))[0.1] ≈ (K*expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+
+        x = axes(P,1)
+        @test (k.(x,x') * expand(P, exp))[0.1] ≈ (k.(x,x') *expand(T, exp))[0.1] ≈ 2.5521805183347417
+
+        y = Inclusion(2..3)
+        @test (k.(y, x') * expand(P, exp))[2.1]  ≈ (k.(y, x') * expand(T, exp))[2.1] ≈ 13.80651898993351
+    end
+
+    @testset "Extension" begin
+        Ex = _BandedBlockBandedMatrix(Ones{Int}((oneto(1),unitblocks(oneto(∞)))), blockedrange(oneto(∞)), (0,0), (0,0))
+        # Ey = _BlockBandedMatrix(mortar(OneElement.(oneto(∞), oneto(∞))), oneto(∞), unitblocks(oneto(∞)), (0,0))
+
+        T = ChebyshevT()
+        T² = RectPolynomial(Fill(T,2))
+        u = T² * (Ex * transform(T, exp))
+        @test u[SVector(0.1,0.2)] ≈ exp(0.1)
+    end
+end
\ No newline at end of file
diff --git a/test/test_rect.jl b/test/test_rect.jl
index ed3490b..852594b 100644
--- a/test/test_rect.jl
+++ b/test/test_rect.jl
@@ -108,6 +108,24 @@ using ContinuumArrays: plotgridvalues
         @test (P²[:,Block.(1:100)] \ f) ≈ f.args[2][Block.(1:100)]
     end
 
+    @testset "Weak Laplacian" begin
+        W = Weighted(Jacobi(1,1))
+        P = Legendre()
+        W² = RectPolynomial(Fill(W, 2))
+        P² = RectPolynomial(Fill(P, 2))
+        𝐱 = axes(P²,1)
+        D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
+        Δ = -((D_x * W²)'*(D_x * W²) + (D_y * W²)'*(D_y * W²))
+
+        f = expand(P² , 𝐱 -> ((x,y) = 𝐱; x^2 + y^2 - 2))
+
+        KR = Block.(Base.OneTo(100))
+        @time 𝐜 = Δ[KR,KR] \ (W²'*f)[KR];
+        @test W²[SVector(0.1,0.2),KR]'*𝐜 ≈ (1-0.1^2)*(1-0.2^2)/2 
+
+        @test \(Δ, (W²'*f); tolerance=1E-15) ≈ [0.5; zeros(∞)]
+    end
+    
     @testset "Show" begin
         @test stringmime("text/plain", KronPolynomial(Legendre(), Chebyshev())) == "Legendre() ⊗ ChebyshevT()"
         @test stringmime("text/plain", KronPolynomial(Legendre(), Chebyshev(), Jacobi(1,1))) == "Legendre() ⊗ ChebyshevT() ⊗ Jacobi(1.0, 1.0)"