diff --git a/diff/.travis.yml b/diff/.travis.yml
new file mode 100644
index 00000000..e0b47f93
--- /dev/null
+++ b/diff/.travis.yml
@@ -0,0 +1,21 @@
+language: go
+
+# Versions of go that are explicitly supported by gonum.
+go:
+ - 1.5.4
+ - 1.6.3
+ - 1.7.3
+
+# Required for coverage.
+before_install:
+ - go get golang.org/x/tools/cmd/cover
+ - go get github.com/mattn/goveralls
+
+# Get deps, build, test, and ensure the code is gofmt'ed.
+# If we are building as gonum, then we have access to the coveralls api key, so we can run coverage as well.
+script:
+ - go get -d -t -v ./...
+ - go build -v ./...
+ - go test -v ./...
+ - test -z "$(gofmt -d .)"
+ - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi
diff --git a/diff/.travis/test-coverage.sh b/diff/.travis/test-coverage.sh
new file mode 100755
index 00000000..7df8aa6a
--- /dev/null
+++ b/diff/.travis/test-coverage.sh
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+PROFILE_OUT=$PWD/profile.out
+ACC_OUT=$PWD/acc.out
+
+testCover() {
+	# set the return value to 0 (succesful)
+	retval=0
+	# get the directory to check from the parameter. Default to '.'
+	d=${1:-.}
+	# skip if there are no Go files here
+	ls $d/*.go &> /dev/null || return $retval
+	# switch to the directory to check
+	pushd $d > /dev/null
+	# create the coverage profile
+	coverageresult=`go test -v -coverprofile=$PROFILE_OUT`
+	# output the result so we can check the shell output
+	echo ${coverageresult}
+	# append the results to acc.out if coverage didn't fail, else set the retval to 1 (failed)
+	( [[ ${coverageresult} == *FAIL* ]] && retval=1 ) || ( [ -f $PROFILE_OUT ] && grep -v "mode: set" $PROFILE_OUT >> $ACC_OUT )
+	# return to our working dir
+	popd > /dev/null
+	# return our return value
+	return $retval
+}
+
+# Init acc.out
+echo "mode: set" > $ACC_OUT
+
+# Run test coverage on all directories containing go files
+find . -maxdepth 10 -type d | while read d; do testCover $d || exit; done
+
+# Upload the coverage profile to coveralls.io
+[ -n "$COVERALLS_TOKEN" ] && goveralls -coverprofile=$ACC_OUT -service=travis-ci -repotoken $COVERALLS_TOKEN
+
diff --git a/diff/README.md b/diff/README.md
new file mode 100644
index 00000000..03ba1b8c
--- /dev/null
+++ b/diff/README.md
@@ -0,0 +1,13 @@
+# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://coveralls.io/repos/gonum/diff/badge.svg?branch=master&service=github)](https://coveralls.io/github/gonum/diff?branch=master) [![GoDoc](https://godoc.org/github.com/gonum/diff?status.svg)](https://godoc.org/github.com/gonum/diff)
+
+This is a package for computing derivatives of functions for the Go language.
+
+## Issues
+
+If you find any bugs, feel free to file an issue on the github issue tracker. Discussions on API changes, added features, code review, or similar requests are preferred on the gonum-dev Google Group.
+
+https://groups.google.com/forum/#!forum/gonum-dev
+
+## License
+
+Please see github.com/gonum/license for general license information, contributors, authors, etc on the Gonum suite of packages.
diff --git a/diff/fd/diff.go b/diff/fd/diff.go
new file mode 100644
index 00000000..3fb45b0a
--- /dev/null
+++ b/diff/fd/diff.go
@@ -0,0 +1,276 @@
+// Copyright ©2014 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package fd provides functions to approximate derivatives using finite differences.
+package fd
+
+import (
+	"math"
+	"runtime"
+	"sync"
+
+	"github.com/gonum/floats"
+)
+
+// A Point is a stencil location in a finite difference formula.
+type Point struct {
+	Loc   float64
+	Coeff float64
+}
+
+// Formula represents a finite difference formula on a regularly spaced grid
+// that approximates the derivative of order k of a function f at x as
+//  d^k f(x) ≈ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
+type Formula struct {
+	// Stencil is the set of sampling Points which are used to estimate the
+	// derivative. The locations will be scaled by Step and are relative to x.
+	Stencil    []Point
+	Derivative int     // The order of the approximated derivative.
+	Step       float64 // Default step size for the formula.
+}
+
+func (f Formula) isZero() bool {
+	return f.Stencil == nil && f.Derivative == 0 && f.Step == 0
+}
+
+// Settings is the settings structure for computing finite differences.
+type Settings struct {
+	// Formula is the finite difference formula used
+	// for approximating the derivative.
+	// Zero value indicates a default formula.
+	Formula Formula
+	// Step is the distance between points of the stencil.
+	// If equal to 0, formula's default step will be used.
+	Step float64
+
+	OriginKnown bool    // Flag that the value at the origin x is known.
+	OriginValue float64 // Value at the origin (only used if OriginKnown is true).
+
+	Concurrent bool // Should the function calls be executed concurrently.
+}
+
+// Derivative estimates the derivative of the function f at the given location.
+// The finite difference formula, the step size, and other options are
+// specified by settings. If settings is nil, the first derivative will be
+// estimated using the Forward formula and a default step size.
+func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
+	if settings == nil {
+		settings = &Settings{}
+	}
+	formula := settings.Formula
+	if formula.isZero() {
+		formula = Forward
+	}
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("fd: bad formula")
+	}
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
+	var deriv float64
+	if !settings.Concurrent || runtime.GOMAXPROCS(0) == 1 {
+		for _, pt := range formula.Stencil {
+			if settings.OriginKnown && pt.Loc == 0 {
+				deriv += pt.Coeff * settings.OriginValue
+				continue
+			}
+			deriv += pt.Coeff * f(x+step*pt.Loc)
+		}
+		return deriv / math.Pow(step, float64(formula.Derivative))
+	}
+
+	wg := &sync.WaitGroup{}
+	mux := &sync.Mutex{}
+	for _, pt := range formula.Stencil {
+		if settings.OriginKnown && pt.Loc == 0 {
+			mux.Lock()
+			deriv += pt.Coeff * settings.OriginValue
+			mux.Unlock()
+			continue
+		}
+		wg.Add(1)
+		go func(pt Point) {
+			defer wg.Done()
+			fofx := f(x + step*pt.Loc)
+			mux.Lock()
+			defer mux.Unlock()
+			deriv += pt.Coeff * fofx
+		}(pt)
+	}
+	wg.Wait()
+	return deriv / math.Pow(step, float64(formula.Derivative))
+}
+
+// Gradient estimates the gradient of the multivariate function f at the
+// location x. If dst is not nil, the result will be stored in-place into dst
+// and returned, otherwise a new slice will be allocated first. Finite
+// difference kernel and other options are specified by settings. If settings is
+// nil, the gradient will be estimated using the Forward formula and a default
+// step size.
+//
+// Gradient panics if the length of dst and x is not equal, or if the derivative
+// order of the formula is not 1.
+func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
+	if dst == nil {
+		dst = make([]float64, len(x))
+	}
+	if len(dst) != len(x) {
+		panic("fd: slice length mismatch")
+	}
+	if settings == nil {
+		settings = &Settings{}
+	}
+
+	formula := settings.Formula
+	if formula.isZero() {
+		formula = Forward
+	}
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("fd: bad formula")
+	}
+	if formula.Derivative != 1 {
+		panic("fd: invalid derivative order")
+	}
+
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
+	expect := len(formula.Stencil) * len(x)
+	nWorkers := 1
+	if settings.Concurrent {
+		nWorkers = runtime.GOMAXPROCS(0)
+		if nWorkers > expect {
+			nWorkers = expect
+		}
+	}
+
+	var hasOrigin bool
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			hasOrigin = true
+			break
+		}
+	}
+	xcopy := make([]float64, len(x)) // So that x is not modified during the call.
+	originValue := settings.OriginValue
+	if hasOrigin && !settings.OriginKnown {
+		copy(xcopy, x)
+		originValue = f(xcopy)
+	}
+
+	if nWorkers == 1 {
+		for i := range xcopy {
+			var deriv float64
+			for _, pt := range formula.Stencil {
+				if pt.Loc == 0 {
+					deriv += pt.Coeff * originValue
+					continue
+				}
+				copy(xcopy, x)
+				xcopy[i] += pt.Loc * step
+				deriv += pt.Coeff * f(xcopy)
+			}
+			dst[i] = deriv / step
+		}
+		return dst
+	}
+
+	sendChan := make(chan fdrun, expect)
+	ansChan := make(chan fdrun, expect)
+	quit := make(chan struct{})
+	defer close(quit)
+
+	// Launch workers. Workers receive an index and a step, and compute the answer.
+	for i := 0; i < nWorkers; i++ {
+		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
+			xcopy := make([]float64, len(x))
+			for {
+				select {
+				case <-quit:
+					return
+				case run := <-sendChan:
+					copy(xcopy, x)
+					xcopy[run.idx] += run.pt.Loc * step
+					run.result = f(xcopy)
+					ansChan <- run
+				}
+			}
+		}(sendChan, ansChan, quit)
+	}
+
+	// Launch the distributor. Distributor sends the cases to be computed.
+	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
+		for i := range x {
+			for _, pt := range formula.Stencil {
+				if pt.Loc == 0 {
+					// Answer already known. Send the answer on the answer channel.
+					ansChan <- fdrun{
+						idx:    i,
+						pt:     pt,
+						result: originValue,
+					}
+					continue
+				}
+				// Answer not known, send the answer to be computed.
+				sendChan <- fdrun{
+					idx: i,
+					pt:  pt,
+				}
+			}
+		}
+	}(sendChan, ansChan)
+
+	for i := range dst {
+		dst[i] = 0
+	}
+	// Read in all of the results.
+	for i := 0; i < expect; i++ {
+		run := <-ansChan
+		dst[run.idx] += run.pt.Coeff * run.result
+	}
+	floats.Scale(1/step, dst)
+	return dst
+}
+
+type fdrun struct {
+	idx    int
+	pt     Point
+	result float64
+}
+
+// Forward represents a first-order accurate forward approximation
+// to the first derivative.
+var Forward = Formula{
+	Stencil:    []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
+	Derivative: 1,
+	Step:       2e-8,
+}
+
+// Backward represents a first-order accurate backward approximation
+// to the first derivative.
+var Backward = Formula{
+	Stencil:    []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
+	Derivative: 1,
+	Step:       2e-8,
+}
+
+// Central represents a second-order accurate centered approximation
+// to the first derivative.
+var Central = Formula{
+	Stencil:    []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
+	Derivative: 1,
+	Step:       6e-6,
+}
+
+// Central2nd represents a secord-order accurate centered approximation
+// to the second derivative.
+var Central2nd = Formula{
+	Stencil:    []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
+	Derivative: 2,
+	Step:       1e-4,
+}
diff --git a/diff/fd/diff_test.go b/diff/fd/diff_test.go
new file mode 100644
index 00000000..35d4f0f3
--- /dev/null
+++ b/diff/fd/diff_test.go
@@ -0,0 +1,145 @@
+// Copyright ©2014 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package fd
+
+import (
+	"math"
+	"testing"
+)
+
+var xSquared = func(x float64) float64 { return x * x }
+
+type testPoint struct {
+	f    func(float64) float64
+	loc  float64
+	fofx float64
+	ans  float64
+}
+
+var testsFirst = []testPoint{
+	{
+		f:    xSquared,
+		loc:  0,
+		fofx: 0,
+		ans:  0,
+	},
+	{
+		f:    xSquared,
+		loc:  5,
+		fofx: 25,
+		ans:  10,
+	},
+	{
+		f:    xSquared,
+		loc:  2,
+		fofx: 4,
+		ans:  4,
+	},
+	{
+		f:    xSquared,
+		loc:  -5,
+		fofx: 25,
+		ans:  -10,
+	},
+}
+
+var testsSecond = []testPoint{
+	{
+		f:    xSquared,
+		loc:  0,
+		fofx: 0,
+		ans:  2,
+	},
+	{
+		f:    xSquared,
+		loc:  5,
+		fofx: 25,
+		ans:  2,
+	},
+	{
+		f:    xSquared,
+		loc:  2,
+		fofx: 4,
+		ans:  2,
+	},
+	{
+		f:    xSquared,
+		loc:  -5,
+		fofx: 25,
+		ans:  2,
+	},
+}
+
+func testDerivative(t *testing.T, formula Formula, tol float64, tests []testPoint) {
+	for i, test := range tests {
+
+		ans := Derivative(test.f, test.loc, &Settings{
+			Formula: formula,
+		})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch serial: expected %v, found %v", i, test.ans, ans)
+		}
+
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:     formula,
+			OriginKnown: true,
+			OriginValue: test.fofx,
+		})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch serial origin known: expected %v, found %v", i, test.ans, ans)
+		}
+
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:    formula,
+			Concurrent: true,
+		})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
+		}
+
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:     formula,
+			OriginKnown: true,
+			OriginValue: test.fofx,
+			Concurrent:  true,
+		})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
+		}
+	}
+}
+
+func TestForward(t *testing.T) {
+	testDerivative(t, Forward, 2e-4, testsFirst)
+}
+
+func TestBackward(t *testing.T) {
+	testDerivative(t, Backward, 2e-4, testsFirst)
+}
+
+func TestCentral(t *testing.T) {
+	testDerivative(t, Central, 1e-6, testsFirst)
+}
+
+func TestCentralSecond(t *testing.T) {
+	testDerivative(t, Central2nd, 1e-3, testsSecond)
+}
+
+// TestDerivativeDefault checks that the derivative works when settings is nil
+// or zero value.
+func TestDerivativeDefault(t *testing.T) {
+	tol := 1e-6
+	for i, test := range testsFirst {
+		ans := Derivative(test.f, test.loc, nil)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch default: expected %v, found %v", i, test.ans, ans)
+		}
+
+		ans = Derivative(test.f, test.loc, &Settings{})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch zero value: expected %v, found %v", i, test.ans, ans)
+		}
+	}
+}
diff --git a/diff/fd/example_test.go b/diff/fd/example_test.go
new file mode 100644
index 00000000..1856ffa6
--- /dev/null
+++ b/diff/fd/example_test.go
@@ -0,0 +1,68 @@
+// Copyright ©2016 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package fd_test
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/gonum/diff/fd"
+	"github.com/gonum/matrix/mat64"
+)
+
+func ExampleDerivative() {
+	f := func(x float64) float64 {
+		return math.Sin(x)
+	}
+	// Compute the first derivative of f at 0 using the default settings.
+	fmt.Println("f'(0) ≈", fd.Derivative(f, 0, nil))
+	// Compute the first derivative of f at 0 using the forward approximation
+	// with a custom step size.
+	df := fd.Derivative(f, 0, &fd.Settings{
+		Formula: fd.Forward,
+		Step:    1e-3,
+	})
+	fmt.Println("f'(0) ≈", df)
+
+	f = func(x float64) float64 {
+		return math.Pow(math.Cos(x), 3)
+	}
+	// Compute the second derivative of f at 0 using
+	// the centered approximation, concurrent evaluation,
+	// and a known function value at x.
+	df = fd.Derivative(f, 0, &fd.Settings{
+		Formula:     fd.Central2nd,
+		Concurrent:  true,
+		OriginKnown: true,
+		OriginValue: f(0),
+	})
+	fmt.Println("f''(0) ≈", df)
+
+	// Output:
+	// f'(0) ≈ 1
+	// f'(0) ≈ 0.9999998333333416
+	// f''(0) ≈ -2.999999981767587
+}
+
+func ExampleJacobian() {
+	f := func(dst, x []float64) {
+		dst[0] = x[0] + 1
+		dst[1] = 5 * x[2]
+		dst[2] = 4*x[1]*x[1] - 2*x[2]
+		dst[3] = x[2] * math.Sin(x[0])
+	}
+	jac := mat64.NewDense(4, 3, nil)
+	fd.Jacobian(jac, f, []float64{1, 2, 3}, &fd.JacobianSettings{
+		Formula:    fd.Central,
+		Concurrent: true,
+	})
+	fmt.Printf("J ≈ %.6v\n", mat64.Formatted(jac, mat64.Prefix("    ")))
+
+	// Output:
+	// J ≈ ⎡       1         0         0⎤
+	//     ⎢       0         0         5⎥
+	//     ⎢       0        16        -2⎥
+	//     ⎣ 1.62091         0  0.841471⎦
+}
diff --git a/diff/fd/gradient_test.go b/diff/fd/gradient_test.go
new file mode 100644
index 00000000..09a335e8
--- /dev/null
+++ b/diff/fd/gradient_test.go
@@ -0,0 +1,187 @@
+// Copyright ©2014 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package fd
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+)
+
+type Rosenbrock struct {
+	nDim int
+}
+
+func (r Rosenbrock) F(x []float64) (sum float64) {
+	deriv := make([]float64, len(x))
+	return r.FDf(x, deriv)
+}
+
+func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
+	for i := range deriv {
+		deriv[i] = 0
+	}
+
+	for i := 0; i < len(x)-1; i++ {
+		sum += math.Pow(1-x[i], 2) + 100*math.Pow(x[i+1]-math.Pow(x[i], 2), 2)
+	}
+	for i := 0; i < len(x)-1; i++ {
+		deriv[i] += -1 * 2 * (1 - x[i])
+		deriv[i] += 2 * 100 * (x[i+1] - math.Pow(x[i], 2)) * (-2 * x[i])
+	}
+	for i := 1; i < len(x); i++ {
+		deriv[i] += 2 * 100 * (x[i] - math.Pow(x[i-1], 2))
+	}
+
+	return sum
+}
+
+func TestGradient(t *testing.T) {
+	rand.Seed(1)
+	for i, test := range []struct {
+		nDim    int
+		tol     float64
+		formula Formula
+	}{
+		{
+			nDim:    2,
+			tol:     2e-4,
+			formula: Forward,
+		},
+		{
+			nDim:    2,
+			tol:     1e-6,
+			formula: Central,
+		},
+		{
+			nDim:    40,
+			tol:     2e-4,
+			formula: Forward,
+		},
+		{
+			nDim:    40,
+			tol:     1e-6,
+			formula: Central,
+		},
+	} {
+		x := make([]float64, test.nDim)
+		for i := range x {
+			x[i] = rand.Float64()
+		}
+		xcopy := make([]float64, len(x))
+		copy(xcopy, x)
+
+		r := Rosenbrock{len(x)}
+		trueGradient := make([]float64, len(x))
+		r.FDf(x, trueGradient)
+
+		// Try with gradient nil.
+		gradient := Gradient(nil, r.F, x, &Settings{
+			Formula: test.formula,
+		})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch in serial with nil. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in serial with nil.", i)
+		}
+
+		// Try with provided gradient.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{
+			Formula: test.formula,
+		})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in serial with non-nil.", i)
+		}
+
+		// Try with known value.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:     test.formula,
+			OriginKnown: true,
+			OriginValue: r.F(x),
+		})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with known origin in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+		// Try with concurrent evaluation.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:    test.formula,
+			Concurrent: true,
+		})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in parallel", i)
+		}
+
+		// Try with concurrent evaluation with origin known.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:     test.formula,
+			Concurrent:  true,
+			OriginKnown: true,
+			OriginValue: r.F(x),
+		})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+		// Try with nil settings.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, nil)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+		// Try with zero-valued settings.
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with zero settings. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+	}
+}
+
+func Panics(fun func()) (b bool) {
+	defer func() {
+		err := recover()
+		if err != nil {
+			b = true
+		}
+	}()
+	fun()
+	return
+}
+
+func TestGradientPanics(t *testing.T) {
+	// Test that it panics
+	if !Panics(func() {
+		Gradient([]float64{0.0}, func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil)
+	}) {
+		t.Errorf("Gradient did not panic with length mismatch")
+	}
+}
diff --git a/diff/fd/jacobian.go b/diff/fd/jacobian.go
new file mode 100644
index 00000000..1aab264c
--- /dev/null
+++ b/diff/fd/jacobian.go
@@ -0,0 +1,203 @@
+// Copyright ©2016 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package fd
+
+import (
+	"runtime"
+	"sync"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+type JacobianSettings struct {
+	Formula     Formula
+	OriginValue []float64
+	Step        float64
+	Concurrent  bool
+}
+
+// Jacobian approximates the Jacobian matrix of a vector-valued function f at
+// the location x and stores the result in-place into dst.
+//
+// Finite difference formula and other options are specified by settings. If
+// settings is nil, the Jacobian will be estimated using the Forward formula and
+// a default step size.
+//
+// The Jacobian matrix J is the matrix of all first-order partial derivatives of f.
+// If f maps an n-dimensional vector x to an m-dimensional vector y = f(x), J is
+// an m×n matrix whose elements are given as
+//  J_{i,j} = ∂f_i/∂x_j,
+// or expanded out
+//      [ ∂f_1/∂x_1 ... ∂f_1/∂x_n ]
+//      [     .  .          .     ]
+//  J = [     .      .      .     ]
+//      [     .          .  .     ]
+//      [ ∂f_m/∂x_1 ... ∂f_m/∂x_n ]
+//
+// dst must be non-nil, the number of its columns must equal the length of x, and
+// the derivative order of the formula must be 1, otherwise Jacobian will panic.
+func Jacobian(dst *mat64.Dense, f func(y, x []float64), x []float64, settings *JacobianSettings) {
+	n := len(x)
+	if n == 0 {
+		panic("jacobian: x has zero length")
+	}
+	m, c := dst.Dims()
+	if c != n {
+		panic("jacobian: mismatched matrix size")
+	}
+
+	if settings == nil {
+		settings = &JacobianSettings{}
+	}
+	if settings.OriginValue != nil && len(settings.OriginValue) != m {
+		panic("jacobian: mismatched OriginValue slice length")
+	}
+
+	formula := settings.Formula
+	if formula.isZero() {
+		formula = Forward
+	}
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("jacobian: bad formula")
+	}
+	if formula.Derivative != 1 {
+		panic("jacobian: invalid derivative order")
+	}
+
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
+	evals := n * len(formula.Stencil)
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			evals -= n - 1
+			break
+		}
+	}
+	nWorkers := 1
+	if settings.Concurrent {
+		nWorkers = runtime.GOMAXPROCS(0)
+		if nWorkers > evals {
+			nWorkers = evals
+		}
+	}
+	if nWorkers == 1 {
+		jacobianSerial(dst, f, x, settings.OriginValue, formula, step)
+	} else {
+		jacobianConcurrent(dst, f, x, settings.OriginValue, formula, step, nWorkers)
+	}
+}
+
+func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64) {
+	m, n := dst.Dims()
+	xcopy := make([]float64, n)
+	y := make([]float64, m)
+	col := make([]float64, m)
+	for j := 0; j < n; j++ {
+		for i := range col {
+			col[i] = 0
+		}
+		for _, pt := range formula.Stencil {
+			if pt.Loc == 0 {
+				if origin == nil {
+					origin = make([]float64, m)
+					copy(xcopy, x)
+					f(origin, xcopy)
+				}
+				floats.AddScaled(col, pt.Coeff, origin)
+			} else {
+				copy(xcopy, x)
+				xcopy[j] += pt.Loc * step
+				f(y, xcopy)
+				floats.AddScaled(col, pt.Coeff, y)
+			}
+		}
+		dst.SetCol(j, col)
+	}
+	dst.Scale(1/step, dst)
+}
+
+func jacobianConcurrent(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64, nWorkers int) {
+	m, n := dst.Dims()
+	for i := 0; i < m; i++ {
+		for j := 0; j < n; j++ {
+			dst.Set(i, j, 0)
+		}
+	}
+
+	var (
+		wg sync.WaitGroup
+		mu = make([]sync.Mutex, n) // Guard access to individual columns.
+	)
+	worker := func(jobs <-chan jacJob) {
+		defer wg.Done()
+		xcopy := make([]float64, n)
+		y := make([]float64, m)
+		yVec := mat64.NewVector(m, y)
+		for job := range jobs {
+			copy(xcopy, x)
+			xcopy[job.j] += job.pt.Loc * step
+			f(y, xcopy)
+			col := dst.ColView(job.j)
+			mu[job.j].Lock()
+			col.AddScaledVec(col, job.pt.Coeff, yVec)
+			mu[job.j].Unlock()
+		}
+	}
+	jobs := make(chan jacJob, nWorkers)
+	for i := 0; i < nWorkers; i++ {
+		wg.Add(1)
+		go worker(jobs)
+	}
+	var hasOrigin bool
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			hasOrigin = true
+			continue
+		}
+		for j := 0; j < n; j++ {
+			jobs <- jacJob{j, pt}
+		}
+	}
+	close(jobs)
+	if hasOrigin && origin == nil {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			origin = make([]float64, m)
+			xcopy := make([]float64, n)
+			copy(xcopy, x)
+			f(origin, xcopy)
+		}()
+	}
+	wg.Wait()
+
+	if hasOrigin {
+		// The formula evaluated at x, we need to add scaled origin to
+		// all columns of dst. Iterate again over all Formula points
+		// because we don't forbid repeated locations.
+
+		originVec := mat64.NewVector(m, origin)
+		for _, pt := range formula.Stencil {
+			if pt.Loc != 0 {
+				continue
+			}
+			for j := 0; j < n; j++ {
+				col := dst.ColView(j)
+				col.AddScaledVec(col, pt.Coeff, originVec)
+			}
+		}
+	}
+
+	dst.Scale(1/step, dst)
+}
+
+type jacJob struct {
+	j  int
+	pt Point
+}
diff --git a/diff/fd/jacobian_test.go b/diff/fd/jacobian_test.go
new file mode 100644
index 00000000..89f1af11
--- /dev/null
+++ b/diff/fd/jacobian_test.go
@@ -0,0 +1,268 @@
+// Copyright ©2016 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package fd
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+func vecFunc13(y, x []float64) {
+	y[0] = 5*x[0] + x[2]*math.Sin(x[1]) + 1
+}
+func vecFunc13Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 5)
+	jac.Set(0, 1, x[2]*math.Cos(x[1]))
+	jac.Set(0, 2, math.Sin(x[1]))
+}
+
+func vecFunc22(y, x []float64) {
+	y[0] = x[0]*x[0]*x[1] + 1
+	y[1] = 5*x[0] + math.Sin(x[1]) + 1
+}
+func vecFunc22Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 2*x[0]*x[1])
+	jac.Set(0, 1, x[0]*x[0])
+	jac.Set(1, 0, 5)
+	jac.Set(1, 1, math.Cos(x[1]))
+}
+
+func vecFunc43(y, x []float64) {
+	y[0] = x[0] + 1
+	y[1] = 5*x[2] + 1
+	y[2] = 4*x[1]*x[1] - 2*x[2] + 1
+	y[3] = x[2]*math.Sin(x[0]) + 1
+}
+func vecFunc43Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 1)
+	jac.Set(0, 1, 0)
+	jac.Set(0, 2, 0)
+	jac.Set(1, 0, 0)
+	jac.Set(1, 1, 0)
+	jac.Set(1, 2, 5)
+	jac.Set(2, 0, 0)
+	jac.Set(2, 1, 8*x[1])
+	jac.Set(2, 2, -2)
+	jac.Set(3, 0, x[2]*math.Cos(x[0]))
+	jac.Set(3, 1, 0)
+	jac.Set(3, 2, math.Sin(x[0]))
+}
+
+func TestJacobian(t *testing.T) {
+	rand.Seed(1)
+
+	// Test with default settings.
+	for tc, test := range []struct {
+		m, n int
+		f    func([]float64, []float64)
+		jac  func(*mat64.Dense, []float64)
+	}{
+		{
+			m:   1,
+			n:   3,
+			f:   vecFunc13,
+			jac: vecFunc13Jac,
+		},
+		{
+			m:   2,
+			n:   2,
+			f:   vecFunc22,
+			jac: vecFunc22Jac,
+		},
+		{
+			m:   4,
+			n:   3,
+			f:   vecFunc43,
+			jac: vecFunc43Jac,
+		},
+	} {
+		const tol = 1e-6
+
+		x := randomSlice(test.n, 10)
+		xcopy := make([]float64, test.n)
+		copy(xcopy, x)
+
+		want := mat64.NewDense(test.m, test.n, nil)
+		test.jac(want, x)
+
+		got := mat64.NewDense(test.m, test.n, nil)
+		fillNaNDense(got)
+		Jacobian(got, test.f, x, nil)
+		if !mat64.EqualApprox(want, got, tol) {
+			t.Errorf("Case %d (default settings): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (default settings): x modified", tc)
+		}
+	}
+
+	// Test with non-default settings.
+	for tc, test := range []struct {
+		m, n    int
+		f       func([]float64, []float64)
+		jac     func(*mat64.Dense, []float64)
+		tol     float64
+		formula Formula
+	}{
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+	} {
+		x := randomSlice(test.n, 10)
+		xcopy := make([]float64, test.n)
+		copy(xcopy, x)
+
+		want := mat64.NewDense(test.m, test.n, nil)
+		test.jac(want, x)
+
+		got := mat64.NewDense(test.m, test.n, nil)
+		fillNaNDense(got)
+		Jacobian(got, test.f, x, &JacobianSettings{
+			Formula: test.formula,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d: unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d: x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, x, &JacobianSettings{
+			Formula:    test.formula,
+			Concurrent: true,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (concurrent): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (concurrent): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		origin := make([]float64, test.m)
+		test.f(origin, x)
+		Jacobian(got, test.f, x, &JacobianSettings{
+			Formula:     test.formula,
+			OriginValue: origin,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (origin): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (origin): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, x, &JacobianSettings{
+			Formula:     test.formula,
+			OriginValue: origin,
+			Concurrent:  true,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (concurrent, origin): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (concurrent, origin): x modified", tc)
+		}
+	}
+}
+
+// randomSlice returns a slice of n elements from the interval [-bound,bound).
+func randomSlice(n int, bound float64) []float64 {
+	x := make([]float64, n)
+	for i := range x {
+		x[i] = 2*bound*rand.Float64() - bound
+	}
+	return x
+}
+
+// fillNaNDense fills the matrix m with NaN values.
+func fillNaNDense(m *mat64.Dense) {
+	r, c := m.Dims()
+	for i := 0; i < r; i++ {
+		for j := 0; j < c; j++ {
+			m.Set(i, j, math.NaN())
+		}
+	}
+}