--- a/diffOp.jl	Fri Jan 25 15:10:41 2019 +0100
+++ b/diffOp.jl	Fri Jan 25 16:47:51 2019 +0100
@@ -33,10 +33,10 @@
 end
 
 # Differential operator for a*d^2/dx^2
-struct Laplace{Dim, T<:Real} <: DiffOp
+struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp
     grid::Grid.EquidistantGrid{Dim,T}
     a::T
-    op::D2{Float64}
+    op::D2{Float64,N,M,K}
 end
 
 # u = L*v
@@ -63,7 +63,7 @@
         @inbounds indices = uview(li,i,:)
         @inbounds apply!(L.op, temp, uview(v, indices), h[2])
         for i ∈ eachindex(indices)
-            @inbounds u[indices[i]] += temp[i]
+            @inbounds u[indices[i]] = temp[i]
         end
     end