diff --git a/examples/tdpmchwt_cq.jl b/examples/tdpmchwt_cq.jl
new file mode 100644
index 00000000..e19e8746
--- /dev/null
+++ b/examples/tdpmchwt_cq.jl
@@ -0,0 +1,81 @@
+using BEAST, LinearAlgebra, CompScienceMeshes, ConvolutionOperators, Printf
+
+# Physical parameters
+ϵ, μ = 1.0, 1.0
+ϵ′, μ′ = 1.0, 1.0
+η, η′ = √(μ/ϵ), √(μ′/ϵ′)
+
+T = BEAST.LaplaceDomainOperator(s::ComplexF64 -> MWSingleLayer3D(s*√(ϵ*μ), -s*√(ϵ*μ), ComplexF64(1/√(ϵ*μ))/s))
+K = BEAST.LaplaceDomainOperator(s::ComplexF64 -> MWDoubleLayer3D(s*√(ϵ*μ)))
+T′ = BEAST.LaplaceDomainOperator(s::ComplexF64 -> MWSingleLayer3D(s*√(ϵ′*μ′), -s*√(ϵ′*μ′), ComplexF64(1/√(ϵ′*μ′))/s))
+K′ = BEAST.LaplaceDomainOperator(s::ComplexF64 -> MWDoubleLayer3D(s*√(ϵ′*μ′)))
+
+Γ = meshsphere(1.0, 0.3)
+X = raviartthomas(Γ)
+Y = buffachristiansen(Γ)
+
+@hilbertspace k l
+@hilbertspace j m
+
+Nt, Δt = 200, 10.0
+(A, b, c) = butcher_tableau_radau_2stages()
+CQ = StagedTimeStep(Δt, Nt, c, A, b, 30, 1.0001)
+
+duration = 40 * Δt                                       
+delay = 60 * Δt                                        
+amplitude = 1.0
+gaussian = creategaussian(duration, delay, amplitude)
+fgaussian = fouriertransform(gaussian)
+polarisation, direction = x̂, ẑ
+E = planewave(polarisation, direction, gaussian, 1.0)
+H = direction × E
+
+BEAST.@defaultquadstrat (T, X⊗CQ, X⊗CQ) BEAST.DoubleNumWiltonSauterQStrat(6, 7, 6, 7, 7, 7, 7, 7)
+BEAST.@defaultquadstrat (T′, X⊗CQ, X⊗CQ) BEAST.DoubleNumWiltonSauterQStrat(6, 7, 6, 7, 7, 7, 7, 7)
+BEAST.@defaultquadstrat (K, X⊗CQ, X⊗CQ) BEAST.DoubleNumWiltonSauterQStrat(6, 7, 6, 7, 7, 7, 7, 7)
+BEAST.@defaultquadstrat (K′, X⊗CQ, X⊗CQ) BEAST.DoubleNumWiltonSauterQStrat(6, 7, 6, 7, 7, 7, 7, 7)
+
+pmchwt = @discretise(η*T[k,j] + η′*T′[k,j] + K[l,j] + K′[l,j] 
+                        - K[k,m] - K′[k,m] + 1/η*T[l,m] + 1/η′*T′[l,m] == E[k] + H[l],
+                        k∈X⊗CQ, l∈X⊗CQ, j∈X⊗CQ, m∈X⊗CQ)
+je = solve(pmchwt)
+
+
+# Exact solution
+N = NCross()
+Nyx = assemble(N, Y, X)
+iNyx = inv(Matrix(Nyx))
+Zr = zeros(size(iNyx))
+
+δ = timebasisdelta(Δt, Nt)	                			                       
+linform = @discretise(H[k] + E[l], k∈Y⊗δ, l∈Y⊗δ)
+rhs = BEAST.assemble(linform.linform, linform.test_space_dict)
+j_ext = [iNyx Zr; Zr iNyx] * rhs
+
+
+# Plot the numerical solutions
+using Plots
+x = 1e-2 * Δt/3 * [1:1:Nt;]
+
+plt = Plots.plot(
+    size = (600, 400),
+    grid = false,
+    xscale = :identity, 
+    # xlims = (0, 1620),
+    # xticks = [400, 800, 1200, 1600],
+    xtickfont = Plots.font(10, "Computer Modern"), 
+    yscale = :log10, 
+    ylims = (1e-15, 1e0), 
+    yticks = [1e-20, 1e-15, 1e-10, 1e-5, 1e0, 1e5],
+    ytickfont = Plots.font(10, "Computer Modern"),
+    xlabel = "Time ({\\mu}s)",
+    ylabel = "Current density intensity (A/m)",
+    titlefont = Plots.font(10, "Computer Modern"),
+    guidefont = Plots.font(11, "Computer Modern"),
+    colorbar_titlefont = Plots.font(10, "Computer Modern"),
+    legendfont = Plots.font(11, "Computer Modern"),
+    legend = :topright,
+    dpi = 300)
+    
+Plots.plot!(x, abs.(j_ext[1, :]), label="Exact solution", linecolor=1, lw=1.2)
+Plots.plot!(x, abs.(je[1, :]), label="TD-PMCHWT", linecolor=2, lw=1.2)
\ No newline at end of file
diff --git a/src/timedomain/rkcq.jl b/src/timedomain/rkcq.jl
index 13ae3b24..6dcfd3ec 100644
--- a/src/timedomain/rkcq.jl
+++ b/src/timedomain/rkcq.jl
@@ -106,3 +106,75 @@ function assemble(rkcq :: RungeKuttaConvolutionQuadrature,
 	return ZC
 
 end
+
+
+struct LaplaceDomainOperator
+	coeff::ComplexF64
+	kernel::Any
+end
+
+scalartype(op::LaplaceDomainOperator) = ComplexF64
+
+*(a::Number, op::LaplaceDomainOperator) = LaplaceDomainOperator(a*op.coeff, op.kernel)
+
+LaplaceDomainOperator(kernel) = LaplaceDomainOperator(Complex(1.0), kernel)
+
+defaultquadstrat(op::LaplaceDomainOperator, tfs::SpaceTimeBasis, bfs::SpaceTimeBasis) = DoubleNumWiltonSauterQStrat(2,3,6,7,5,5,4,3)
+
+function assemble(op::LaplaceDomainOperator, testfns::SpaceTimeBasis, trialfns::SpaceTimeBasis; quadstrat=defaultquadstrat(op, testfns, trialfns))
+	laplaceKernel = op.kernel
+
+	A = testfns.time.A
+	b = testfns.time.b
+	Δt = testfns.time.Δt
+	Q = testfns.time.zTransformedTermCount
+	rho = testfns.time.contourRadius
+	p = length(b) # stage count
+
+	test_spatial_basis  = testfns.space
+	trial_spatial_basis = trialfns.space
+
+	# Compute the Z transformed sequence.
+	# Assume that the operator applied on the conjugate of s is the same as the
+	# conjugate of the operator applied on s,
+	# so that only half of the values are computed
+	Qmax = Q>>1+1
+	M = numfunctions(test_spatial_basis)
+	N = numfunctions(trial_spatial_basis)
+	Tz = ComplexF64
+	Zz = Vector{Array{Tz,2}}(undef,Qmax)
+	blocksEigenvalues = Vector{Array{Tz,2}}(undef,p)
+	tmpDiag = Vector{Tz}(undef,p)
+
+	for q = 0:Qmax-1
+		# Build a temporary matrix for each eigenvalue
+		s = laplace_to_z(rho, q, Q, Δt, A, b)
+		sFactorized = diagonalizedmatrix(s)
+		for (i,sD) in enumerate(sFactorized.D)
+			blocksEigenvalues[i] = op.coeff * assemble(laplaceKernel(sD), test_spatial_basis, trial_spatial_basis, quadstrat=quadstrat)
+		end
+
+		# Compute the Z transformed matrix by block
+		Zz[q+1] = zeros(Tz, M*p, N*p)
+		for m = 1:M
+			for n = 1:N
+				for i = 1:p
+					tmpDiag[i] = blocksEigenvalues[i][m,n]
+				end
+				D = SVector{p,Tz}(tmpDiag);
+				Zz[q+1][(m-1)*p.+(1:p),(n-1)*p.+(1:p)] = sFactorized.H * diagm(D) * sFactorized.invH
+			end
+		end
+	end
+
+	# return the inverse Z transform
+	kmax = Q
+	T = real(Tz)
+	Z = zeros(T, M*p, N*p, kmax)
+	for q = 0:kmax-1
+		Z[:,:,q+1] = real_inverse_z_transform(q, rho, Q, Zz)
+	end
+	ZC = ConvolutionOperators.DenseConvOp(Z)
+
+	return ZC
+end
\ No newline at end of file
diff --git a/src/timedomain/tdtimeops.jl b/src/timedomain/tdtimeops.jl
index 61bf39d7..cdb3e99e 100644
--- a/src/timedomain/tdtimeops.jl
+++ b/src/timedomain/tdtimeops.jl
@@ -72,7 +72,7 @@ function allocatestorage(op::TensorOperator, test_functions, trial_functions,
     tail = zeros(T, M, N)
 
     Nt = numfunctions(temporalbasis(trial_functions))
-    Z = ConvolutionOperators.ConvOp(data, K0, K1, tail, Nt)
+    Z = ConvolutionOperators.ConvOp(data, K0, K1, tail, bandwidth)
     function store1(v,m,n,k)
 		if Z.k0[m,n] ≤ k ≤ Z.k1[m,n]
 			Z.data[k - Z.k0[m,n] + 1,m,n] += v
diff --git a/test/test_td_tensoroperator.jl b/test/test_td_tensoroperator.jl
index db480dac..3bddb66d 100644
--- a/test/test_td_tensoroperator.jl
+++ b/test/test_td_tensoroperator.jl
@@ -28,7 +28,7 @@ using Test
     op = Id ⊗ Id
     Z = assemble(op, V, U)
     A = BEAST.ConvolutionOperators.ConvOpAsArray(Z)
-    @test size(A) == (1,1,10)
+    @test size(A) == (1,1,2)
     for i in 2:10
         @test A[1,1,i] ≈ 0 atol=sqrt(eps(eltype(A)))
     end