Better README and added docstrings

README.md

@@ -2,7 +2,7 @@
 
 [![codecov](https://codecov.io/github/QuEraComputing/DormandPrince.jl/graph/badge.svg?token=qYZ4V7m0JY)](https://codecov.io/github/QuEraComputing/DormandPrince.jl)
 
-Julia-native implementation of the Dormand-Prince 5th order solver
+Julia-native implementation of the Dormand-Prince 5th and 8th order solvers
 
 ## Installation
 
@@ -22,6 +22,66 @@ DormandPrince is a  
 pkg> add DormandPrince
 ```
 
+## Usage 
+
+### Single End-Time
+
+```julia
+julia> using DormandPrince
+
+# define your ODE, in this case, dy/dx = 0.85y
+julia> function fcn(x, y, f)
+            f[1] = 0.85y[1]
+        end
+
+# Create a solver object. We will use the 5th order solver and
+# start integrating at t = 0.0 with initial value of 19.0
+julia> solver = DP5Solver(fcn, 0.0, [19.0])
+
+# Begin integration up to t = 2.1
+julia> integrate!(solver, 2.1)
+
+# get_current_state will return the answer
+julia> get_current_state(solver)
+```
+
+Both the `DP5Solver` and `DP8Solver`'s are stateful allowing for memory-efficient integration to future end times from the last integrated end point (e.g. if you chose t = 1.0 as your endpoint you can call `integrate!` again with t=2.0 and it will "carry forward" the work starting from t = 1.0 instead of requiring you to set things up all over again).
+
+### Multiple End-Times
+
+```julia
+julia> using DormandPrince
+
+# Define ODE
+julia> function fcn(x, y, f)
+            f[1] = 0.85y[1]
+        end
+
+# Define times of interest to analyze/perform actions on the solution
+julia> times = [1.0, 1.1, 1.9, 2.4]
+
+# Create the solver object, with integration starting from t = 0.0 and initial value of 19.0
+julia> solver = DP5Solver(fcn, 0.0, [19.0])
+
+# Empty array to store intermediate values
+julia> intermediate_values = []
+
+# Use callback to store intermediate values. The solver object will also be mutated to store the solution 
+# found at the last time point.
+julia> integrate!(solver, times) do time, val
+            push!(intermediate_values, val[])
+        end
+
+```
+
+You may also create a `SolverIterator` that can then be iterated over. Note that this will require a fresh solver
+
+```julia
+for (time, value) in SolverIterator(solver, times)
+    println("Time: ", time, " ", "Value: ", value[])
+end
+```
+
 ## License
 
 Apache License 2.0

src/interface.jl

@@ -1,4 +1,28 @@
+"""
+    struct SolverIterator{T <: Real}
+    SolverIterator(solver, times)
+
+Given a solver and a vector of times, 
+this iterator will return the state of the solver at each time.
 
+The solver will be mutated to hold the solution from the last time given.
+
+# Examples
+
+```julia
+julia> solver = DP5Solver(fcn, 0.0, [0.0])
+
+julia> times = [1.0, 2.0, 3.0]
+
+julia> intermediate_vals = []
+
+julia> for (time, value) in SolverIterator(solver, times)
+            push!(intermediate_values, value[])
+        end
+
+```
+
+"""
 struct SolverIterator{T <: Real}
     solver::AbstractDPSolver{T}
     times::AbstractVector{T}
@@ -32,12 +56,69 @@ end
 get_current_state(::AbstractDPSolver) = error("not implemented")
 integrate_core!(::AbstractDPSolver{T}, ::T) where T = error("not implemented")
 
+"""
+    integrate!(solver, time)
+
+Integrate the ODE problem that is part of solver to the end time `time`.
+`solver` can be a `DP5Solver` or a `DP8Solver` type.
+
+The `solver` is mutated to hold the solution and `solver` state at the time `time`, allowing for 
+efficient integration to future end times of interest through subsequent calls to `integrate!` with 
+later times.
+
+# Examples
+
+```julia
+julia> function fcn(x, y, f)
+            f[1] = 0.85y[1]
+        end
+
+julia> solver = DP5Solver(fcn, 0.0, [19.0])
+
+# Integrate from initial time of 0.0 to 1.0
+julia> integrate!(solver, 1.0)
+
+# integrate from the last time of 1.0 to 2.0
+julia> integrate!(solver, 2.0)
+
+# Solver object now holds solution and state at 2.0
+julia> get_current_state(solver)
+```
+"""
 function integrate!(solver::AbstractDPSolver{T}, time::T) where T <: Real
     report = integrate_core!(solver, time)
     if report.idid != COMPUTATION_SUCCESSFUL
         throw(DPException(report))
     end
 end
+
+"""
+    integrate!(callback, solver, times)
+
+Integrate the ODE problem that is part of `solver` to the end times `times` and apply the callback on 
+the time and solution at each time. 
+`solver` can be a `DP5Solver` or a `DP8Solver` type.
+
+At the end of all the times the solver holds the solution and solver state at the last time in `times`.
+
+# Examples
+
+```julia
+julia> function fcn(x, y, f)
+            f[1] = 0.85y[1]
+        end
+
+julia> times = [1.0, 1.1, 1.9, 2.4]
+
+julia> solver = DP5Solver(fcn, 0.0, [19.0])
+
+julia> intermediate_values = []
+
+julia> integrate!(solver, times) do time, val
+            push!(intermediate_values, val[])
+        end
+```
+"""
 function integrate!(callback, solver::AbstractDPSolver{T}, times::AbstractVector{T}; sort_times::Bool = true) where {T <: Real}
     times = sort_times ? sort(collect(times)) : times