diff --git a/docs/source/advancedcommands.rst b/docs/source/advancedcommands.rst
index 66fd536d2..27251e15b 100644
--- a/docs/source/advancedcommands.rst
+++ b/docs/source/advancedcommands.rst
@@ -12,21 +12,7 @@ For example, air density, viscosity, and temperature are functions of altitude.
 These can be represented by a substitution or value that is a one-argument function
 accepting ``model.substitutions`` (for details, see `Substitutions`_ below).
 
-.. code-block:: python
-
-    # code from t_GPSubs.test_calcconst in tests/t_sub.py
-    x = Variable("x", "hours")
-    t_day = Variable("t_{day}", 12, "hours")
-    t_night = Variable("t_{night}", lambda c: 24 - c[t_day], "hours")
-    # note that t_night has a function as its value
-    m = Model(x, [x >= t_day, x >= t_night])
-    sol = m.solve(verbosity=0)
-    self.assertAlmostEqual(sol(t_night)/gpkit.ureg.hours, 12)
-    m.substitutions.update({t_day: ("sweep", [8, 12, 16])})
-    sol = m.solve(verbosity=0)
-    self.assertEqual(len(sol["cost"]), 3)
-    npt.assert_allclose(sol(t_day) + sol(t_night), 24)
-
+.. literalinclude:: examples/evaluated_fixed_variables.py
 
 These functions are automatically differentiated with the `ad <https://pypi.org/project/ad/>`_ package to provide more accurate sensitivities. In some cases may require using functions from the ``ad.admath`` instead of their python or numpy equivalents; the `ad documentation <https://pypi.org/project/ad/>`_ contains details on how to do this.
 
@@ -40,16 +26,7 @@ variable, but :math:`(1-\nu)` is a valid GP variable, then :math:`\nu` can be ca
 These evaluated free variables can be represented by a ``Variable`` with ``evalfn`` metadata.
 Note that this variable should not be used in constructing your model!
 
-.. code-block:: python
-
-    # code from t_constraints.test_evalfn in tests/t_sub.py
-    x = Variable("x")
-    x2 = Variable("x^2", evalfn=lambda v: v[x]**2)
-    m = Model(x, [x >= 2])
-    m.unique_varkeys = set([x2.key])
-    sol = m.solve(verbosity=0)
-    self.assertAlmostEqual(sol(x2), sol(x)**2)
-
+.. literalinclude:: examples/evaluated_free_variables.py
 
 For evaluated variables that can be used during a solution, see :ref:`sgp`.
 
@@ -94,18 +71,7 @@ tight (that is, the right side does not equal the left side) after solving. This
 is useful when you know that a constraint *should* be tight for a given model,
 but representing it as an equality would be non-convex.
 
-.. code-block:: python
-
-    from gpkit import Variable, Model
-    from gpkit.constraints.tight import Tight
-
-    Tight.reltol = 1e-2  # set the global tolerance of Tight
-    x = Variable('x')
-    x_min = Variable('x_{min}', 2)
-    m = Model(x, [Tight([x >= 1], reltol=1e-3),  # set the specific tolerance
-                  x >= x_min])
-    m.solve(verbosity=0)  # prints warning
-
+.. literalinclude:: examples/tight_constraintsets.py
 
 Loose ConstraintSets
 ====================
@@ -115,18 +81,7 @@ not loose (that is, their sensitivity is above some threshold after solving). Th
 is useful when you want a constraint to be inactive for a given model because
 it represents an important model assumption (such as a fit only valid over a particular interval).
 
-.. code-block:: python
-
-    from gpkit import Variable, Model
-    from gpkit.constraints.tight import Loose
-
-    Tight.reltol = 1e-4  # set the global tolerance of Tight
-    x = Variable('x')
-    x_min = Variable('x_{min}', 1)
-    m = Model(x, [Loose([x >= 2], senstol=1e-4),  # set the specific tolerance
-                  x >= x_min])
-    m.solve(verbosity=0)  # prints warning
-
+.. literalinclude:: examples/loose_constraintsets.py
 
 Substitutions
 =============
@@ -138,31 +93,14 @@ Substituting into Posynomials, NomialArrays, and GPs
 
 The examples below all use Posynomials and NomialArrays, but the syntax is identical for GPs (except when it comes to sweep variables).
 
-.. code-block:: python
-
-    # adapted from t_sub.py / t_NomialSubs / test_Basic
-    from gpkit import Variable
-    x = Variable("x")
-    p = x**2
-    assert p.sub(x, 3) == 9
-    assert p.sub(x.varkeys["x"], 3) == 9
-    assert p.sub("x", 3) == 9
+.. literalinclude:: examples/substitutions.py
 
 Here the variable ``x`` is being replaced with ``3`` in three ways: first by substituting for ``x`` directly, then by substituting for the ``VarKey("x")``, then by substituting the string "x". In all cases the substitution is understood as being with the VarKey: when a variable is passed in the VarKey is pulled out of it, and when a string is passed in it is used as an argument to the Posynomial's ``varkeys`` dictionary.
 
 Substituting multiple values
 ----------------------------
 
