Solvned project from 2020

Note that Go didn't have generics back then.

32 files changed, 3091 insertions, 0 deletions

diff --git a/COPYING.LESSER b/COPYING.LESSER
new file mode 100644
index 0000000..5357f69
--- /dev/null
+++ b/COPYING.LESSER
@@ -0,0 +1,166 @@
+                   GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.
+
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..3e32a67
--- /dev/null
+++ b/README.md
@@ -0,0 +1,10 @@
+# Solvned: A numerical ODE framework for Go
+
+Solvned is a framework for solving ordinary differential equations in Go. It
+provides a variety of methods to solve these equations, from the simple Euler
+method to the sophisticated Runge-Kutta-Fehlberg method, as well as the means
+to define new methods and operate with them in a general way, with 
+interpolation and an event mechanism for iteration until a given point.
+
+Check out the visualization examples at the `examples` directory and the 
+[pkg.go.dev documentation](https://pkg.go.dev/github.com/JwanMan/mned).
diff --git a/caching.go b/caching.go
new file mode 100644
index 0000000..1d31b94
--- /dev/null
+++ b/caching.go
@@ -0,0 +1,301 @@
+// Caching of solution points.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+import "sort"
+
+// A CacheSolver is like a dynamic version of a DenseSolution. It stores
+// points of the solution of a problem in a conceptually maximal interval that
+// is evaluated lazily; that is, it starts with the solutions in a range that
+// only includes the initial point and, when a point is asked for a Time value
+// outside the range, the range is expanded to cover that point.
+type CacheSolver struct {
+	backward []Point // Invariant: Decreasing times, times before forward's.
+	forward  []Point // Invariant: Non-empty, increasing times.
+	backStep Stepper // Invariant: Decreasing times from end of backward.
+	forStep  Stepper // Invariant: Increasing times from end of forward.
+	evs      []Event
+	interp   Interpolator
+}
+
+// Create a CacheSolver to solve the given IVP with the given solving Method,
+// with the given interpolator to calculate points between steps and taking into
+// account the events provided.
+func CacheSolve(
+	m Method, ivp *IVP, interp Interpolator, evs ...Event,
+) CacheSolver {
+	return CacheSolver{
+		forward:  []Point{ivp.Start.Clone()},
+		backward: []Point{},
+		forStep:  m.Forward(ivp),
+		backStep: m.Backward(ivp),
+		evs:      evs,
+		interp:   interp,
+	}
+}
+
+// Get the time of the stored solution point with the lowest time.
+func (c *CacheSolver) Start() float64 {
+	if len(c.backward) > 0 {
+		return c.backward[len(c.backward)-1].Time
+	} else {
+		return c.forward[0].Time
+	}
+}
+
+// Get the time of the stored solution point with the greatest time.
+func (c *CacheSolver) End() float64 {
+	return c.forward[len(c.forward)-1].Time
+}
+
+func (c *CacheSolver) getActions(
+	current *Point, next *Point, vals []float64,
+) requiredActions {
+	var actions requiredActions = make([]requiredAction, 0)
+	for i := 0; i < len(vals); i++ {
+		newVal := c.evs[i].Cross(next)
+		if differentSign(vals[i], newVal) {
+			point := c.evs[i].FindPoint(c.interp, current, next)
+			actions = append(
+				actions,
+				requiredAction{index: i, point: point},
+			)
+		}
+		vals[i] = newVal
+	}
+	sort.Sort(actions)
+	return actions
+}
+
+func (c *CacheSolver) expandBackward(t float64) bool {
+	if c.backStep == nil {
+		return false
+	}
+
+	var current Point
+	if len(c.backward) == 0 {
+		current = c.forward[0]
+	} else {
+		current = c.backward[len(c.backward)-1]
+	}
+
+	vals := make([]float64, len(c.evs))
+	for i := 0; i < len(vals); i++ {
+		vals[i] = c.evs[i].Cross(&current)
+	}
+
+	for t < current.Time {
+		next, ok := c.backStep.Next()
+		if !ok {
+			// TODO Consider bisection for reasonable last point
+			c.backStep = nil
+			return false
+		}
+		actions := c.getActions(&current, next, vals)
+		for i := len(actions) - 1; i >= 0; i-- {
+			if !c.evs[actions[i].index].Action(&actions[i].point) {
+				return false
+			}
+		}
+		current = next.Clone()
+		c.backward = append(c.backward, current)
+	}
+	return true
+}
+
+func (c *CacheSolver) expandForward(t float64) bool {
+	if c.forStep == nil {
+		return false
+	}
+
+	current := c.forward[len(c.forward)-1]
+	vals := make([]float64, len(c.evs))
+	for i := 0; i < len(vals); i++ {
+		vals[i] = c.evs[i].Cross(&current)
+	}
+
+	for t < current.Time {
+		next, ok := c.forStep.Next()
+		if !ok {
+			// TODO Consider bisection for reasonable last point
+			c.forStep = nil
+			return false
+		}
+		actions := c.getActions(&current, next, vals)
+		for i := 0; i < len(actions); i++ {
+			if !c.evs[actions[i].index].Action(&actions[i].point) {
+				return false
+			}
+		}
+		current = next.Clone()
+		c.forward = append(c.forward, current)
+	}
+	return true
+}
+
+// Get the value of the solution for a given time.
+//
+// If the value is outside of the bounds of what's stored, values are added
+// until reaching the given time. If ok is false, the value is out of bounds for
+// the problem or some event and points have only been added as far as it was
+// possible.
+func (c *CacheSolver) Get(t float64) (x []float64, ok bool) {
+	if t < c.Start() {
+		if ok := c.expandBackward(t); !ok {
+			return nil, false
+		}
+	}
+	if t > c.End() {
+		if ok := c.expandForward(t); !ok {
+			return nil, false
+		}
+	}
+	if t >= c.forward[0].Time {
+		i := 1
+		for t > c.forward[i].Time {
+			i++
+		}
+		return c.interp.FindValue(
+			&c.forward[i-1], &c.forward[i], t,
+		), true
+	}
+
+	i := 0
+	for t > c.backward[i].Time {
+		i++
+	}
+	if i == 0 {
+		return c.interp.FindValue(
+			&c.backward[0], &c.forward[0], t,
+		), true
+	}
+	return c.interp.FindValue(&c.backward[i], &c.backward[i-1], t), true
+}
+
+// Get the value at `c.Start() - h` as in `c.Get` except that, if the underlying
+// backward Stepper is a ConfigurableStepper, its method `NextStep(-h)` is used.
+// The result value should *NOT* be mutated.
+//
+// This is most useful with fixed step methods to specify the step.
+func (c *CacheSolver) StepBackward(h float64) ([]float64, bool) {
+	if cs, ok := c.backStep.(ConfigurableStepper); ok {
+		if point, ok := cs.NextStep(-h); ok {
+			if h > 0 {
+				c.backward = append(c.backward, point.Clone())
+			}
+			return point.Value, true
+		} else {
+			return nil, false
+		}
+	}
+	return c.Get(c.Start() - h)
+}
+
+// Get the value at `c.End() + h` as in `c.Get` except that, if the underlying
+// forward Stepper is a ConfigurableStepper, its method `NextStep(h)` is used.
+// The result value should *NOT* be mutated.
+//
+// This is most useful with fixed step methods to specify the step.
+func (c *CacheSolver) StepForward(h float64) ([]float64, bool) {
+	if cs, ok := c.forStep.(ConfigurableStepper); ok {
+		if point, ok := cs.NextStep(h); ok {
+			if h > 0 {
+				c.forward = append(c.forward, point.Clone())
+			}
+			return point.Value, true
+		} else {
+			return nil, false
+		}
+	}
+	return c.Get(c.End() + h)
+}
+
+// Get the points calculated from the initial values to earlier times. The
+// points **MUST NOT** be modified.
+func (c *CacheSolver) BackwardPoints() []Point {
+	return c.backward
+}
+
+// Get the points calculated from the initial values to later times. The points
+// **MUST NOT** be modified.
+func (c *CacheSolver) ForwardPoints() []Point {
+	return c.forward
+}
+
+// Rearrange the stored solution points in a result such that the i-th point has
+// time `result[0][i]` and value `(result[1][i],...,result[size][i])`, where
+// `size` is the number of dimensions of a value. This is useful for plotting.
+func (c *CacheSolver) PointCoords() [][]float64 {
+	size := len(c.forward[0].Value)
+	result := make([][]float64, size+1)
+	backs := len(c.backward)
+	elems := backs + len(c.forward)
+
+	for i := 0; i <= size; i++ {
+		result[i] = make([]float64, elems)
+	}
+	for i := 0; i < backs; i++ {
+		result[0][i] = c.backward[backs-i-1].Time
+		for j := 0; j < size; j++ {
+			result[j+1][i] = c.backward[backs-i-1].Value[j]
+		}
+	}
+	for i := backs; i < elems; i++ {
+		result[0][i] = c.forward[i-backs].Time
+		for j := 0; j < size; j++ {
+			result[j+1][i] = c.forward[i-backs].Value[j]
+		}
+	}
+	return result
+}
+
+// Force the generation of values to the left until the domain ends or an event
+// tells the stepping to stop.
+func (c *CacheSolver) StepToBeginning() {
+	var init Point
+	if len(c.backward) > 0 {
+		init = c.backward[len(c.backward)-1]
+	} else {
+		init = c.forward[0]
+	}
+	StepUntil(
+		c.backStep,
+		init,
+		c.interp,
+		func(point *Point) {
+			c.backward = append(c.backward, point.Clone())
+		},
+		c.evs...,
+	)
+}
+
+// Force the generation of values to the right until the domain ends or an event
+// tells the stepping to stop.
+func (c *CacheSolver) StepToEnd() {
+	StepUntil(
+		c.forStep,
+		c.forward[len(c.forward)-1],
+		c.interp,
+		func(point *Point) {
+			c.forward = append(c.forward, point.Clone())
+		},
+		c.evs...,
+	)
+}
diff --git a/common.go b/common.go
new file mode 100644
index 0000000..9b41888
--- /dev/null
+++ b/common.go
@@ -0,0 +1,75 @@
+// Core data type definitions.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+/*
+This package implements several methods to compute approximations of
+solutions of initial value problems (IVPs) with first-order ODEs. These problems
+have the form
+```
+x'(t) = f(t, x(t)),
+x(t0) = x0,
+```
+where `t0 : R` and `x0 : R^n` are the **initial values**,
+`f : Ω ⊆ R x R^n -> R^n` is the **derivative function**, and `x : I ⊆ R -> R^n`
+is the unknown of the equation. Here `I` is an open interval of `R` containing
+`t0` and `Ω` is an open subset of `R x R^n` containing `(t0, x0)`.
+
+We refer to `t` as the **independent variable** and to `x(t)` as the `dependent
+variable`. We know that `x` is uniquely defined in a neighbourhood of `t0`.
+Given a point `t` in such a neighbourhood, we say that `(t,x(t))` is a point of
+the solution of the initial value problem.
+
+The results calculated will be approximations, and judgement is required to get
+approximations suitable for a particular problem; nevertheless, for the purpose
+of explanation, we will use the same terminology for the approximations than for
+the actual functions.
+*/
+
+// A Point is an element of the solution, given by the values of the independent
+// and dependent variables.
+type Point struct {
+	Time  float64   // The independent variable.
+	Value []float64 // The dependent variable.
+}
+
+// Deep-copy a point.
+func (p *Point) Clone() Point {
+	value := make([]float64, len(p.Value))
+	copy(value, p.Value)
+	return Point{Time: p.Time, Value: value}
+}
+
+// An IVP is an initial value problem given by a first-order ODE. It's given by
+// the Derivative function, which is a pure function, and the initial values.
+//
+// The second return value of Derivative indicates if the Point was in the
+// domain of the function. If it wasn't the first return value should not be
+// used. Ideally, the domain should be restricted to the connected component of
+// the initial value, as otherwise a solving method could "jump" to another
+// component and the results from there on would be invalid.
+type IVP struct {
+	Derivative func(Point) ([]float64, bool)
+	Start      Point
+}
+
+func (ivp *IVP) Clone() IVP {
+	return IVP{Derivative: ivp.Derivative, Start: ivp.Start.Clone()}
+}
diff --git a/dense.go b/dense.go
new file mode 100644
index 0000000..8bfe9c5
--- /dev/null
+++ b/dense.go
@@ -0,0 +1,174 @@
+// Dense solutions.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+// A DenseSolution is a structure with precalculated points of a solution within
+// a given interval. It allows getting the value of any point within that
+// interval via interpolation. The zero value is meaningless.
+type DenseSolution struct {
+	points []Point // This can't be empty. Points are ordered increasingly.
+	interp Interpolator
+}
+
+// Get the earliest time of the points in the solution.
+func (d *DenseSolution) Start() float64 {
+	return d.points[0].Time
+}
+
+// Get the latest time of the points in the solution.
+func (d *DenseSolution) End() float64 {
+	return d.points[len(d.points)-1].Time
+}
+
+func (d *DenseSolution) search(time float64, min int, max int) []float64 {
+	if min == max {
+		return d.points[min].Value
+	}
+	if max-min == 1 {
+		switch time {
+		case d.points[min].Time:
+			return d.points[min].Clone().Value
+		case d.points[max].Time:
+			return d.points[max].Clone().Value
+		}
+		return d.interp.FindValue(&d.points[min], &d.points[max], time)
+	}
+	mid := (max - min) / 2
+	midtime := d.points[mid].Time
+	switch {
+	case time == midtime:
+		return d.points[mid].Value
+	case time > midtime:
+		return d.search(time, mid, max)
+	default:
+		return d.search(time, min, mid)
+	}
+}
+
+// Get the value for a point in the solution with a given time between
+// `d.Start()` and `d.End()`. If the second return value is false, the point
+// was not in the interval and the first return value is meaningless.
+func (d *DenseSolution) Get(time float64) ([]float64, bool) {
+	if time < d.Start() {
+		return d.points[0].Value, false
+	}
+	if time > d.End() {
+		return d.points[len(d.points)-1].Value, false
+	}
+	return d.search(time, 0, len(d.points)-1), true
+}
+
+// Return the list of all explicitly-calculated points. This slice *MUST NOT*
+// be modified, nor any of the points it contains.
+func (d *DenseSolution) InnerPoints() []Point {
+	return d.points
+}
+
+// Like CacheSolver.PointCoords but for a DenseSolution.
+func (d *DenseSolution) PointCoords() [][]float64 {
+	size := len(d.points[0].Value)
+	result := make([][]float64, size+1)
+	for i := range result {
+		result[i] = make([]float64, len(d.points))
+	}
+	for i, p := range d.points {
+		result[0][i] = p.Time
+		for j, v := range p.Value {
+			result[j+1][i] = v
+		}
+	}
+	return result
+}
+
+// Solve an initial value problem in an interval `[start, end]` with the given
+// method. The points are calculated eagerly and, when a point is requested not
+// given by the method, interpolation is used with the given interpolator.
+//
+// If the second return value is false, the method didn't get to generate any
+// point and the first return value is meaningless. If it's true, the first
+// return value contains a solution for an interval that contains at most two
+// explicitly-stored points outside the interval requested (the two extremes)
+// but which might not contain the whole interval requested if the method
+// stopped.
+func DenseSolve(
+	method Method,
+	ivp *IVP,
+	start float64,
+	end float64,
+	interp Interpolator,
+) (DenseSolution, bool) {
+	forward := method.Forward(ivp)
+	backward := method.Backward(ivp)
+	points := make([]Point, 0)
+	current := &ivp.Start
+	pendingAddInitial := true
+	var prev Point
+	var ok bool
+
+	if end < current.Time { // Special case: end is on the left
+		for end < current.Time {
+			prev = current.Clone()
+			if current, ok = backward.Next(); !ok {
+				return DenseSolution{}, false
+			}
+		}
+		if current.Time < end {
+			points = append(points, prev)
+		}
+		points = append(points, current.Clone())
+		pendingAddInitial = false
+	}
+	for start < current.Time { // Move backwards until the start
+		if current, ok = backward.Next(); !ok {
+			break
+		}
+		points = append(points, current.Clone())
+	}
+	size := len(points) // Reverse backwards points
+	for i := 0; i < size/2; i++ {
+		points[i], points[size-i] = points[size-i], points[i]
+	}
+
+	current = &ivp.Start
+	if start > current.Time { // Special case: start is on the right
+		for start > current.Time {
+			prev = current.Clone()
+			if current, ok = forward.Next(); !ok {
+				return DenseSolution{}, false
+			}
+		}
+		if current.Time > start {
+			points = append(points, prev)
+		}
+		points = append(points, current.Clone())
+		pendingAddInitial = false
+	}
+	if pendingAddInitial {
+		points = append(points, ivp.Start)
+	}
+	for end > current.Time { // Move forwards until the end
+		if current, ok = forward.Next(); !ok {
+			break
+		}
+		points = append(points, current.Clone())
+	}
+
+	return DenseSolution{points: points, interp: interp}, true
+}
diff --git a/event.go b/event.go
new file mode 100644
index 0000000..c4a687b
--- /dev/null
+++ b/event.go
@@ -0,0 +1,95 @@
+// Events mechanism for early iteration ending.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+import "math"
+
+// An Event is a condition that can happen in a point in the solution of an
+// initial value problem and an associated action to take. The event happens
+// at the point of the solution where a given function changes its sign.
+//
+// The Cross function is a continuous function whose zeroes are the points
+// where the event happens. The Tolerance is the margin of error allowed; the
+// maximum absolute value of `Cross(p)` such that `p` is considered to be close
+// enough to an event to be passed to Action. The Action is to be called when
+// an occurrence of the event is found; it can use the point in some way and it
+// returns a boolean indicating whether the calculation should continue.
+type Event struct {
+	Cross     func(*Point) float64
+	Tolerance float64
+	Action    func(*Point) bool
+}
+
+// If e.Cross has a zero between p1 and p2, find such a zero or a point `p`
+// close enough to the zero that `|e.Cross(p)| < e.Tolerance`. To get the
+// intermediante points, the given interpolator is used.
+func (e *Event) FindPoint(i Interpolator, p1 *Point, p2 *Point) Point {
+	var min, max, mid Point
+	if p1.Time < p2.Time {
+		min, max = *p1, *p2
+	} else {
+		max, min = *p1, *p2
+	}
+	downwards := e.Cross(p1) > 0
+	for {
+		mid.Time = (min.Time + max.Time) / 2
+		mid.Value = i.FindValue(p1, p2, mid.Time)
+		value := e.Cross(&mid)
+		if math.Abs(value) < e.Tolerance {
+			return mid
+		}
+		if downwards {
+			if value > 0 {
+				min = mid
+			} else {
+				max = mid
+			}
+		} else {
+			if value > 0 {
+				max = mid
+			} else {
+				min = mid
+			}
+		}
+	}
+}
+
+type requiredAction struct {
+	index int
+	point Point
+}
+
+type requiredActions []requiredAction
+
+func (r requiredActions) Len() int {
+	return len(r)
+}
+
+func (r requiredActions) Less(i, j int) bool {
+	return r[i].point.Time < r[j].point.Time
+}
+
+func (r requiredActions) Swap(i, j int) {
+	r[i], r[j] = r[j], r[i]
+}
+
+func differentSign(f1 float64, f2 float64) bool {
+	return (f1 <= 0 && f2 >= 0) || (f1 >= 0 && f2 <= 0)
+}
diff --git a/examples/arenstorf.go b/examples/arenstorf.go
new file mode 100644
index 0000000..6506a6a
--- /dev/null
+++ b/examples/arenstorf.go
@@ -0,0 +1,74 @@
+// Example visualization of the Arenstorf orbit.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"github.com/Arafatk/glot"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/ivp"
+	"github.com/JwanMan/mned/method"
+	"os"
+)
+
+func main() {
+	var backed bool = false
+	problem := ivp.ArenstorfOrbit()
+	solution := mned.CacheSolve(
+		method.RKFehlberg(0.01, 0.01, 1e-7, 0.1),
+		&problem,
+		mned.HermiteForIVP(&problem),
+		mned.Event{
+			Cross: func(p *mned.Point) float64 {
+				return p.Value[1]
+			},
+			Tolerance: 0.01,
+			Action: func(p *mned.Point) bool {
+				if p.Value[0] < 0 {
+					backed = true
+				}
+				return true
+			},
+		},
+		mned.Event{
+			Cross: func(p *mned.Point) float64 {
+				return p.Value[1]
+			},
+			Tolerance: 0.01,
+			Action: func(p *mned.Point) bool {
+				return p.Value[0] <= 0.9 || !backed
+			},
+		},
+	)
+	solution.StepToEnd()
+	coords := solution.PointCoords()
+	fmt.Printf("Took %v steps.\n", len(coords[0]))
+	plot, err := glot.NewPlot(2, true, false)
+	if err != nil {
+		fmt.Printf("Couldn't create plot: %v\n", err)
+		os.Exit(1)
+	}
+	if err := plot.AddPointGroup("Orbit", "lines", [][]float64{
+		coords[1], coords[2],
+	}); err != nil {
+		fmt.Printf("Couldn't draw the points: %v\n", err)
+		os.Exit(1)
+	}
+}
diff --git a/examples/arenstorf.png b/examples/arenstorf.png
new file mode 100644
index 0000000..22cf99b
--- /dev/null
+++ b/examples/arenstorf.png
diff --git a/examples/arenstorf.txt b/examples/arenstorf.txt
new file mode 100644
index 0000000..b9d3337
--- /dev/null
+++ b/examples/arenstorf.txt
@@ -0,0 +1 @@
+Took 357 steps.
diff --git a/examples/exp.go b/examples/exp.go
new file mode 100644
index 0000000..1021be6
--- /dev/null
+++ b/examples/exp.go
@@ -0,0 +1,74 @@
+// Example of trapezium method for exponential solutions.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"math"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/method"
+	"github.com/JwanMan/mned/mstep"
+)
+
+type methodDescriptor struct {
+	name string
+	m    mned.Method
+}
+
+const step = 0.01
+const tol = 0.01
+
+var df mstep.YDerivative = func(t float64, v float64) float64 { return -1 }
+var methods = [5]methodDescriptor{
+	{name: "Euler", m: method.Euler(0.01)},
+	{name: "RK4", m: method.RK4(0.1)},
+	{name: "Implicit Euler", m: mstep.ImplicitEuler(0.01, 0.01, df)},
+	{name: "Trapezium", m: mstep.Trapezium(0.01, 0.01, df)},
+	{name: "RKF", m: method.RKFehlberg(0.01, 0.01, 1e-7, 0.1)},
+}
+
+func main() {
+	problem := mned.IVP{
+		Derivative: func(p mned.Point) ([]float64, bool) {
+			return []float64{-p.Value[0]}, true
+		},
+		Start: mned.Point{Time: 0, Value: []float64{1}},
+	}
+	for _, m := range methods {
+		solution, _ := mned.DenseSolve(
+			m.m, &problem, 0, 10, mned.HermiteForIVP(&problem),
+		)
+		points := solution.InnerPoints()
+		maxerr := float64(0)
+		for _, p := range points {
+			realVal := math.Exp(-p.Time)
+			err := math.Abs(p.Value[0] / realVal)
+			if err < 1 {
+				err = 1 / err
+			}
+			err -= 1
+			if err > maxerr {
+				maxerr = err
+			}
+		}
+		fmt.Printf("%v method got %v steps, max rel err = %v.\n",
+			m.name, len(points), maxerr)
+	}
+}
diff --git a/examples/exp.txt b/examples/exp.txt
new file mode 100644
index 0000000..3c3eddf
--- /dev/null
+++ b/examples/exp.txt
@@ -0,0 +1,5 @@
+Euler method got 1002 steps, max rel err = 0.05167716448732684.
+RK4 method got 102 steps, max rel err = 9.149053072698976e-06.
+Implicit Euler method got 1002 steps, max rel err = 0.05097553729673554.
+Trapezium method got 1002 steps, max rel err = 8.342139748207522e-05.
+RKF method got 103 steps, max rel err = 1.5683599816629368e-06.
diff --git a/examples/hooke.go b/examples/hooke.go
new file mode 100644
index 0000000..36229c1
--- /dev/null
+++ b/examples/hooke.go
@@ -0,0 +1,77 @@
+// Example visualization of Hooke's law.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"github.com/Arafatk/glot"
+	"mned"
+	"mned/ivp"
+	"mned/method"
+)
+
+func spring1(speed float64) {
+	spring := ivp.HookeSpring{
+		Mass:       1,
+		Spring:     1.5,
+		Length:     0.7,
+		Friction:   0.3,
+		Amplitude:  0.4,
+		AngleSpeed: speed,
+	}
+	problem := spring.ToIVP()
+	plot, err := glot.NewPlot(2, true, false)
+	if err != nil {
+		fmt.Printf("Creating the plot: %v\n", err)
+		return
+	}
+	if err := plot.SetXLabel(
+		fmt.Sprintf("Time (s) -- w = %v", speed),
+	); err != nil {
+		fmt.Printf("Drawing the X label: %v\n", err)
+		return
+	}
+
+	solution, _ := mned.DenseSolve(
+		method.RKFehlberg(0.0001, 0.01, 1e-7, 0.5),
+		&problem,
+		0,
+		40,
+		mned.HermiteForIVP(&problem),
+	)
+	coord := solution.PointCoords()
+	if err := plot.AddPointGroup("Position (m)", "lines", [][]float64{
+		coord[0], coord[1],
+	}); err != nil {
+		fmt.Printf("Drawing the position: %v\n", err)
+	}
+	if err := plot.AddPointGroup("Velocity (m/s)", "lines", [][]float64{
+		coord[0], coord[2],
+	}); err != nil {
+		fmt.Printf("Drawing the velocity: %v\n", err)
+	}
+}
+
+func main() {
+	spring1(0)
+	spring1(1.3)
+	spring1(2.4)
+	spring1(3.5)
+}
diff --git a/examples/hooke.png b/examples/hooke.png
new file mode 100644
index 0000000..f8bbd32
--- /dev/null
+++ b/examples/hooke.png
diff --git a/examples/hooke_check.go b/examples/hooke_check.go
new file mode 100644
index 0000000..c51339c
--- /dev/null
+++ b/examples/hooke_check.go
@@ -0,0 +1,91 @@
+// Example comparing the formula and numeric approximation for Hooke's law
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"math"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/ivp"
+	"github.com/JwanMan/mned/method"
+)
+
+type methodDescriptor struct {
+	name string
+	m    mned.Method
+}
+
+const step = 0.01
+const tol = 0.01
+
+var methods = [4]methodDescriptor{
+	{name: "Euler", m: method.Euler(step)},
+	{name: "RK4", m: method.RK4(step)},
+	{name: "Adaptive RK4", m: method.AdaptiveRK4(tol, step, 1e-6, 1)},
+	{name: "RK-Fehlberg", m: method.RKFehlberg(tol, step, 1e-6, 1)},
+}
+
+func normdiff(p1 []float64, p2 []float64) float64 {
+	var norm float64 = 0
+	for i, v := range p1 {
+		norm += (v - p2[i]) * (v - p2[i])
+	}
+	return norm
+}
+
+func maxerr(points []mned.Point, sol func(float64) []float64) float64 {
+	var max float64 = 0
+	for _, p := range points {
+		actual := sol(p.Time)
+		norm := normdiff(actual, p.Value)
+		if norm > max {
+			max = norm
+		}
+	}
+	return max
+}
+
+func main() {
+	spring := ivp.HookeSpring{
+		Mass:   1,
+		Spring: 0.7,
+		X0:     1,
+		Length: 0.7,
+	}
+	problem := spring.ToIVP()
+	amplitude := spring.X0 - spring.Length
+	sqratio := math.Sqrt(spring.Spring / spring.Mass)
+	realSol := func(t float64) []float64 {
+		return []float64{
+			amplitude*math.Cos(sqratio*t) + spring.Length,
+			-amplitude * sqratio * math.Sin(sqratio*t),
+		}
+	}
+	interp := mned.HermiteForIVP(&problem)
+
+	for _, m := range methods {
+		solution, _ := mned.DenseSolve(
+			m.m, &problem, 0, 40, interp,
+		)
+		points := solution.InnerPoints()
+		fmt.Printf("%v: Error of %v in %v steps.\n",
+			m.name, maxerr(points, realSol), len(points))
+	}
+}
diff --git a/examples/hooke_check.txt b/examples/hooke_check.txt
new file mode 100644
index 0000000..8081664
--- /dev/null
+++ b/examples/hooke_check.txt
@@ -0,0 +1,4 @@
+Euler: Error of 0.0017828128093108954 in 4001 steps.
+RK4: Error of 1.6399864728279775e-19 in 4001 steps.
+Adaptive RK4: Error of 7.496803458130451e-07 in 45 steps.
+RK-Fehlberg: Error of 7.835831771324868e-05 in 45 steps.
diff --git a/examples/parabolic.go b/examples/parabolic.go
new file mode 100644
index 0000000..f45872b
--- /dev/null
+++ b/examples/parabolic.go
@@ -0,0 +1,110 @@
+// Example visualization of parabolic throw with air friction.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"github.com/Arafatk/glot"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/ivp"
+	"github.com/JwanMan/mned/method"
+	"os"
+)
+
+func addPlot(
+	plot *glot.Plot, name string, arrays [][]float64, idx int, idy int,
+) {
+	if err := plot.AddPointGroup(
+		name,
+		"lines",
+		[][]float64{arrays[idx], arrays[idy]}); err != nil {
+		fmt.Printf("Error adding %v: %w\n", name, err)
+		os.Exit(1)
+	}
+}
+
+func main() {
+	var end mned.Point
+	var err error
+
+	throw := ivp.ParabolicThrow{
+		Height:     300,
+		V0:         [2]float64{100, 0},
+		Mass:       1,
+		Resistance: 0.01,
+		Gravity:    ivp.EARTH_GRAVITY,
+	}
+	problem := throw.ToIVP()
+	solution := mned.CacheSolve(
+		method.Euler(0.01),
+		&problem,
+		mned.LinearInterpolator{},
+		mned.Event{
+			Cross: func(p *mned.Point) float64 {
+				return p.Value[1]
+			},
+			Tolerance: 0.01,
+			Action: func(p *mned.Point) bool {
+				end = *p
+				return false
+			},
+		})
+	solution.StepToEnd()
+	fmt.Printf(
+		"Reached ground after %v s at %v m, with (%v, %v) m/s.\n",
+		end.Time,
+		end.Value[0],
+		end.Value[2],
+		end.Value[3],
+	)
+
+	arrays := solution.PointCoords()
+	byTime, err := glot.NewPlot(2, true, false)
+	if err != nil {
+		fmt.Printf("Creating plot window: %w\n")
+		os.Exit(1)
+	}
+	defer byTime.Close()
+	if err := byTime.SetXLabel("Time (s)"); err != nil {
+		fmt.Printf("Adding X label: %w\n")
+		os.Exit(1)
+	}
+	addPlot(byTime, "X position (m)", arrays, 0, 1)
+	addPlot(byTime, "Y position (m)", arrays, 0, 2)
+	addPlot(byTime, "X speed (m/s)", arrays, 0, 3)
+	addPlot(byTime, "Y speed (m/s)", arrays, 0, 4)
+
+	byTraj, err := glot.NewPlot(2, true, false)
+	if err != nil {
+		fmt.Printf("Creating plot window: %w\n")
+		os.Exit(1)
+	}
+	defer byTraj.Close()
+	if err := byTraj.SetXLabel("X coordinate"); err != nil {
+		fmt.Printf("Adding X label: %w\n")
+		os.Exit(1)
+	}
+	if err := byTraj.SetYLabel("Y coordinate"); err != nil {
+		fmt.Printf("Adding Y label: %w\n")
+		os.Exit(1)
+	}
+	addPlot(byTraj, "Trajectory (m)", arrays, 1, 2)
+	addPlot(byTraj, "Speed trajectory (m/s)", arrays, 3, 4)
+}
diff --git a/examples/parabolic.png b/examples/parabolic.png
new file mode 100644
index 0000000..a8dffb6
--- /dev/null
+++ b/examples/parabolic.png
diff --git a/examples/parabolic_check.go b/examples/parabolic_check.go
new file mode 100644
index 0000000..4189fe9
--- /dev/null
+++ b/examples/parabolic_check.go
@@ -0,0 +1,108 @@
+// Example comparing the formula and numeric approximation for parabolic throw.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package main
+
+import (
+	"fmt"
+	"math"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/ivp"
+	"github.com/JwanMan/mned/method"
+)
+
+type methodDescriptor struct {
+	name string
+	m    mned.Method
+}
+
+const step = 0.01
+
+var methods = [3]methodDescriptor{
+	{name: "Euler", m: method.Euler(step)},
+	{name: "Modified Euler", m: method.ModifiedEuler(step)},
+	{name: "RK4", m: method.RK4(step)},
+}
+
+func normdiff(p1 []float64, p2 []float64) float64 {
+	var norm float64 = 0
+	for i, v := range p1 {
+		norm += (v - p2[i]) * (v - p2[i])
+	}
+	return norm
+}
+
+func maxerr(points []mned.Point, sol func(float64) []float64) float64 {
+	var max float64 = 0
+	for _, p := range points {
+		actual := sol(p.Time)
+		norm := normdiff(actual, p.Value)
+		if norm > max {
+			max = norm
+		}
+	}
+	return max
+}
+
+func main() {
+	// Here we test the precision of different methods against the known
+	// problem of parabolic throw without air resistance. The problem is
+	// x''(t) = -gj, so x(t) = -g(t^2/2)j + v0*t + hj, thus
+	// x(t)[1] = 0 when t = (v0[1] + sqrt(v0[1]^2 + 2hg)) / g
+	var fallen float64
+	throw := ivp.ParabolicThrow{
+		Height:  300,
+		V0:      [2]float64{100, 0},
+		Mass:    1,
+		Gravity: ivp.EARTH_GRAVITY,
+	}
+	problem := throw.ToIVP()
+	realSolution := func(t float64) []float64 {
+		return []float64{
+			throw.V0[0] * t,
+			throw.Height + throw.V0[1]*t - throw.Gravity*t*t/2,
+			throw.V0[0],
+			throw.V0[1] - throw.Gravity*t,
+		}
+	}
+	realEnd := (throw.V0[1] + math.Sqrt(
+		throw.V0[1]*throw.V0[1]+2*throw.Height*throw.Gravity,
+	)) / throw.Gravity * throw.V0[0]
+
+	for _, method := range methods {
+		solution := mned.CacheSolve(
+			method.m, &problem, mned.LinearInterpolator{},
+			mned.Event{
+				Cross: func(p *mned.Point) float64 {
+					return p.Value[1]
+				},
+				Tolerance: 0.01,
+				Action: func(p *mned.Point) bool {
+					fallen = p.Value[0]
+					return false
+				},
+			})
+
+		solution.StepToEnd()
+		soldiff := maxerr(solution.ForwardPoints(), realSolution)
+		enddiff := math.Abs(fallen - realEnd)
+		fmt.Printf("%v: Point diff @ %v, end diff @ %v.\n",
+			method.name, soldiff, enddiff)
+	}
+}
diff --git a/examples/parabolic_check.txt b/examples/parabolic_check.txt
new file mode 100644
index 0000000..373b723
--- /dev/null
+++ b/examples/parabolic_check.txt
@@ -0,0 +1,3 @@
+Euler: Point diff @ 0.1470067554878007, end diff @ 0.49737244764885347.
+Modified Euler: Point diff @ 1.9920977079185658e-22, end diff @ 0.0026275523512602206.
+RK4: Point diff @ 1.9926361270567946e-22, end diff @ 0.0026275523512602206.
diff --git a/interpolate.go b/interpolate.go
new file mode 100644
index 0000000..075a534
--- /dev/null
+++ b/interpolate.go
@@ -0,0 +1,97 @@
+// Interpolation routines.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+// An Interpolator tries to approximate a function based on a finite set of
+// known points in it. We'll focus just on interpolation with two points since
+// this is precise enough for our purposes, although the interpolator might be
+// constructed with some knowledge about the underlying function.
+type Interpolator interface {
+	// Approximate the value of a function in a point `t` given two
+	// **different** points `p1` and `p2` of the function such that `t`
+	// is between `p1.Time` and `p2.Time`.
+	FindValue(p1 *Point, p2 *Point, t float64) []float64
+}
+
+// A LinearInterpolator is an interpolator that assumes that there's a straight
+// line between the two points given. It can be constructed with `new`.
+type LinearInterpolator struct{}
+
+func (li LinearInterpolator) FindValue(
+	p1 *Point, p2 *Point, t float64,
+) []float64 {
+	// The line between (t1,x1) and (t2,x2) is x1 + (x2-x1)*(t-t1)/(t2-t1).
+	// If ratio:=(t-t1)/(t2-t1), this is x2*ratio + x1*(1-ratio)
+	size := len(p1.Value)
+	result := make([]float64, size)
+	ratio := (t - p1.Time) / (p2.Time - p1.Time)
+	for i := 0; i < size; i++ {
+		result[i] = p1.Value[i]*(1-ratio) + p2.Value[i]*ratio
+	}
+	return result
+}
+
+// A HermiteInterpolator is an interpolator that takes into account the
+// derivative of the solution function in the end points.
+type HermiteInterpolator struct {
+	F func(Point) ([]float64, bool) // The derivative function.
+}
+
+func (h HermiteInterpolator) FindValue(
+	p1 *Point, p2 *Point, t float64,
+) []float64 {
+	diff := p2.Time - p1.Time
+	off := t - p1.Time
+	d01, ok := h.F(*p1)
+	if !ok { // Can't derive: fallback
+		return LinearInterpolator{}.FindValue(p1, p2, t)
+	}
+	d23, ok := h.F(*p2)
+	if !ok {
+		return LinearInterpolator{}.FindValue(p1, p2, t)
+	}
+	d12 := make([]float64, len(d01))
+	for i := range d12 {
+		d12[i] = (p2.Value[i] - p1.Value[i]) / diff
+	}
+	d02 := make([]float64, len(d01))
+	for i := range d02 {
+		d02[i] = (d12[i] - d01[i]) / diff
+	}
+	d13 := make([]float64, len(d01))
+	for i := range d13 {
+		d13[i] = (d23[i] - d12[i]) / diff
+	}
+	d03 := make([]float64, len(d01))
+	for i := range d03 {
+		d03[i] = (d13[i] - d02[i]) / diff
+	}
+	result := make([]float64, len(d01))
+	for i := range result {
+		result[i] = p1.Value[i] + off*
+			(d01[i]+off*(d02[i]+(off-diff)*d03[i]))
+	}
+	return result
+}
+
+// Create a Hermite interpolator suitable for the given IVP.
+func HermiteForIVP(ivp *IVP) HermiteInterpolator {
+	return HermiteInterpolator{F: ivp.Derivative}
+}
diff --git a/ivp/arenstorf.go b/ivp/arenstorf.go
new file mode 100644
index 0000000..05b1527
--- /dev/null
+++ b/ivp/arenstorf.go
@@ -0,0 +1,68 @@
+// Arenstorf orbit problem.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package ivp
+
+import (
+	"math"
+	"github.com/JwanMan/mned"
+)
+
+var arenstorfOrbit = mned.IVP{
+	Derivative: func(p mned.Point) ([]float64, bool) {
+		const μ = 0.012277471
+		const μ2 = 1 - μ
+		xμ := p.Value[0] + μ
+		xμ2 := p.Value[0] - μ2
+		D1 := xμ*xμ + p.Value[1]*p.Value[1]
+		D1 *= math.Sqrt(D1)
+		D2 := xμ2*xμ2 + p.Value[1]*p.Value[1]
+		D2 *= math.Sqrt(D2)
+		return []float64{
+			p.Value[2],
+			p.Value[3],
+			p.Value[0] + 2*p.Value[3] - μ2*xμ/D1 - μ*xμ2/D2,
+			p.Value[1]*(1-μ2/D1-μ/D2) - 2*p.Value[2],
+		}, true
+	},
+	Start: mned.Point{
+		Time:  0,
+		Value: []float64{0.994, 0, 0, -2.001585106},
+	},
+}
+
+// The Arenstorf orbit is a stable orbit of a satellite around the Earth taking
+// into account the gravity of the Earth and the Moon and considering the
+// satellite to have a negligible mass compared to both. The orbit was
+// discovered by Richard Arenstorf and is given by the following formulae:
+//
+// ```
+// μ = 0.012277471, μ' = 1 - μ;
+// D1 = ((x+μ)^2 + y^2)^(3/2), D2 = ((x-μ')^2 + y^2)^(3/2);
+// x" = x + 2y' - μ'(x+μ)/D1 - μ(x-μ')/D2;
+// y" = y - 2x' - μ'y/D1 - μy/D2;
+// x(0) = 0.994, y(0) = 0, x'(0) = 0, y'(0) = -2.001585106.
+// ```
+//
+// The orbit has three big loops at one side and a small one at the other, and
+// varying μ changes the number of loops. Because of the unstability, this
+// problem is often used to test different ODE solving methods.
+func ArenstorfOrbit() mned.IVP {
+	return arenstorfOrbit.Clone()
+}
diff --git a/ivp/spring.go b/ivp/spring.go
new file mode 100644
index 0000000..4c8798d
--- /dev/null
+++ b/ivp/spring.go
@@ -0,0 +1,79 @@
+// Spring movement problem based on Hooke's law.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package ivp
+
+import (
+	"math"
+	"github.com/JwanMan/mned"
+)
+
+// A HookeSpring is a spring moving following Hooke's law.
+//
+// The spring is laid out horizontally over a table and has one end attached to
+// a body which receives a certain friction from the table and such that the
+// mass of the spring is negligible compared to the mass of the body. The other
+// end of the spring is attached to a harmonic oscillator. The position of the
+// body is measured as the (signed) distance to the rest position of the
+// oscillator, and the velocity is measured in terms of such position.
+//
+// Only Mass and Spring are required, although the result of only specifying
+// that would be really boring. Note also that, since we're following Hooke's
+// law, the length of the spring can get negative.
+type HookeSpring struct {
+	Mass       float64 // Mass of the body (kg).
+	Spring     float64 // Hooke's constant of the spring (kg/s^2).
+	X0         float64 // Initial position of the body (m).
+	V0         float64 // Initial velocity of the body (m/s).
+	Length     float64 // Length of the spring at rest (m).
+	Friction   float64 // Friction between the body and the table (kg/s).
+	Amplitude  float64 // Amplitude of the oscillator's force (N).
+	AngleSpeed float64 // Angle speed of the oscillator (rad/s).
+}
+
+// Get the IVP associated to the problem. The solution values have the form
+// (length, speed), time starts at 0, and the harmonic oscillator starts at
+// its natural position.
+//
+// The ODE is `m * x''(t) = A*sin(wt) - k*(x(t)-l) - rx'(t)`, where m is the
+// mass of the body, x is its position, A is the amplitude of the oscillator,
+// w is its angular velocity, k is the spring's constant, l is its length at
+// rest, and r is the friction coefficient.
+func (h *HookeSpring) ToIVP() mned.IVP {
+	adjSpring := h.Spring / h.Mass
+	adjFriction := h.Friction / h.Mass
+	adjAmplitude := h.Amplitude / h.Mass
+	w := h.AngleSpeed
+	length := h.Length
+
+	return mned.IVP{
+		Derivative: func(p mned.Point) ([]float64, bool) {
+			return []float64{
+				p.Value[1],
+				adjAmplitude*math.Sin(w*p.Time) -
+					adjSpring*(p.Value[0]-length) -
+					adjFriction*p.Value[1],
+			}, true
+		},
+		Start: mned.Point{
+			Time:  0,
+			Value: []float64{h.X0, h.V0},
+		},
+	}
+}
diff --git a/ivp/throw.go b/ivp/throw.go
new file mode 100644
index 0000000..e91d989
--- /dev/null
+++ b/ivp/throw.go
@@ -0,0 +1,68 @@
+// Parabolic throw problem with air friction.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package ivp
+
+import (
+	"math"
+	"github.com/JwanMan/mned"
+)
+
+// Gravity acceleration in the Earth surface
+const EARTH_GRAVITY float64 = 9.806
+
+// A ParabolicThrow is the result of throwing an object, which we consider to
+// be punctual, from a given height and with some initial velocity. We assume
+// that Newton mechanics apply and that air resistance has a force proportional
+// to the square of the velocity, and opposite to it.
+//
+// The resulting ODE is `x'' = -gj - (k/m)*x'*|x'|`, where `x` is the position,
+// `g` is the gravity acceleration, `j` is the vertical unit vector (upwards),
+// `k` is the air resistance, and `m` is the mass of the object.
+//
+// The comments on fields assume standard SI units, but other units may be used.
+type ParabolicThrow struct {
+	Height     float64    // Initial height over the floor (m).
+	Mass       float64    // Mass of the object (kg). *Must* be positive.
+	V0         [2]float64 // Initial velocity (x,y) (m/s).
+	Resistance float64    // Air resistance (kg/m).
+	Gravity    float64    // Gravity acceleration (m/s^2).
+}
+
+// Create an IVP from the data on this problem.
+//
+// Result values have the form (x,y,vx,vy), where (x,y) is the position and
+// (vx,vy) is the velocity. The initial x and time are 0.
+func (p *ParabolicThrow) ToIVP() mned.IVP {
+	ratio := p.Resistance / p.Mass
+	gravity := p.Gravity
+	return mned.IVP{
+		Derivative: func(p mned.Point) ([]float64, bool) {
+			x := p.Value
+			air := -ratio * math.Sqrt(x[2]*x[2]+x[3]*x[3])
+			return []float64{
+				x[2], x[3], air * x[2], air*x[3] - gravity,
+			}, true
+		},
+		Start: mned.Point{
+			Time:  0,
+			Value: []float64{0, p.Height, p.V0[0], p.V0[1]},
+		},
+	}
+}
diff --git a/method/euler.go b/method/euler.go
new file mode 100644
index 0000000..8370ae9
--- /dev/null
+++ b/method/euler.go
@@ -0,0 +1,79 @@
+// Euler method and variants.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package method
+
+import "github.com/JwanMan/mned"
+
+func eulerStep(ivp *mned.IVP, h float64) bool {
+	deriv, ok := ivp.Derivative(ivp.Start)
+	if !ok {
+		return false
+	}
+	for i := 0; i < len(deriv); i++ {
+		ivp.Start.Value[i] += h * deriv[i]
+	}
+	return true
+}
+
+// Make an instance of the Euler method with the given step, which must be
+// positive. For a problem `x'(t)=f(t,x(t)), x(t0)=x0`, a step of the Euler
+// method is given by `x(t0+step) = x(t0) + step*f(t0,x(t0))`.
+func Euler(step float64) mned.Method {
+	return mned.FixedStepMethod{Step: step, Method: eulerStep}
+}
+
+// Make an instance of the adaptive version of the Euler method given the
+// tolerance, the initial step, and the minimum and maximum steps.
+func AdaptiveEuler(
+	tolerance float64, step float64, hmin float64, hmax float64,
+) mned.Method {
+	return &mned.AdaptiveStepMethod{
+		Step:      step,
+		Hmin:      hmin,
+		Hmax:      hmax,
+		Tolerance: tolerance,
+		Method:    eulerStep,
+		Order:     1,
+	}
+}
+
+func modEulerStep(ivp *mned.IVP, h float64) bool {
+	deriv, ok := ivp.Derivative(ivp.Start)
+	if !ok {
+		return false
+	}
+	point := ivp.Start.Clone()
+	point.Time += h
+	for i, v := range deriv {
+		point.Value[i] += h * v
+	}
+	deriv2, ok := ivp.Derivative(point)
+	for i, _ := range deriv {
+		ivp.Start.Value[i] += h * (deriv[i] + deriv2[i]) / 2
+	}
+	return true
+}
+
+// Make an instance of the modified Euler method with the given step, which must
+// be positive. This is a FixedStepMethod given by
+// `x(t+h)=x(t) + h/2 * (f(t,x(t)) + f(t+h,x(t)+h*f(t,x(t))))`.
+func ModifiedEuler(step float64) mned.Method {
+	return mned.FixedStepMethod{Step: step, Method: modEulerStep}
+}
diff --git a/method/fehlberg.go b/method/fehlberg.go
new file mode 100644
index 0000000..3ea14d5
--- /dev/null
+++ b/method/fehlberg.go
@@ -0,0 +1,208 @@
+// Runge-Kutta-Fehlberg method.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package method
+
+import (
+	"math"
+	"github.com/JwanMan/mned"
+)
+
+type rkFehlbergStepper struct {
+	current mned.Point
+	deriv   func(mned.Point) ([]float64, bool)
+	step    float64
+	tol     float64
+	hmin    float64
+	hmax    float64
+}
+
+func (s *rkFehlbergStepper) Next() (*mned.Point, bool) {
+	for s.step >= s.hmin {
+		k1, ok := s.deriv(s.current)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k1 {
+			k1[i] *= s.step
+		}
+
+		p2 := s.current.Clone()
+		p2.Time += s.step / 4
+		for i := range p2.Value {
+			p2.Value[i] += k1[i] / 4
+		}
+		k2, ok := s.deriv(p2)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k2 {
+			k2[i] *= s.step
+		}
+
+		p3 := s.current.Clone()
+		p3.Time += s.step * 3 / 8
+		for i := range p3.Value {
+			p3.Value[i] += (3*k1[i] + 9*k2[i]) / 32
+		}
+		k3, ok := s.deriv(p3)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k3 {
+			k3[i] *= s.step
+		}
+
+		p4 := s.current.Clone()
+		p4.Time += s.step * 12 / 13
+		for i := range p4.Value {
+			p4.Value[i] +=
+				(1932*k1[i] - 7200*k2[i] + 7296*k3[i]) / 2197
+		}
+		k4, ok := s.deriv(p4)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k4 {
+			k4[i] *= s.step
+		}
+
+		p5 := s.current.Clone()
+		p5.Time += s.step
+		for i := range p5.Value {
+			p5.Value[i] += 439*k1[i]/216 - 8*k2[i] +
+				3680*k3[i]/513 - 845*k4[i]/4104
+		}
+		k5, ok := s.deriv(p5)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k5 {
+			k5[i] *= s.step
+		}
+
+		p6 := s.current.Clone()
+		p6.Time += s.step / 2
+		for i := range p6.Value {
+			p6.Value[i] +=
+				-8*k1[i]/27 + 2*k2[i] - 3544*k3[i]/2656 +
+					1859*k4[i]/4104 - 11*k5[i]/40
+		}
+		k6, ok := s.deriv(p6)
+		if !ok {
+			s.step *= 0.5
+			continue
+		}
+		for i := range k6 {
+			k6[i] *= s.step
+		}
+
+		sqerr := float64(0)
+		// Get \tilde{w}_{i+1}-w_{i+1} directly:
+		// * (- 16/135 25/216)
+		// 1/360
+		// * (- 6656/12825 1408/2565)
+		// -128/4275
+		// * (- 28561/56430 2197/4104)
+		// -2197/75240
+		// * (- -9/50 -1/5)
+		// 1/50
+		for i := range s.current.Value {
+			v := k1[i]/360 - 128*k3[i]/4275 -
+				2197*k4[i]/75240 + k5[i]/50 + 2*k6[i]/55
+			sqerr += v * v
+		}
+		err := math.Sqrt(sqerr)
+
+		if err > s.tol {
+			s.step = mned.AdaptiveUpdateStep(
+				s.step, s.tol, err, s.hmax, 4,
+			)
+			continue
+		}
+
+		for i := range s.current.Value {
+			s.current.Value[i] += 25*k1[i]/216 + 1408*k3[i]/2565 +
+				2197*k4[i]/4104 - k5[i]/5
+		}
+		s.current.Time += s.step
+		s.step = mned.AdaptiveUpdateStep(s.step, s.tol, err, s.hmax, 4)
+		return &s.current, true
+	}
+	return nil, false
+}
+
+type rkFehlbergMethod struct {
+	step float64
+	tol  float64
+	hmin float64
+	hmax float64
+}
+
+func (m *rkFehlbergMethod) withStep(step float64, ivp *mned.IVP) mned.Stepper {
+	return &rkFehlbergStepper{
+		step:    step,
+		tol:     m.tol,
+		hmin:    m.hmin,
+		hmax:    m.hmax,
+		current: ivp.Start.Clone(),
+		deriv:   ivp.Derivative,
+	}
+}
+
+func (m *rkFehlbergMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(m.step, ivp)
+}
+
+func (m *rkFehlbergMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-m.step, ivp)
+}
+
+// Create an instance of the Runge-Kutta-Fehlberg method, an adaptive,
+// fourth-order method, with the given parameter.
+//
+// This method is given by the following parameters:
+// ```
+// k1 = hf(t, x(t))
+// k2 = hf(t+h/4, x(t)+k1/4)
+// k3 = hf(t+3h/8, x(t)+3k1/32+9k2/32)
+// k4 = hf(t+12h/13, x(t)+1932k1/2197-7200k2/2197+7296k3/2197)
+// k5 = hf(t+h, x(t)+439k1/216-8k2+3680k3/513-845k4/4104)
+// k6 = hf(t+h/2, x(t)-8k1/27+2k2-3544k3/2656+1859k4/4104-11k5/40)
+// e = k1/360-128k3/4275-2197k4/75240+k5/50+2k6/55
+// x(t+h) ~ x(t)+25k1/216+1408k3/2565+2197k4/4104-k5/5
+// ```
+// If `|e|` is greater than `tol`, we let `q:=(h*tol/(2*|e|))^(1/4)`, saturated
+// into `[0.1,4]`, and redo the operations with `h*q`, saturating over `hmax`
+// and stopping below `hmin`.
+func RKFehlberg(
+	tolerance float64, step float64, hmin float64, hmax float64,
+) mned.Method {
+	return &rkFehlbergMethod{
+		step: step,
+		tol:  tolerance,
+		hmin: hmin,
+		hmax: hmax,
+	}
+}
diff --git a/method/rk.go b/method/rk.go
new file mode 100644
index 0000000..391fcb8
--- /dev/null
+++ b/method/rk.go
@@ -0,0 +1,107 @@
+// Runge-Kutta method and adaptive variant.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package method
+
+import "github.com/JwanMan/mned"
+
+func RK4Step(ivp *mned.IVP, h float64) ([]float64, bool) {
+	k1, ok := ivp.Derivative(ivp.Start)
+	if !ok {
+		return nil, false
+	}
+	for i, _ := range k1 {
+		k1[i] *= h
+	}
+
+	p1 := ivp.Start.Clone()
+	p1.Time += h / 2
+	for i, v := range k1 {
+		p1.Value[i] += v / 2
+	}
+	k2, ok := ivp.Derivative(p1)
+	if !ok {
+		return nil, false
+	}
+	for i, _ := range k2 {
+		k2[i] *= h
+	}
+
+	p2 := ivp.Start.Clone()
+	p2.Time += h / 2
+	for i, v := range k2 {
+		p2.Value[i] += v / 2
+	}
+	k3, ok := ivp.Derivative(p2)
+	if !ok {
+		return nil, false
+	}
+	for i, _ := range k3 {
+		k3[i] *= h
+	}
+
+	p3 := ivp.Start.Clone()
+	p3.Time += h
+	for i, v := range k3 {
+		p3.Value[i] += v
+	}
+	k4, ok := ivp.Derivative(p3)
+	if !ok {
+		return nil, false
+	}
+	for i, _ := range k4 {
+		k4[i] *= h
+	}
+
+	for i, _ := range ivp.Start.Value {
+		ivp.Start.Value[i] += (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) / 6
+	}
+	return k1, true
+}
+
+func rk4Step(ivp *mned.IVP, h float64) bool {
+	_, ok := RK4Step(ivp, h)
+	return ok
+}
+
+// Initialize the 4th-order Runge-Kutta method, a.k.a. RK4. This is a fixed step
+// method given by the following formulae:
+// ```
+// k1 = h*f(t, x(t))
+// k2 = h*f(t+h/2, x(t)+k1/2)
+// k3 = h*f(t+h/2, x(t)+k2/2)
+// k4 = h*f(t+h, x(t)+k3)
+// x(t+h) = x(t) + (k1 + 2*k2 + 2*k3 + k4) / 6
+// ```
+func RK4(step float64) mned.Method {
+	return mned.FixedStepMethod{Step: step, Method: rk4Step}
+}
+
+func AdaptiveRK4(
+	tolerance float64, step float64, hmin float64, hmax float64,
+) mned.Method {
+	return &mned.AdaptiveStepMethod{
+		Step:      step,
+		Hmin:      hmin,
+		Hmax:      hmax,
+		Tolerance: tolerance,
+		Method:    rk4Step,
+		Order:     4,
+	}
+}
diff --git a/mstep/adams.go b/mstep/adams.go
new file mode 100644
index 0000000..e74269c
--- /dev/null
+++ b/mstep/adams.go
@@ -0,0 +1,105 @@
+// Adams stepping methods.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mstep
+
+import (
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/method"
+)
+
+func adamsBashfordStep(steps []Step, h float64) []float64 {
+	return Implicit(steps, []float64{1}, []float64{
+		float64(55) / 24, float64(-59) / 24,
+		float64(37) / 24, float64(-9) / 24,
+	}, h)
+}
+
+func adamsMoultonStep(steps []Step, h float64) []float64 {
+	return Implicit(steps, []float64{0, 1}, []float64{
+		float64(9) / 24, float64(19) / 24,
+		float64(-5) / 24, float64(1) / 24,
+	}, h)
+}
+
+type adamsBashfordStepper struct {
+	steps   [4]Step
+	h       float64
+	problem mned.IVP
+	init    uint8
+}
+
+func (s *adamsBashfordStepper) NextStep(h float64) (*mned.Point, bool) {
+	var ok bool
+	if s.h == 0 { // Error indicator
+		return nil, false
+	}
+	if s.init != 0 {
+		deriv, ok := method.RK4Step(&s.problem, h)
+		if !ok {
+			return nil, false
+		}
+		s.problem.Start.Time += h
+		s.steps[s.init].Deriv = deriv
+		s.init--
+		s.steps[s.init].Value = make([]float64, len(deriv))
+		copy(s.steps[s.init].Value, s.problem.Start.Value)
+		return &s.problem.Start, true
+	}
+
+	s.steps[0].Deriv, ok = s.problem.Derivative(s.problem.Start)
+	if !ok {
+		s.h = 0
+		return nil, false
+	}
+	s.problem.Start.Value = adamsBashfordStep(s.steps[:], h)
+	s.problem.Start.Time += h
+	Shift(s.steps[:])
+	s.steps[0].Value = s.problem.Start.Value
+	return &s.problem.Start, true
+}
+
+func (s *adamsBashfordStepper) Next() (*mned.Point, bool) {
+	return s.NextStep(s.h)
+}
+
+type adamsBashfordMethod float64
+
+func (m adamsBashfordMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper {
+	result := adamsBashfordStepper{
+		h:       h,
+		problem: ivp.Clone(),
+		init:    3,
+	}
+	result.steps[3].Value = make([]float64, len(ivp.Start.Value))
+	copy(result.steps[3].Value, ivp.Start.Value)
+	return &result
+}
+
+func (m adamsBashfordMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(float64(m), ivp)
+}
+
+func (m adamsBashfordMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-float64(m), ivp)
+}
+
+func AdamsBashford(step float64) mned.Method {
+	return adamsBashfordMethod(step)
+}
diff --git a/mstep/common.go b/mstep/common.go
new file mode 100644
index 0000000..2ecff09
--- /dev/null
+++ b/mstep/common.go
@@ -0,0 +1,53 @@
+// Common code to define stepping methods.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mstep
+
+type Step struct {
+	Value []float64
+	Deriv []float64
+}
+
+func Implicit(steps []Step, a []float64, b []float64, step float64) []float64 {
+	if len(steps) == 0 {
+		panic("There must be at least one step!\n")
+	}
+	if len(steps) < len(a) || len(steps) < len(b) {
+		panic("There must be at least as many steps as coefficients.\n")
+	}
+
+	result := make([]float64, len(steps[0].Value))
+	for i, v := range a {
+		for j := range result {
+			result[j] += v * steps[i].Value[j]
+		}
+	}
+	for i, v := range b {
+		for j := range result {
+			result[j] += step * v * steps[i].Deriv[j]
+		}
+	}
+	return result
+}
+
+func Shift(steps []Step) {
+	for i := len(steps) - 1; i > 0; i-- {
+		steps[i] = steps[i-1]
+	}
+}
diff --git a/mstep/implicit.go b/mstep/implicit.go
new file mode 100644
index 0000000..228a3b8
--- /dev/null
+++ b/mstep/implicit.go
@@ -0,0 +1,164 @@
+// Implicit and trapezium stepping methods.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mstep
+
+import "github.com/JwanMan/mned"
+
+type YDerivative func(float64, float64) float64
+
+type implicitEulerStepper struct {
+	f         func(mned.Point) ([]float64, bool)
+	dfy       YDerivative
+	step      float64
+	tolerance float64
+	current   mned.Point
+}
+
+func (s *implicitEulerStepper) NextStep(h float64) (*mned.Point, bool) {
+	inner := s.f
+	time := s.current.Time + s.step
+	value := s.current.Value[0]
+	deriv := s.dfy
+	fn := func(w float64) (float64, bool) {
+		v, oki := inner(mned.Point{Time: time, Value: []float64{w}})
+		if oki {
+			return w - value - h*v[0], true
+		} else {
+			return 0, false
+		}
+	}
+	result, ok := NewtonRoot1D(
+		s.current.Value[0], s.tolerance, fn,
+		func(w float64) (float64, bool) {
+			return 1 - h*deriv(time, w), true
+		})
+	if !ok {
+		return nil, false
+	}
+	s.current.Value[0] = result
+	s.current.Time += s.step
+	return &s.current, true
+}
+
+func (s *implicitEulerStepper) Next() (*mned.Point, bool) {
+	return s.NextStep(s.step)
+}
+
+type implicitEulerMethod struct {
+	dfy       YDerivative
+	h         float64
+	tolerance float64
+}
+
+func (m *implicitEulerMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper {
+	return &implicitEulerStepper{
+		f:         ivp.Derivative,
+		current:   ivp.Start.Clone(),
+		step:      h,
+		dfy:       m.dfy,
+		tolerance: m.tolerance,
+	}
+}
+
+func (m *implicitEulerMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(m.h, ivp)
+}
+
+func (m *implicitEulerMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-m.h, ivp)
+}
+
+func ImplicitEuler(step, tolerance float64, df YDerivative) mned.Method {
+	return &implicitEulerMethod{h: step, tolerance: tolerance, dfy: df}
+}
+
+type trapeziumStepper struct {
+	f         func(mned.Point) ([]float64, bool)
+	dfy       YDerivative
+	step      float64
+	tolerance float64
+	current   mned.Point
+	prevf     float64
+}
+
+func (s *trapeziumStepper) NextStep(h float64) (*mned.Point, bool) {
+	inner := s.f
+	time := s.current.Time + s.step
+	pv, ok := s.f(mned.Point{
+		Time:  s.current.Time,
+		Value: []float64{s.current.Value[0]},
+	})
+	if !ok {
+		return nil, false
+	}
+	value := s.current.Value[0] + h/2*pv[0]
+	fn := func(w float64) (float64, bool) {
+		v, oki := inner(mned.Point{Time: time, Value: []float64{w}})
+		if oki {
+			return w - value - h/2*v[0], true
+		} else {
+			return 0, false
+		}
+	}
+	deriv := s.dfy
+	result, ok := NewtonRoot1D(
+		s.current.Value[0], s.tolerance, fn,
+		func(w float64) (float64, bool) {
+			return 1 - h/2*deriv(time, w), true
+		})
+	if !ok {
+		return nil, false
+	}
+	s.current.Value[0] = result
+	s.current.Time += s.step
+	return &s.current, true
+}
+
+func (s *trapeziumStepper) Next() (*mned.Point, bool) {
+	return s.NextStep(s.step)
+}
+
+type trapeziumMethod struct {
+	dfy       YDerivative
+	h         float64
+	tolerance float64
+}
+
+func (m *trapeziumMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper {
+	return &trapeziumStepper{
+		f:         ivp.Derivative,
+		current:   ivp.Start.Clone(),
+		step:      h,
+		dfy:       m.dfy,
+		tolerance: m.tolerance,
+	}
+}
+
+func (m *trapeziumMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(m.h, ivp)
+}
+
+func (m *trapeziumMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-m.h, ivp)
+}
+
+func Trapezium(step, tolerance float64, df YDerivative) mned.Method {
+	return &trapeziumMethod{h: step, tolerance: tolerance, dfy: df}
+}
diff --git a/mstep/pc.go b/mstep/pc.go
new file mode 100644
index 0000000..b6147bb
--- /dev/null
+++ b/mstep/pc.go
@@ -0,0 +1,300 @@
+// Correction routines for optimal step estimation.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mstep
+
+import (
+	"math"
+	"github.com/JwanMan/mned"
+	"github.com/JwanMan/mned/method"
+)
+
+type pcStepper struct {
+	steps   [5]Step
+	problem mned.IVP
+	h       float64
+	init    uint8
+}
+
+func (s *pcStepper) NextStep(h float64) (*mned.Point, bool) {
+	var ok bool
+	if s.h == 0 {
+		return nil, false
+	}
+	if s.init != 0 {
+		deriv, ok := method.RK4Step(&s.problem, h)
+		if !ok {
+			return nil, false
+		}
+		s.problem.Start.Time += h
+		s.steps[s.init+1].Deriv = deriv
+		s.steps[s.init].Value = make([]float64, len(deriv))
+		copy(s.steps[s.init].Value, s.problem.Start.Value)
+		return &s.problem.Start, true
+	}
+
+	s.steps[1].Deriv, ok = s.problem.Derivative(s.problem.Start)
+	if !ok {
+		s.h = 0
+		return nil, false
+	}
+	s.steps[0].Value = adamsBashfordStep(s.steps[1:], h)
+	s.problem.Start.Time += h
+	s.problem.Start.Value = s.steps[0].Value
+	s.steps[0].Deriv, ok = s.problem.Derivative(s.problem.Start)
+	if !ok {
+		s.problem.Start.Time -= h
+		s.problem.Start.Value = s.steps[1].Value
+		return nil, false
+	}
+	s.steps[0].Value = adamsMoultonStep(s.steps[:], h)
+	s.problem.Start.Value = s.steps[0].Value
+	Shift(s.steps[:])
+	s.steps[0].Value = nil
+	s.steps[0].Deriv = nil
+	s.steps[1].Deriv = nil
+	return &s.problem.Start, true
+}
+
+func (s *pcStepper) Next() (*mned.Point, bool) {
+	return s.NextStep(s.h)
+}
+
+type pcMethod float64
+
+func (m pcMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper {
+	result := pcStepper{
+		problem: ivp.Clone(),
+		h:       h,
+		init:    3,
+	}
+	result.steps[4].Value = make([]float64, len(ivp.Start.Value))
+	copy(result.steps[4].Value, ivp.Start.Value)
+	return &result
+}
+
+func (m pcMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(float64(m), ivp)
+}
+
+func (m pcMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-float64(m), ivp)
+}
+
+func PredictorCorrector(step float64) mned.Method {
+	return pcMethod(step)
+}
+
+type adaptivePCStepper struct {
+	steps     [5]Step
+	lasttime  float64
+	deriv     func(mned.Point) ([]float64, bool)
+	h         float64
+	hmin      float64
+	hmax      float64
+	tolerance float64
+	pos       uint8
+}
+
+func (s *adaptivePCStepper) tryStep() (errNorm float64, ok bool) {
+	s.steps[1].Deriv, ok = s.deriv(mned.Point{
+		Time:  s.lasttime,
+		Value: s.steps[1].Value,
+	})
+	if !ok {
+		return 0, false
+	}
+	s.steps[0].Value = adamsBashfordStep(s.steps[1:], s.h)
+
+	s.steps[0].Deriv, ok = s.deriv(mned.Point{
+		Time:  s.lasttime + s.h,
+		Value: s.steps[0].Value,
+	})
+	if !ok {
+		return 0, false
+	}
+	bash := s.steps[0].Value
+	moulton := adamsMoultonStep(s.steps[:], s.h)
+	s.steps[0].Value = moulton
+	s.steps[0].Deriv = nil
+
+	sqerr := float64(0)
+	for i, v := range s.steps[0].Value {
+		diff := v - bash[i]
+		sqerr += diff * diff
+	}
+	err := (math.Sqrt(sqerr) * 19) / 270
+	return err, true
+}
+
+func (s *adaptivePCStepper) adjustStep(err float64) {
+	q := math.Pow((s.tolerance*math.Abs(s.h))/(2*err), 0.25)
+	if q < 0.1 {
+		q = 0.1
+	}
+	if q > 4 {
+		q = 4
+	}
+	s.h *= q
+	if math.Abs(s.h) > s.hmax {
+		if s.h > 0 {
+			s.h = s.hmax
+		} else {
+			s.h = -s.hmax
+		}
+	}
+}
+
+func (s *adaptivePCStepper) init() bool {
+	var ok bool
+	size := len(s.steps[4].Value)
+	ivp := mned.IVP{Derivative: s.deriv}
+	for math.Abs(s.h) >= s.hmin {
+		ivp.Start.Time = s.lasttime
+		ivp.Start.Value = make([]float64, size)
+		copy(ivp.Start.Value, s.steps[4].Value)
+
+		for i := 4; i > 1; i-- {
+			s.steps[i].Deriv, ok = method.RK4Step(&ivp, s.h)
+			if !ok {
+				return false
+			}
+			ivp.Start.Time += s.h
+			s.steps[i-1].Value = make([]float64, size)
+			copy(s.steps[i-1].Value, ivp.Start.Value)
+		}
+		s.lasttime = ivp.Start.Time
+		errn, ok := s.tryStep()
+		s.lasttime += s.h
+		if !ok {
+			s.lasttime -= 4 * s.h
+			return false
+		}
+		tol := s.tolerance * s.h
+		if errn > tol {
+			s.lasttime -= 4 * s.h
+			s.adjustStep(errn)
+			continue
+		}
+		s.pos = 4
+		return true
+	}
+	return false
+}
+
+func (s *adaptivePCStepper) tryInit() bool {
+	for math.Abs(s.h) >= s.hmin {
+		if s.init() {
+			Shift(s.steps[:])
+			return true
+		}
+		for i := 0; i < 4; i++ {
+			s.steps[i].Value = nil
+			s.steps[i].Deriv = nil
+		}
+		s.steps[4].Deriv = nil
+		s.h /= 2
+	}
+	return false
+}
+
+func (s *adaptivePCStepper) Next() (*mned.Point, bool) {
+	var result mned.Point
+
+	if s.pos == 5 {
+		if !s.tryInit() {
+			return nil, false
+		}
+	}
+	for s.h >= s.hmin {
+		if s.pos > 0 {
+			result.Time = s.lasttime - s.h*float64(s.pos-1)
+			result.Value = s.steps[s.pos].Value
+			s.pos--
+			return &result, true
+		}
+
+		errn, ok := s.tryStep()
+		if !ok {
+			s.h /= 2
+			continue
+		}
+		tol := s.tolerance * s.h
+		if errn > tol {
+			s.steps[4] = s.steps[1]
+			s.adjustStep(errn)
+			if !s.tryInit() {
+				return nil, false
+			}
+			continue
+		}
+		s.lasttime += s.h
+		Shift(s.steps[:])
+		result.Time = s.lasttime
+		result.Value = s.steps[1].Value
+		if 10*errn < tol {
+			s.steps[4] = s.steps[1]
+			s.adjustStep(errn)
+			s.pos = 5
+		}
+		return &result, true
+	}
+	return nil, false
+}
+
+type adaptivePCMethod struct {
+	h         float64
+	hmax      float64
+	hmin      float64
+	tolerance float64
+}
+
+func (m *adaptivePCMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper {
+	result := adaptivePCStepper{
+		h:         h,
+		hmax:      m.hmax,
+		hmin:      m.hmin,
+		tolerance: m.tolerance,
+		lasttime:  ivp.Start.Time,
+		deriv:     ivp.Derivative,
+		pos:       5,
+	}
+	result.steps[4].Value = make([]float64, len(ivp.Start.Value))
+	copy(result.steps[4].Value, ivp.Start.Value)
+	return &result
+}
+
+func (m *adaptivePCMethod) Forward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(m.h, ivp)
+}
+
+func (m *adaptivePCMethod) Backward(ivp *mned.IVP) mned.Stepper {
+	return m.withStep(-m.h, ivp)
+}
+
+func AdaptivePredictorCorrector(
+	tolerance, step, hmin, hmax float64,
+) mned.Method {
+	return &adaptivePCMethod{
+		h:         step,
+		hmax:      hmax,
+		hmin:      hmin,
+		tolerance: tolerance,
+	}
+}
diff --git a/mstep/root.go b/mstep/root.go
new file mode 100644
index 0000000..20758f1
--- /dev/null
+++ b/mstep/root.go
@@ -0,0 +1,71 @@
+// Root-finding algorithms.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mstep
+
+import "math"
+
+func NewtonRoot1D(
+	x0, tolerance float64, f, df func(float64) (float64, bool),
+) (float64, bool) {
+	prev := x0
+	for {
+		fprev, ok := f(prev)
+		if !ok {
+			return 0, false
+		}
+		dfprev, ok := df(prev)
+		if !ok {
+			return 0, false
+		}
+		off := fprev / dfprev
+		next := prev - off
+		if math.Abs(off) <= tolerance {
+			return next, true
+		}
+		prev = next
+	}
+}
+
+func SecantRoot1D(
+	x0, x1, tolerance float64, f func(float64) (float64, bool),
+) (float64, bool) {
+	prev := x0
+	fprev, ok := f(prev)
+	if !ok {
+		return 0, false
+	}
+	cur := x1
+	fcur, ok := f(cur)
+	if !ok {
+		return 0, false
+	}
+	for {
+		off := fcur * (prev - cur) / (fprev - fcur)
+		next := cur - off
+		if math.Abs(off) <= tolerance {
+			return next, true
+		}
+		fnext, ok := f(next)
+		if !ok {
+			return 0, false
+		}
+		prev, fprev, cur, fcur = cur, fcur, next, fnext
+	}
+}
diff --git a/step.go b/step.go
new file mode 100644
index 0000000..b4a3657
--- /dev/null
+++ b/step.go
@@ -0,0 +1,324 @@
+// Stepping routines and method interface.
+//
+// Copyright (C) 2020 Juan Marín Noguera
+//
+// This file is part of Solvned.
+//
+// Solvned is free software: you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any 
+// later version.
+//
+// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY 
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with Solvned. If not, see <https://www.gnu.org/licenses/>.
+
+package mned
+
+import (
+	"math"
+	"sort"
+)
+
+// A Stepper is an iterator over points of the solution of an IVP. ODE solver
+// methods generally produce steppers that go from the initial values of the
+// problems to lower or greater values.
+type Stepper interface {
+	// Go to the next point and return a pointer to it. If the second value
+	// is false, there are no more points in the direction the Stepper is
+	// following and the value of Point should not be trusted.
+	//
+	// The point returned might change after the next call to next() on the
+	// same Stepper.
+	Next() (*Point, bool)
+}
+
+// Starting in the point init, step until the domain of the function ends or
+// one of the events says that the stepping should stop, calculating the
+// intermediante points with the given interpolator and, if `callback` is not
+// nil, calling the given callback with each new point. If ok is false, the
+// iteration ended because the stepper returned false; otherwise index contains
+// the index of the event that made the iteration stop.
+func StepUntil(
+	s Stepper,
+	init Point,
+	interp Interpolator,
+	callback func(*Point),
+	evs ...Event,
+) (index int, ok bool) {
+	var prev Point = init
+	vs := make([]float64, len(evs))
+	for i := 0; i < len(evs); i++ {
+		vs[i] = evs[i].Cross(&prev)
+	}
+	for {
+		point, ok := s.Next()
+		if !ok {
+			return 0, false
+		}
+		actions := make(requiredActions, 0)
+		for i := 0; i < len(evs); i++ {
+			newVal := evs[i].Cross(point)
+			if differentSign(vs[i], newVal) {
+				actions = append(actions, requiredAction{
+					index: i,
+					point: evs[i].FindPoint(
+						interp, &prev, point,
+					)})
+			}
+			vs[i] = newVal
+		}
+
+		sort.Sort(actions)
+		size := len(actions)
+		if point.Time < prev.Time {
+			for i := 0; i < size/2; i++ {
+				actions[i], actions[size-i] =
+					actions[size-i], actions[i]
+			}
+		}
+		for _, a := range actions {
+			if !evs[a.index].Action(&a.point) {
+				return a.index, true
+			}
+		}
+
+		callback(point)
+		prev = point.Clone()
+	}
+}
+
+// A ConfigurableStepper is a Stepper which allows specifying the delta of the
+// next step; that is, the difference between the time of a point and that of
+// the previous point.
+type ConfigurableStepper interface {
+	Stepper
+	// Like Stepper.Next but the time of the next point is specified to be
+	// the time of the latest point plus the value given.
+	NextStep(float64) (*Point, bool)
+}
+
+// A Method is a way of approximating points of the solution to an IVP. It works
+// as a factory of Steppers that step forward and backward from the initial
+// values of the problem.
+type Method interface {
+	// Make a stepper over the solution of the given IVP that starts at the
+	// initial values of the problem and goes to ever increasing values for
+	// the independent variable.
+	Forward(*IVP) Stepper
+	// Make a stepper over the solution of the given IVP that starts at the
+	// initial values of the problem and goes to ever decreasing values for
+	// the independent variable.
+	Backward(*IVP) Stepper
+}
+
+type fixedStepper struct {
+	state  IVP
+	step   float64
+	method func(*IVP, float64) bool
+}
+
+func (s *fixedStepper) NextStep(h float64) (*Point, bool) {
+	ok := s.method(&s.state, h)
+	if !ok {
+		return nil, false
+	}
+	s.state.Start.Time += h
+	return &s.state.Start, true
+}
+
+func (s *fixedStepper) Next() (*Point, bool) {
+	return s.NextStep(s.step)
+}
+
+// Create a fixed step method, one where all steps increment or decrement the
+// time by the same amount, with the added condition that it shouldn't depend
+// on extra state. The Step is the fixed increment/decrement of time, which
+// must be positive. The steppers returned by this method are
+// `ConfigurableStepper`s.
+//
+// The Method is a function that takes an IVP and a step and stores in
+// IVP.Start.Value the value for IVP.Start.Time+step, while leaving
+// IVP.Start.Time untouched. The return value of `method` is true if the
+// operation completed successfully or false if there was a problem, such as
+// going out of the domain of the IVP.Derivative, in which case the resulting
+// state of the IVP is unspecified.
+type FixedStepMethod struct {
+	Step   float64
+	Method func(*IVP, float64) bool
+}
+
+func (m FixedStepMethod) Forward(ivp *IVP) Stepper {
+	return &fixedStepper{
+		state:  ivp.Clone(),
+		step:   m.Step,
+		method: m.Method}
+}
+
+func (m FixedStepMethod) Backward(ivp *IVP) Stepper {
+	return &fixedStepper{
+		state:  ivp.Clone(),
+		step:   -m.Step,
+		method: m.Method}
+}
+
+type adaptiveStepper struct {
+	state     IVP
+	step      float64
+	hmin      float64
+	hmax      float64
+	tolerance float64
+	method    func(*IVP, float64) bool
+	order     uint8
+}
+
+// Return the new step that should be taken by an AdaptiveStepMethod of the
+// given order and with a given tolerance took a step of a given size yielding
+// the given amount of error, with max as the maximum size of a step.
+func AdaptiveUpdateStep(
+	step float64, tol float64, err float64, max float64, order uint8,
+) float64 {
+	var q float64
+	adjust := tol * math.Abs(step) / (2 * err)
+	if order == 1 { // Fast path
+		q = adjust
+	} else {
+		q = math.Pow(adjust, 1/float64(order))
+	}
+	switch {
+	case q < 0.1:
+		q = 0.1
+	case q > 4:
+		q = 4
+	}
+	newStep := q * step
+	if math.Abs(newStep) > max {
+		if newStep > 0 {
+			newStep = max
+		} else {
+			newStep = -max
+		}
+	}
+	return newStep
+}
+
+func (s *adaptiveStepper) adaptStep(err float64) {
+	s.step = AdaptiveUpdateStep(s.step, s.tolerance, err, s.hmax, s.order)
+}
+
+func (s *adaptiveStepper) getSteps(full []float64, half []float64) bool {
+	orig := s.state.Start.Value
+	copy(full, orig)
+	s.state.Start.Value = full
+	if !s.method(&s.state, s.step) {
+		return false
+	}
+	copy(half, orig)
+	s.state.Start.Value = half
+	if !s.method(&s.state, s.step/2) {
+		return false
+	}
+	if !s.method(&s.state, s.step/2) {
+		return false
+	}
+	s.state.Start.Value = orig
+	return true
+}
+
+func (s *adaptiveStepper) Next() (*Point, bool) {
+	size := len(s.state.Start.Value)
+	full := make([]float64, size)
+	half := make([]float64, size)
+
+	for s.step >= s.hmin {
+		for !s.getSteps(full, half) {
+			if s.step *= 0.5; s.step < s.hmin {
+				return nil, false
+			}
+		}
+
+		orderm1 := 1 / (math.Ldexp(1, int(s.order)) - 1)
+		multiplier := 1 + orderm1
+		sqerr := float64(0)
+		for i, v := range full {
+			diff := v - half[i]
+			sqerr += diff * diff
+		}
+		err := math.Sqrt(sqerr) * multiplier
+		tol := s.tolerance * math.Abs(s.step)
+		if err > tol {
+			s.adaptStep(err)
+			continue
+		}
+
+		for i, v := range full {
+			s.state.Start.Value[i] =
+				multiplier*half[i] - orderm1*v
+		}
+
+		s.state.Start.Time += s.step
+		s.adaptStep(err)
+		return &s.state.Start, true
+	}
+
+	return nil, false // s.step < s.hmin
+}
+
+// An AdaptiveStepMethod is a kind of adaptive method (that is, one that changes
+// the step over time for greater precision and performance) that can be adapted
+// from a fixed-step method.
+//
+// Step, Hmin, and Hmax are the initial, minimum, and maximum step. If the
+// step must go lower than Hmin to satisfy the tolerance, the stepper ends;
+// if it could get greater than Hmax, it saturates to Hmax. The Tolerance is
+// such that the Euclidean norm of the estimated error of any step of size `h`
+// must be lower than `h*Tolerance`.
+//
+// Method works the same way as in `FixedStepMethod` except that its error for a
+// given interval of time must be on the order of `O(h^Order)`, with the
+// variable `h` being the step. Order must be positive.
+//
+// The algorithm used is the following, where `(t,x(t))` are the initial values,
+// `y(h,t,x(t))` is the result of applying the inner method with step `h` from
+// `(t,x(t))`, `k` is the order of the method and `h` is the step.
+// 1. First, `F:=y(h,t,x(t))` and `H:=y(h/2,t+h/2,y(h/2,t,x(t)))` are obtained,
+//    and the error is estimated as `e:=|(1+1/(2^k-1))*(H-F)|`.
+// 2. If that error is lower than `|h|*Tolerance`, the step is accepted.
+//    Otherwise we set `q:=(|h|*Tolerance)/(2*e)`, saturate so that `q` is in
+//    `[0.1, 4]`, set `h:=q*h`, check bounds on `h` and return to the start.
+// 3. The next value is taken to be `(t+h, (2^k*H - H)/(2^k - 1))`, and the
+//    next step size is calculated by changing `h` as in step 2.
+//
+// The steppers are *not* `ConfigurableStepper`s.
+type AdaptiveStepMethod struct {
+	Step      float64
+	Hmin      float64
+	Hmax      float64
+	Tolerance float64
+	Method    func(*IVP, float64) bool
+	Order     uint8
+}
+
+func (m *AdaptiveStepMethod) toStepper(ivp *IVP, step float64) Stepper {
+	return &adaptiveStepper{
+		state:     ivp.Clone(),
+		step:      step,
+		hmin:      m.Hmin,
+		hmax:      m.Hmax,
+		tolerance: m.Tolerance,
+		method:    m.Method,
+		order:     m.Order,
+	}
+}
+
+func (m *AdaptiveStepMethod) Forward(ivp *IVP) Stepper {
+	return m.toStepper(ivp, m.Step)
+}
+
+func (m *AdaptiveStepMethod) Backward(ivp *IVP) Stepper {
+	return m.toStepper(ivp, -m.Step)
+}