--- a/d2_2nd.txt	Thu Jan 10 15:49:44 2019 +0100
+++ b/d2_2nd.txt	Thu Jan 10 15:53:44 2019 +0100
@@ -7,7 +7,7 @@
 1 -2 1
 
 e
-1 0
+1
 
 d1
 -3/2 2 -1/2

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/d2_4th.txt	Thu Jan 10 15:53:44 2019 +0100
@@ -0,0 +1,16 @@
+# D2 order 4
+
+boundary_stencils
+     2    -5       4      -1     0     0
+     1    -2       1       0     0     0
+ -4/43 59/43 -110/43   59/43 -4/43     0
+ -1/49     0   59/49 -118/49 64/49 -4/49
+
+inner_stencil
+-1/12     4/3  -5/2   4/3 -1/12
+
+e
+1
+
+d1
+-11/6 3 -3/2 1/3
\ No newline at end of file

--- a/grid.jl	Thu Jan 10 15:49:44 2019 +0100
+++ b/grid.jl	Thu Jan 10 15:53:44 2019 +0100
@@ -1,5 +1,5 @@
 module grid
-using Plots
+using PyPlot
 
 abstract type Grid end
 
@@ -41,7 +41,7 @@
         return new(nPointsPerDim, lims)
     end
     # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
-    function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Integer})
+    function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
         return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
     end
 
@@ -94,26 +94,42 @@
     end
     dx̄ = Tuple(dx̄)
 
-    nPoints = numberOfPoints(grid)
-    points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, nPoints)
+    points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid))
     # Compute the points based on their Cartesian indices and the signed
     # grid spacings
     cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ 1:nPoints
+    for i ∈ 1:numberOfPoints(grid)
         ci = Tuple(cartesianIndices[i]) .-1
         points[i] = grid.limits[1] .+ dx̄.*ci
     end
+    # TBD: Keep? this? How do we want to represent points in 1D?
+    if numberOfDimensions(grid) == 1
+        points = broadcast(x -> x[1], points)
+    end
     return points
 end
 
-function plotOnGrid(grid::EquidistantGrid,v::Vector)
-    dim = numberOfDimensions(grid)
-    x = points(grid)
+function pointsalongdim(grid::EquidistantGrid, dim::Integer)
+    @assert dim<=numberOfDimensions(grid)
+    @assert dim>0
+    points = range(grid.limits[1][dim],stop=grid.limits[2][dim],length=grid.numberOfPointsPerDim[dim])
+end
 
-    if dim ==1
-        plot(x,v)
+function plotgridfunction(grid::EquidistantGrid, gridfunction)
+    if numberOfDimensions(grid) == 1
+        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
+    elseif numberOfDimensions(grid) == 2
+        x = pointsalongdim(grid,1)
+        X = repeat(x,1,grid.numberOfPointsPerDim[2])
+        y = pointsalongdim(grid,2)
+        Y = repeat(y,1,grid.numberOfPointsPerDim[1])'
+    elseif numberOfDimensions(grid) == 3
+        x = pointsalongdim(grid,1)
+        X = repeat(x,1,grid.numberOfPointsPerDim[2])
+        y = pointsalongdim(grid,2)
+        Y = repeat(y,1,grid.numberOfPointsPerDim[1])'
     else
-        error(string("Plot not implemented for dim =", string(dim)))
+        error(string("Plot not implemented for dimension =", string(dim)))
     end
 end
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/h_4th.txt	Thu Jan 10 15:53:44 2019 +0100
@@ -0,0 +1,7 @@
+# H order 4
+
+closure
+17/48 59/48 43/48 49/48
+
+inner_stencil
+1
\ No newline at end of file

--- a/sbpD2.jl	Thu Jan 10 15:49:44 2019 +0100
+++ b/sbpD2.jl	Thu Jan 10 15:53:44 2019 +0100
@@ -39,31 +39,33 @@
     h = readSectionedFile(Hfn)
 
     # Create inner stencil
-    innerStencilWeights = stringToVector(Float64, d["inner_stencil"])
+    innerStencilWeights = stringToVector(Float64, d["inner_stencil"][1])
     width = length(innerStencilWeights)
-    r = (-width//2, width//2)
+    r = (-div(width,2), div(width,2))
+
     innerStencil = Stencil(r, innerStencilWeights)
 
     # Create boundary stencils
     boundarySize = length(d["boundary_stencils"])
     closureStencils = Vector{Stencil}()
+
     for i ∈ 1:boundarySize
         stencilWeights = stringToVector(Float64, d["boundary_stencils"][i])
-
+        width = length(stencilWeights)
+        r = (1-i,width-i)
+        push!(closureStencils,Stencil(r, stencilWeights))
     end
 
     d2 = D2(
-        stringToVector(Float64, h["closure"]),
+        stringToVector(Float64, h["closure"][1]),
         innerStencil,
         closureStencils,
-        stringToVector(Float64, d["e"]),
-        stringToVector(Float64, d["d1"]),
+        stringToVector(Float64, d["e"][1]),
+        stringToVector(Float64, d["d1"][1]),
         even
     )
 
-    # Return d2!
-
-    return nothing
+    return d2
 end
 
 
@@ -91,6 +93,6 @@
     return sections
 end
 
-function stringToVector(T::DataType, s::String; delimiter = " ")
-    return parse(T, split(s, delimiter))
+function stringToVector(T::DataType, s::String)
+    return T.(eval.(Meta.parse.(split(s))))
 end