--- a/EquidistantGrid.jl	Mon Feb 25 10:19:11 2019 -0800
+++ b/EquidistantGrid.jl	Mon Apr 15 16:15:04 2019 +0200
@@ -52,5 +52,5 @@
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
     @assert dim<=dimension(grid)
     @assert dim>0
-    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
+    points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim]))
 end

--- a/diffOp.jl	Mon Feb 25 10:19:11 2019 -0800
+++ b/diffOp.jl	Mon Apr 15 16:15:04 2019 +0200
@@ -82,10 +82,11 @@
 using TiledIteration
 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
     ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
-    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2)) # TBD: Is this the right way, the right size?
+    # TODO: Pass Tilesize to function
+    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
         for j ∈ tileaxs[2], i ∈ tileaxs[1]
             I = ri[i,j]
-            u[i,j] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
+            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
         end
     end
     return nothing

--- a/index.jl	Mon Feb 25 10:19:11 2019 -0800
+++ b/index.jl	Mon Apr 15 16:15:04 2019 +0200
@@ -39,6 +39,10 @@
 
 IndexTuple(t::Vararg{Tuple{T, DataType}}) where T<:Integer = Index.(t)
 
+# TODO: Use the values of the region structs, e.g. Lower(), for the region parameter instead of the types.
+# For example the following works:
+#   (Lower(),Upper()) isa NTuple{2, Region} -> true
+#   typeof((Lower(),Upper()))               -> Tuple{Lower,Upper}
 function regionindices(gridsize::NTuple{Dim,Integer}, closuresize::Integer, region::NTuple{Dim,DataType}) where Dim
     return regionindices(gridsize, ntuple(x->closuresize,Dim), region)
 end

--- a/plotDerivative2d.jl	Mon Feb 25 10:19:11 2019 -0800
+++ b/plotDerivative2d.jl	Mon Apr 15 16:15:04 2019 +0200
@@ -1,16 +1,16 @@
+include("sbpPlot.jl")
+
 g = sbp.Grid.EquidistantGrid((100,75), (0.0, 0.0), (2pi, 3/2*pi))
 op = sbp.readOperator("d2_4th.txt","h_4th.txt")
 Laplace = sbp.Laplace(g, 1.0, op)
 
 init(x,y) = sin(x) + cos(y)
 v = sbp.Grid.evalOn(g,init)
-
-u = zeros(length(v))
+u = zero(v)
 
 sbp.apply!(Laplace,u,v)
 
 #@show u
 #@show u'*u
 
-sbp.Grid.plotgridfunction(g,u)
-
+plotgridfunction(g,u)

--- a/sbpPlot.jl	Mon Feb 25 10:19:11 2019 -0800
+++ b/sbpPlot.jl	Mon Apr 15 16:15:04 2019 +0200
@@ -1,17 +1,12 @@
-module sbpPlot
-using PyPlot, PyCall
-
-function plotgridfunction(grid::EquidistantGrid, gridfunction)
-    if dimension(grid) == 1
-        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
-    elseif dimension(grid) == 2
-        mx = grid.size[1]
-        my = grid.size[2]
-        X = repeat(pointsalongdim(grid,1),1,my)
-        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
-        plot_surface(X,Y,reshape(gridfunction,mx,my));
+include("sbp.jl")
+using Makie
+import .sbp.Grid
+function plotgridfunction(grid::sbp.Grid.EquidistantGrid, gridfunction::AbstractArray)
+    if sbp.Grid.dimension(grid) == 1
+        plot(sbp.Grid.pointsalongdim(grid,1), gridfunction)
+    elseif sbp.Grid.dimension(grid) == 2
+        scene = surface(sbp.Grid.pointsalongdim(grid,1),sbp.Grid.pointsalongdim(grid,2), gridfunction)
     else
         error(string("Plot not implemented for dimension ", string(dimension(grid))))
     end
 end
-end