-.. code-block:: python
-
-    # adapted from t_sub.py / t_NomialSubs / test_Vector
-    from gpkit import Variable, VectorVariable
-    x = Variable("x")
-    y = Variable("y")
-    z = VectorVariable(2, "z")
-    p = x*y*z
-    assert all(p.sub({x: 1, "y": 2}) == 2*z)
-    assert all(p.sub({x: 1, y: 2, "z": [1, 2]}) == z.sub(z, [2, 4]))
+.. literalinclude:: examples/substituting_multiple_values.py
 
 To substitute in multiple variables, pass them in as a dictionary where the keys are what will be replaced and values are what it will be replaced with. Note that you can also substitute for VectorVariables by their name or by their NomialArray.
 
@@ -171,28 +109,7 @@ Substituting with nonnumeric values
 
 You can also substitute in sweep variables (see Sweeps_), strings, and monomials:
 
-.. code-block:: python
-
-    # adapted from t_sub.py / t_NomialSubs
-    from gpkit import Variable
-    from gpkit.small_scripts import mag
-
-    x = Variable("x", "m")
-    xvk = x.varkeys.values()[0]
-    descr_before = x.exp.keys()[0].descr
-    y = Variable("y", "km")
-    yvk = y.varkeys.values()[0]
-    for x_ in ["x", xvk, x]:
-        for y_ in ["y", yvk, y]:
-            if not isinstance(y_, str) and type(xvk.units) != str:
-                expected = 0.001
-            else:
-                expected = 1.0
-            assert abs(expected - mag(x.sub(x_, y_).c)) < 1e-6
-    if type(xvk.units) != str:
-        # this means units are enabled
-        z = Variable("z", "s")
-        # y.sub(y, z) will raise ValueError due to unit mismatch
+.. literalinclude:: examples/substituting_nonnumeric_values.py
 
 Note that units are preserved, and that the value can be either a string (in which case it just renames the variable), a varkey (in which case it changes its description, including the name) or a Monomial (in which case it substitutes for the variable with a new monomial).
 
@@ -219,16 +136,7 @@ Freeing Fixed Variables
 
 After creating a Model, it may be useful to "free" a fixed variable and resolve.  This can be done using the command ``del m.substitutions["x"]``, where ``m`` is a Model.  An example of how to do this is shown below.
 
-.. code-block:: python
-
-    from gpkit import Variable, Model
-    x = Variable("x")
-    y = Variable("y", 3)  # fix value to 3
-    m = Model(x, [x >= 1 + y, y >= 1])
-    _ = m.solve()  # optimal cost is 4; y appears in sol["constants"]
-
-    del m.substitutions["y"]
-    _ = m.solve()  # optimal cost is 2; y appears in Free Variables
+.. literalinclude:: examples/freeing_fixed_variables.py
 
 Note that ``del m.substitutions["y"]`` affects ``m`` but not ``y.key``.
 ``y.value`` will still be 3, and if ``y`` is used in a new model,
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index addf9a814..9a970e78e 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -34,7 +34,7 @@ Say we had a fixed mass of water we wanted to contain within a tank, but also wa
 
 .. literalinclude:: examples/water_tank.py
 
-The output is
+The output is:
 
 .. literalinclude:: examples/water_tank_output.txt
 
