--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/TimeStepper.jl	Mon Dec 17 14:50:23 2018 +0100
@@ -0,0 +1,60 @@
+abstract type TimeStepper end
+
+# Returns v and t
+function getState(ts::TimeStepper)
+	error("not implemented")
+end
+
+
+function step(ts::TimeStepper)
+	error("not implemented")
+end
+
+function stepN(ts::TimeStepper,N::Int)
+	for i ∈ 1:N
+		ts.step()
+	end
+end
+
+function stepTo(ts::TimeStepper)
+	error("Not yet implemented")
+end
+
+function evolve(ts::TimeStepper)
+	error("Not yet implemented")
+end
+
+
+mutable struct Rk4 <: TimeStepper
+	F::Function
+	k::Real
+	v::Vector
+	t::Real
+	n::UInt
+
+	function Rk4(F::Function,k::Real,v0::Vector,t0::Real)
+		# TODO: Check that F has two inputs and one output
+		v = v0
+		t = t0
+		n = 0
+		return new(F,k,v,t,n)
+	end
+end
+
+function getState(ts::Rk4)
+	return ts.t, ts.v
+end
+
+function step(ts::Rk4)
+    k1 = ts.F(ts.v,ts.t)
+	k2 = ts.F(ts.v+0.5*ts.k*k1,ts.t+0.5*ts.k)
+	k3 = ts.F(ts.v+0.5*ts.k*k2,ts.t+0.5*ts.k)
+    k4 = ts.F(ts.v+    ts.k*k3,ts.t+    ts.k)
+    ts.v  = ts.v + (1/6)*(k1+2*(k2+k3)+k4)*ts.k
+
+	ts.n = ts.n + 1
+	ts.t = ts.t + ts.k
+
+	return nothing
+end
+

--- a/diffOp.jl	Mon Dec 17 14:45:10 2018 +0100
+++ b/diffOp.jl	Mon Dec 17 14:50:23 2018 +0100
@@ -0,0 +1,51 @@
+abstract type DiffOp end
+
+function apply(D::DiffOp, v::AbstractVector)
+    error("not implemented")
+end
+
+function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real
+    error("not implemented")
+end
+
+function matrixRepresentation(D::DiffOp)
+    error("not implemented")
+end
+
+function boundaryCondition(D::DiffOp)
+    error("not implemented")
+end
+
+function interface(Du::DiffOp, Dv::DiffOp, b::BoundaryID; type)
+    error("not implemented")
+end
+
+
+# Differential operator for a*d^2/dx^2
+struct Laplace1D <: DiffOp
+    grid
+    a
+    op
+end
+
+# u = L*v
+function apply(L::Laplace1D, u::AbstractVector, v::AbstractVector)
+    N = closureSize(L.op)
+    M = length(v)
+
+    h = scaling(L.grid)
+
+    for i ∈ 1:N
+        u[i] = apply(L.op.closureStencils[i], v, i)/h^2
+    end
+
+    for i ∈ N+1:M-N
+        u[i] = apply(L.op.innerStencil, i)/h^2
+    end
+
+    for i ∈ M:-1:M-N+1
+        u[i] = apply(flip(L.op.closureStencils[M-i+1]), v, i)/h^2
+    end
+
+    return nothing
+end

--- a/sbp.jl	Mon Dec 17 14:45:10 2018 +0100
+++ b/sbp.jl	Mon Dec 17 14:50:23 2018 +0100
@@ -1,5 +1,7 @@
 module sbp
 include("grid.jl")
 include("diffOp.jl")
+include("stencil.jl")
 include("sbpD2.jl")
+include("TimeStepper.jl")
 end  # module

--- a/sbpD2.jl	Mon Dec 17 14:45:10 2018 +0100
+++ b/sbpD2.jl	Mon Dec 17 14:50:23 2018 +0100
@@ -1,1 +1,12 @@
 
+struct D2{T}
+    quadratureClosure::Vector{T}
+    innerStencil::Stencil
+    closureStencils::Vector{Stencil} # TBD: Should this be a tuple?
+    eClosure::Vector{T}
+    dClosure::Vector{T}
+end
+
+function closureSize(D::D2)::Int
+    return length(quadratureClosure)
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stencil.jl	Mon Dec 17 14:50:23 2018 +0100
@@ -0,0 +1,33 @@
+struct Stencil
+    range::NTuple{2,Int}
+    weights::Vector # TBD: Should this be a tuple?
+    function Stencil(range, weights)
+        width = range[2]-range[1]+1
+        if width != length(weights)
+            error("The width and the number of weights must be the same")
+        end
+        new(range, weights)
+    end
+end
+
+function flip(s::Stencil)
+    range = (-s.range[2], -s.range[1])
+    s = Stencil(range, s.weights[end:-1:1])
+end
+
+# Provides index into the Stencil based on offset for the root element
+function Base.getindex(s::Stencil, i::Int)
+    if s.range[1] <= i <= s.range[2]
+        return s.weights[1 + i - s.range[1]]
+    else
+        return 0
+    end
+end
+
+function apply(s::Stencil, v::AbstractVector, i::Int)
+    w = zero(v[0])
+    for j ∈ i+(s.range[1]:s.range[2])
+        w += v[j]
+    end
+    return w
+end