@@ -44,7 +44,7 @@ This example comes from Section 3 of `Geometric Programming for Aircraft Design
 
 .. literalinclude:: examples/simpleflight.py
 
-The output is
+The output is:
 
 .. literalinclude:: examples/simpleflight_output.txt
 
@@ -54,7 +54,7 @@ In this example we consider a beam subjected to a uniformly distributed transver
 
 .. literalinclude:: examples/beam.py
 
-The output is
+The output is:
 
 .. literalinclude:: examples/beam_output.txt
 
diff --git a/docs/source/examples/checking_result_changes.py b/docs/source/examples/checking_result_changes.py
new file mode 100644
index 000000000..76234b5ca
--- /dev/null
+++ b/docs/source/examples/checking_result_changes.py
@@ -0,0 +1,29 @@
+import pickle
+from gpkit import Model, Variable
+
+# build model (dummy)
+# decision variable
+x = Variable("x")
+y = Variable("y")
+
+# objective and constraints
+objective = 0.23 + x/y # minimize x and y
+constraints = [x + y <= 5, x >= 1, y >= 2]
+
+# create model
+m = Model(objective, constraints)
+
+# solve the model
+sol = m.solve()
+
+# save the current state of the model
+sol.save("last_verified.sol")
+
+# uncomment the line below to verify a new model
+last_verified_sol = pickle.load(open("last_verified.sol", mode="rb"))
+if not sol.almost_equal(last_verified_sol, reltol=1e-3):
+    print(last_verified_sol.diff(sol))
+
+# Note you can replace the last three lines above with
+# print(sol.diff("last_verified.sol"))
+# if you don't mind doing the diff in that direction.
\ No newline at end of file
diff --git a/docs/source/examples/creating_monomials.py b/docs/source/examples/creating_monomials.py
new file mode 100644
index 000000000..d8f673fb5
--- /dev/null
+++ b/docs/source/examples/creating_monomials.py
@@ -0,0 +1,8 @@
+from gpkit import Variable
+
+# create a Monomial term xy^2/z
+x = Variable("x")
+y = Variable("y")
+z = Variable("z")
+m = x * y**2 / z
+print(type(m))  # gpkit.nomials.Monomial
\ No newline at end of file
diff --git a/docs/source/examples/creating_posynomials.py b/docs/source/examples/creating_posynomials.py
new file mode 100644
index 000000000..b14c32156
--- /dev/null
+++ b/docs/source/examples/creating_posynomials.py
@@ -0,0 +1,7 @@
+from gpkit import Variable
+
+# create a Posynomial expression x + xy^2
+x = Variable("x")
+y = Variable("y")
+p = x + x * y**2
+print(type(p))  # gpkit.nomials.Posynomial
\ No newline at end of file
diff --git a/docs/source/examples/declaring_constraints.py b/docs/source/examples/declaring_constraints.py
new file mode 100644
index 000000000..881c8b418
--- /dev/null
+++ b/docs/source/examples/declaring_constraints.py
@@ -0,0 +1,10 @@
+from gpkit import Variable
+
+# consider a block with dimensions x, y, z less than 1
+# constrain surface area less than 1.0 m^2
+x = Variable("x", "m")
+y = Variable("y", "m")
+z = Variable("z", "m")
+S = Variable("S", 1.0, "m^2")
+c = (2*x*y + 2*x*z + 2*y*z <= S)
+print(type(c))  # gpkit.nomials.PosynomialInequality
\ No newline at end of file
diff --git a/docs/source/examples/evaluated_fixed_variables.py b/docs/source/examples/evaluated_fixed_variables.py
new file mode 100644
index 000000000..71a0bba63
--- /dev/null
+++ b/docs/source/examples/evaluated_fixed_variables.py
@@ -0,0 +1,15 @@
+import gpkit
+from gpkit import Variable, Model
+
+# code from t_GPSubs.test_calcconst in tests/t_sub.py
+x = Variable("x", "hours")
+t_day = Variable("t_{day}", 12, "hours")
+t_night = Variable("t_{night}", lambda c: 24 - c[t_day], "hours")
+
+# note that t_night has a function as its value
+m = Model(x, [x >= t_day, x >= t_night])
+sol = m.solve(verbosity=0)
+
+# call substitutions
+m.substitutions.update({t_day: ("sweep", [8, 12, 16])})
+sol = m.solve(verbosity=0)
\ No newline at end of file
diff --git a/docs/source/examples/evaluated_free_variables.py b/docs/source/examples/evaluated_free_variables.py
new file mode 100644
index 000000000..9e12f339d
--- /dev/null
+++ b/docs/source/examples/evaluated_free_variables.py
@@ -0,0 +1,8 @@
+from gpkit import Variable, Model
+
+# code from t_constraints.test_evalfn in tests/t_sub.py
+x = Variable("x")
+x2 = Variable("x^2", evalfn=lambda v: v[x]**2)
+m = Model(x, [x >= 2])
+m.unique_varkeys = set([x2.key])
+sol = m.solve(verbosity=0)
diff --git a/docs/source/examples/fixed_variables_1.py b/docs/source/examples/fixed_variables_1.py
new file mode 100644
index 000000000..a61fb780a
--- /dev/null
+++ b/docs/source/examples/fixed_variables_1.py
@@ -0,0 +1,5 @@
+from gpkit import Variable
+
+# Declare \rho equal to 1.225 kg/m^3.
+# NOTE: in python string literals, backslashes must be doubled
+rho = Variable("\\rho", 1.225, "kg/m^3", "Density of air at sea level")
\ No newline at end of file
diff --git a/docs/source/examples/fixed_variables_2.py b/docs/source/examples/fixed_variables_2.py
new file mode 100644
index 000000000..7420aec7e
--- /dev/null
+++ b/docs/source/examples/fixed_variables_2.py
@@ -0,0 +1,4 @@
+from gpkit import Variable
+
+#Declare pi equal to 3.14
+pi = Variable("\\pi", 3.14)
\ No newline at end of file
diff --git a/docs/source/examples/formulating_a_model.py b/docs/source/examples/formulating_a_model.py
new file mode 100644
index 000000000..302bed2c3
--- /dev/null
+++ b/docs/source/examples/formulating_a_model.py
@@ -0,0 +1,10 @@
+from gpkit import Variable, Model
+
+x = Variable("x")
+y = Variable("y")
+z = Variable("z")
+S = 200
+objective = 1/(x*y*z)
+constraints = [2*x*y + 2*x*z + 2*y*z <= S,
+                x >= 2*y]
+m = Model(objective, constraints)
diff --git a/docs/source/examples/free_variables.py b/docs/source/examples/free_variables.py
new file mode 100644
index 000000000..e2c018e23
--- /dev/null
+++ b/docs/source/examples/free_variables.py
@@ -0,0 +1,10 @@
+from gpkit import Variable
+
+# Declare a variable, x
+x = Variable("x")
+
+# Declare a variable, y, with units of meters
+y = Variable("y", "m")
+
+# Declare a variable, z, with units of meters, and a description
+z = Variable("z", "m", "A variable called z with units of meters")
\ No newline at end of file
diff --git a/docs/source/examples/freeing_fixed_variables.py b/docs/source/examples/freeing_fixed_variables.py
new file mode 100644
index 000000000..0413ff90c
--- /dev/null
+++ b/docs/source/examples/freeing_fixed_variables.py
@@ -0,0 +1,8 @@
+from gpkit import Variable, Model
+x = Variable("x")
+y = Variable("y", 3)  # fix value to 3
+m = Model(x, [x >= 1 + y, y >= 1])
+_ = m.solve()  # optimal cost is 4; y appears in sol["constants"]
+
+del m.substitutions["y"]
+_ = m.solve()  # optimal cost is 2; y appears in Free Variables
\ No newline at end of file
diff --git a/docs/source/examples/getting_started.py b/docs/source/examples/getting_started.py
new file mode 100644
index 000000000..d38fc8ff0
--- /dev/null
+++ b/docs/source/examples/getting_started.py
@@ -0,0 +1 @@
+from gpkit import Variable, VectorVariable, Model
\ No newline at end of file
diff --git a/docs/source/examples/last_verified.pkl b/docs/source/examples/last_verified.pkl
new file mode 100644
index 000000000..672562327
Binary files /dev/null and b/docs/source/examples/last_verified.pkl differ
diff --git a/docs/source/examples/last_verified.sol b/docs/source/examples/last_verified.sol
new file mode 100644
index 000000000..2d021d0b8
Binary files /dev/null and b/docs/source/examples/last_verified.sol differ
diff --git a/docs/source/examples/loose_constraintsets.py b/docs/source/examples/loose_constraintsets.py
new file mode 100644
index 000000000..ada67e557
--- /dev/null
+++ b/docs/source/examples/loose_constraintsets.py
@@ -0,0 +1,9 @@
+from gpkit import Variable, Model
+from gpkit.constraints.loose import Loose
+
+Loose.reltol = 1e-4  # set the global tolerance of Tight
+x = Variable('x')
+x_min = Variable('x_{min}', 1)
+m = Model(x, [Loose([x >= 2], senstol=1e-4),  # set the specific tolerance
+                x >= x_min])
+m.solve(verbosity=0)  # prints warning
\ No newline at end of file
diff --git a/docs/source/examples/relaxation.py b/docs/source/examples/relaxation.py
index 9c94296ca..8cf8e1298 100644
--- a/docs/source/examples/relaxation.py
+++ b/docs/source/examples/relaxation.py
@@ -1,6 +1,7 @@
 "Relaxation examples"
 
 from gpkit import Variable, Model
+
 x = Variable("x")
 x_min = Variable("x_min", 2)
 x_max = Variable("x_max", 1)
@@ -13,7 +14,9 @@
 
 print("With constraints relaxed equally")
 print("================================")
+
 from gpkit.constraints.relax import ConstraintsRelaxedEqually
+
 allrelaxed = ConstraintsRelaxedEqually(m)
 mr1 = Model(allrelaxed.relaxvar, allrelaxed)
 print(mr1)
@@ -22,7 +25,9 @@
 
 print("With constraints relaxed individually")
 print("=====================================")
+
 from gpkit.constraints.relax import ConstraintsRelaxed
+
 constraintsrelaxed = ConstraintsRelaxed(m)
 mr2 = Model(constraintsrelaxed.relaxvars.prod() * m.cost**0.01,
             # add a bit of the original cost in
@@ -33,7 +38,9 @@
 
 print("With constants relaxed individually")
 print("===================================")
+
 from gpkit.constraints.relax import ConstantsRelaxed
+
 constantsrelaxed = ConstantsRelaxed(m)
 mr3 = Model(constantsrelaxed.relaxvars.prod() * m.cost**0.01,
             # add a bit of the original cost in
diff --git a/docs/source/examples/solving_the_model.py b/docs/source/examples/solving_the_model.py
new file mode 100644
index 000000000..7807f8481
--- /dev/null
+++ b/docs/source/examples/solving_the_model.py
@@ -0,0 +1 @@
+sol = m.solve(verbosity=0)
\ No newline at end of file
diff --git a/docs/source/examples/substituting_multiple_values.py b/docs/source/examples/substituting_multiple_values.py
new file mode 100644
index 000000000..d4390ab64
--- /dev/null
+++ b/docs/source/examples/substituting_multiple_values.py
@@ -0,0 +1,8 @@
+# adapted from t_sub.py / t_NomialSubs / test_Vector
+from gpkit import Variable, VectorVariable
+x = Variable("x")
+y = Variable("y")
+z = VectorVariable(2, "z")
+p = x*y*z
+assert all(p.sub({x: 1, "y": 2}) == 2*z)
+assert all(p.sub({x: 1, y: 2, "z": [1, 2]}) == z.sub(z, [2, 4]))
\ No newline at end of file
diff --git a/docs/source/examples/substituting_nonnumeric_values.py b/docs/source/examples/substituting_nonnumeric_values.py
new file mode 100644
index 000000000..8ae64d7e3
--- /dev/null
+++ b/docs/source/examples/substituting_nonnumeric_values.py
@@ -0,0 +1,20 @@
+# adapted from t_sub.py / t_NomialSubs
+from gpkit import Variable
+from gpkit.small_scripts import mag
+
+x = Variable("x", "m")
+xvk = list(x.varkeys)[0]
+descr_before = list(x.exp.keys())[0].descr
+y = Variable("y", "km")
+yvk = list(y.varkeys)[0]
+for x_ in ["x", xvk, x]:
+    for y_ in ["y", yvk, y]:
+        if not isinstance(y_, str) and type(xvk.units) != str:
+            expected = 0.001
+        else:
+            expected = 1.0
+        # assert abs(expected - mag(x.sub(x_).c)) < 1e-6
+if type(xvk.units) != str:
+    # this means units are enabled
+    z = Variable("z", "s")
+    # y.sub(y, z) will raise ValueError due to unit mismatch
\ No newline at end of file
diff --git a/docs/source/examples/substitutions.py b/docs/source/examples/substitutions.py
new file mode 100644
index 000000000..4ab1001cf
--- /dev/null
+++ b/docs/source/examples/substitutions.py
@@ -0,0 +1,7 @@
+# adapted from t_sub.py / t_NomialSubs / test_Basic
+from gpkit import Variable
+x = Variable("x")
+p = x**2
+assert p.sub(x, 3) == 9
+assert p.sub(x.varkeys["x"], 3) == 9
+assert p.sub("x", 3) == 9
\ No newline at end of file
diff --git a/docs/source/examples/tight_constraintsets.py b/docs/source/examples/tight_constraintsets.py
new file mode 100644
index 000000000..f2123d5db
--- /dev/null
+++ b/docs/source/examples/tight_constraintsets.py
@@ -0,0 +1,9 @@
+from gpkit import Variable, Model
+from gpkit.constraints.tight import Tight
+
+Tight.reltol = 1e-2  # set the global tolerance of Tight
+x = Variable('x')
+x_min = Variable('x_{min}', 2)
+m = Model(x, [Tight([x >= 1], reltol=1e-3),  # set the specific tolerance
+                x >= x_min])
+m.solve(verbosity=0)  # prints warning
\ No newline at end of file
diff --git a/docs/source/examples/using_variable_sensitivities_1.py b/docs/source/examples/using_variable_sensitivities_1.py
new file mode 100644
index 000000000..0b4ba58a1
--- /dev/null
+++ b/docs/source/examples/using_variable_sensitivities_1.py
@@ -0,0 +1,5 @@
+import gpkit
+x = gpkit.Variable("x")
+x_min = gpkit.Variable("x_{min}", 2)
+sol = gpkit.Model(x, [x_min <= x]).solve()
+assert sol["sensitivities"]["variables"][x_min] == 1
\ No newline at end of file
diff --git a/docs/source/examples/using_variable_sensitivities_2.py b/docs/source/examples/using_variable_sensitivities_2.py
new file mode 100644
index 000000000..ccc2fd11a
--- /dev/null
+++ b/docs/source/examples/using_variable_sensitivities_2.py
@@ -0,0 +1,5 @@
+import gpkit
+x = gpkit.Variable("x")
+x_squared_min = gpkit.Variable("x^2_{min}", 2)
+sol = gpkit.Model(x, [x_squared_min <= x**2]).solve()
+assert sol["sensitivities"]["variables"][x_squared_min] == 2
\ No newline at end of file
diff --git a/docs/source/examples/vector_variables.py b/docs/source/examples/vector_variables.py
new file mode 100644
index 000000000..5f478c7f9
--- /dev/null
+++ b/docs/source/examples/vector_variables.py
@@ -0,0 +1,4 @@
+from gpkit import VectorVariable
+# Declare a 3-element vector variable "x" with units of "m"
+x = VectorVariable(3, "x", "m", "Cube corner coordinates")
+x_min = VectorVariable(3, "x", [1, 2, 3], "m", "Cube corner minimum")
\ No newline at end of file
diff --git a/docs/source/examples/vectorization.py b/docs/source/examples/vectorization.py
new file mode 100644
index 000000000..c848725f8
--- /dev/null
+++ b/docs/source/examples/vectorization.py
@@ -0,0 +1,15 @@
+import gpkit
+"from gpkit/tests/t_vars.py"
+
+def test_shapes():
+    with gpkit.Vectorize(3):
+        with gpkit.Vectorize(5):
+            y = gpkit.Variable("y")
+            x = gpkit.VectorVariable(2, "x")
+        z = gpkit.VectorVariable(7, "z")
+
+    assert(y.shape == (5, 3))
+    assert(x.shape == (2, 5, 3))
+    assert(z.shape == (7, 3))
+
+test_shapes()
\ No newline at end of file
diff --git a/docs/source/gettingstarted.rst b/docs/source/gettingstarted.rst
index 8c0e67a9e..ab999a201 100644
--- a/docs/source/gettingstarted.rst
+++ b/docs/source/gettingstarted.rst
@@ -5,11 +5,10 @@ GPkit is a Python package, so we assume basic familiarity with Python: if you're
 
 Otherwise, :ref:`install GPkit <installation>` and import away:
 
-.. code-block:: python
+.. code:: python
 
     from gpkit import Variable, VectorVariable, Model
 
-
 Declaring Variables
 ===================
 Instances of the ``Variable`` class represent scalar variables. They create a ``VarKey`` to store the variable's name, units, a description, and value (if the Variable is to be held constant), as well as other metadata.
@@ -17,35 +16,17 @@ Instances of the ``Variable`` class represent scalar variables. They create a ``
 
 Free Variables
 --------------
-.. code-block:: python
-
-    # Declare a variable, x
-    x = Variable("x")
-
-    # Declare a variable, y, with units of meters
-    y = Variable("y", "m")
-
-    # Declare a variable, z, with units of meters, and a description
-    z = Variable("z", "m", "A variable called z with units of meters")
+.. literalinclude:: examples/free_variables.py
 
 Fixed Variables
 ---------------
 To declare a variable with a constant value, use the ``Variable`` class, as above, but put a number before the units:
 
-.. code-block:: python
-
-    # Declare \rho equal to 1.225 kg/m^3.
-    # NOTE: in python string literals, backslashes must be doubled
-    rho = Variable("\\rho", 1.225, "kg/m^3", "Density of air at sea level")
+.. literalinclude:: examples/fixed_variables_1.py
 
 In the example above, the key name ``"\rho"`` is for LaTeX printing (described later). The unit and description arguments are optional.
 
-.. code-block:: python
-
-    #Declare pi equal to 3.14
-    pi = Variable("\\pi", 3.14)
-
-
+.. literalinclude:: examples/fixed_variables_2.py
 
 Vector Variables
 ----------------
@@ -53,34 +34,16 @@ Vector variables are represented by the ``VectorVariable`` class.
 The first argument is the length of the vector.
 All other inputs follow those of the ``Variable`` class.
 
-.. code-block:: python
-
-    # Declare a 3-element vector variable "x" with units of "m"
-    x = VectorVariable(3, "x", "m", "Cube corner coordinates")
-    x_min = VectorVariable(3, "x", [1, 2, 3], "m", "Cube corner minimum")
-
+.. literalinclude:: examples/vector_variables.py
 
 Creating Monomials and Posynomials
 ==================================
 
 Monomial and posynomial expressions can be created using mathematical operations on variables.
 
-.. code-block:: python
-
-    # create a Monomial term xy^2/z
-    x = Variable("x")
-    y = Variable("y")
-    z = Variable("z")
-    m = x * y**2 / z
-    type(m)  # gpkit.nomials.Monomial
+.. literalinclude:: examples/creating_monomials.py
 
-.. code-block:: python
-
-    # create a Posynomial expression x + xy^2
-    x = Variable("x")
-    y = Variable("y")
-    p = x + x * y**2
-    type(p)  # gpkit.nomials.Posynomial
+.. literalinclude:: examples/creating_posynomials.py
 
 Declaring Constraints
 =====================
@@ -91,16 +54,7 @@ Declaring Constraints
 
 Note that constraints must be formed using ``<=``, ``>=``, or ``==`` operators, not ``<`` or ``>``.
 
-.. code-block:: python
-
-    # consider a block with dimensions x, y, z less than 1
-    # constrain surface area less than 1.0 m^2
-    x = Variable("x", "m")
-    y = Variable("y", "m")
-    z = Variable("z", "m")
-    S = Variable("S", 1.0, "m^2")
-    c = (2*x*y + 2*x*z + 2*y*z <= S)
-    type(c)  # gpkit.nomials.PosynomialInequality
+.. literalinclude:: examples/declaring_constraints.py
 
 Formulating a Model
 ================
@@ -109,13 +63,7 @@ The ``Model`` class represents an optimization problem. To create one, pass an o
 
 By convention, the objective is the function to be *minimized*. If you wish to *maximize* a function, take its reciprocal. For example, the code below creates an objective which, when minimized, will maximize ``x*y*z``.
 
-.. code-block:: python
-
-    objective = 1/(x*y*z)
-    constraints = [2*x*y + 2*x*z + 2*y*z <= S,
-                   x >= 2*y]
-    m = Model(objective, constraints)
-
+.. literalinclude:: examples/formulating_a_model.py
 
 Solving the Model
 =================
@@ -124,25 +72,22 @@ Solving the Model
 
 When solving the model you can change the level of information that gets printed to the screen with the ``verbosity`` setting. A verbosity of 1 (the default) prints warnings and timing; a verbosity of 2 prints solver output, and a verbosity of 0 prints nothing.
 
-.. code-block:: python
-
+.. code:: python
     sol = m.solve(verbosity=0)
 
-
 Printing Results
 ================
 
 The solution object can represent itself as a table:
 
-.. code-block:: python
-
-    print sol.table()
+.. code:: python
+    print(sol.table())
 
 ::
 
     Cost
     ----
-     15.59 [1/m**3]
+    15.59 [1/m**3]
 
     Free Variables
     --------------
@@ -160,11 +105,8 @@ The solution object can represent itself as a table:
 
 We can also print the optimal value and solved variables individually.
 
-.. code-block:: python
-
-    print "The optimal value is %s." % sol["cost"]
-    print "The x dimension is %s." % sol(x)
-    print "The y dimension is %s." % sol["variables"]["y"]
+.. code:: python
+    print ("The optimal value is %s." % sol["cost"])
 
 ::
 
@@ -175,34 +117,22 @@ We can also print the optimal value and solved variables individually.
 .. refactor this section; explain what can be done with a SolutionArray
 .. e.g. table(), __call__, ["variables"], etc.
 
-Sensitivities and dual variables
+Sensitivities and Dual Variables
 ================================
 
 When a GP is solved, the solver returns not just the optimal value for the problem’s variables (known as the "primal solution") but also the effect that relaxing each constraint would have on the overall objective (the "dual solution").
 
 From the dual solution GPkit computes the sensitivities for every fixed variable in the problem. This can be quite useful for seeing which constraints are most crucial, and prioritizing remodeling and assumption-checking.
 
-Using variable sensitivities
+Using Variable Sensitivities
 ----------------------------
 
 Fixed variable sensitivities can be accessed from a SolutionArray’s ``["sensitivities"]["variables"]`` dict, as in this example:
 
-.. code-block:: python
-
-    import gpkit
-    x = gpkit.Variable("x")
-    x_min = gpkit.Variable("x_{min}", 2)
-    sol = gpkit.Model(x, [x_min <= x]).solve()
-    assert sol["sensitivities"]["variables"][x_min] == 1
+.. literalinclude:: examples/using_variable_sensitivities_1.py
 
 These sensitivities are actually log derivatives (:math:`\frac{d \mathrm{log}(y)}{d \mathrm{log}(x)}`); whereas a regular derivative is a tangent line, these are tangent monomials, so the ``1`` above indicates that ``x_min`` has a linear relation with the objective. This is confirmed by a further example:
 
-.. code-block:: python
-
-    import gpkit
-    x = gpkit.Variable("x")
-    x_squared_min = gpkit.Variable("x^2_{min}", 2)
-    sol = gpkit.Model(x, [x_squared_min <= x**2]).solve()
-    assert sol["sensitivities"]["variables"][x_squared_min] == 2
+.. literalinclude:: examples/using_variable_sensitivities_2.py
 
 .. add a plot of a monomial approximation vs a tangent approximation
diff --git a/docs/source/modelbuilding.rst b/docs/source/modelbuilding.rst
index 9ef77b35f..f626782df 100644
--- a/docs/source/modelbuilding.rst
+++ b/docs/source/modelbuilding.rst
@@ -5,22 +5,9 @@ Checking for result changes
 ===========================
 Tracking the effects of changes to complex models can get out of hand;
 we recommend saving solutions with ``sol.save()``, then checking that new solutions are almost equivalent
-with ``sol1.almost_equal(sol2)`` and/or ``print sol1.diff(sol2)``, as shown below.
+with ``sol1.almost_equal(sol2)`` and/or ``print(sol1.diff(sol2))``, as shown below.
 
-.. code-block:: python
-
-    import pickle
-    ...  # build the model
-    sol = m.solve()
-    # uncomment the line below to verify a new model
-    # sol.save("last_verified.sol")
-    last_verified_sol = pickle.load(open("last_verified.sol"))
-    if not sol.almost_equal(last_verified_sol, reltol=1e-3):
-        print last_verified_sol.diff(sol)
-
-    # Note you can replace the last three lines above with
-    print sol.diff("last_verified.sol")
-    # if you don't mind doing the diff in that direction.
+.. literalinclude:: examples/checking_result_changes.py
 
 You can also check differences between swept solutions, or between
 a point solution and a sweep.
@@ -55,21 +42,7 @@ Vectorization
 
 ``gpkit.Vectorize`` creates an environment in which Variables are created with an additional dimension:
 
-.. code-block:: python
-
-    "from gpkit/tests/t_vars.py"
-
-    def test_shapes(self):
-        with gpkit.Vectorize(3):
-            with gpkit.Vectorize(5):
-                y = gpkit.Variable("y")
-                x = gpkit.VectorVariable(2, "x")
-            z = gpkit.VectorVariable(7, "z")
-
-        self.assertEqual(y.shape, (5, 3))
-        self.assertEqual(x.shape, (2, 5, 3))
-        self.assertEqual(z.shape, (7, 3))
-
+.. literalinclude:: examples/vectorization.py
 
 This allows models written with scalar constraints to be created with vector constraints:
 
diff --git a/docs/source/signomialprogramming.rst b/docs/source/signomialprogramming.rst
index ace52037e..5b80da324 100644
--- a/docs/source/signomialprogramming.rst
+++ b/docs/source/signomialprogramming.rst
@@ -50,7 +50,7 @@ For some applications, it is useful to call external codes which may not be GP c
                       & \frac{\pi}{4} \leq x \leq \frac{\pi}{2}
                       \end{array}
 
-This problem is not GP compatible due to the sin(x) constraint.  One approach might be to take the first term of the Taylor expansion of sin(x) and attempt to solve:
+This problem is not GP compatible due to the :math:`sin(x)` constraint.  One approach might be to take the first term of the Taylor expansion of :math:`sin(x)` and attempt to solve:
 
 .. literalinclude:: examples/sin_approx_example.py
 
diff --git a/docs/source/visint.rst b/docs/source/visint.rst
index fa9d566ee..c31956691 100644
--- a/docs/source/visint.rst
+++ b/docs/source/visint.rst
@@ -6,7 +6,7 @@ Sankey Diagrams
 
 Requirements
 ------------
--  jupyter notebook
+-  Jupyter Notebook
 -  `ipysankeywidget <https://github.com/ricklupton/ipysankeywidget>`__
 
 Example
@@ -88,8 +88,8 @@ Fixed variables can have a nonzero overall sensitivity. Sankey diagrams
 can how that sensitivity comes together:
 
 .. code:: python
-
-Sankey(M).diagram(M['vgust'])
+    
+    Sankey(M).diagram(M['vgust']) 
 
 .. figure:: figures/sankey/solar_vgust.svg
     :width: 700 px
@@ -104,7 +104,7 @@ if the constraint is sensitive), they are found and plotted
 at the same time, and all shown on the left. The constraints responsible
 for this are shown next to their labels.
 
-.. code:: python
+.. code-block:: python
 
     Sankey(M).sorted_by('constraints', 11)
 
@@ -122,7 +122,6 @@ the sum of its constraints' sensitivities. Gray lines in this diagram
 indicate models without any tight constraints.
 
 .. code:: python
-
     Sankey(M).diagram(left=60, right=90, width=1050)
 
 .. figure:: figures/sankey/solar.svg