diff --git a/doc/source/algorithms/critical_curves.ipynb b/doc/source/algorithms/critical_curves.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..c90ede679cad8a43259e0e535e94a3874ff33471
--- /dev/null
+++ b/doc/source/algorithms/critical_curves.ipynb
@@ -0,0 +1,38 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "8218498b",
+   "metadata": {},
+   "source": [
+    "# Critical curves\n",
+    "\n",
+    "A pure fluid has a single vapor-liquid critical point, but mixtures are different:\n",
+    "\n",
+    "* They may have multiple (or zero!) critical points for given mixture composition\n",
+    "* The critical curves may not emanate from the pure fluid endpoints"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/algorithms/index.rst b/doc/source/algorithms/index.rst
index 57169ed566c6e16605cea6f9d6dea353106dde89..13225b9ff81375465bb8140e44e6e13a58b07862 100644
--- a/doc/source/algorithms/index.rst
+++ b/doc/source/algorithms/index.rst
@@ -5,4 +5,10 @@ Algorithms
    :maxdepth: 2
    :caption: Contents:
 
-   VLE
\ No newline at end of file
+   VLE
+   critical_curves
+
+Information
+-----------
+
+The algorithms are written in a very generic way; they take an instance of a thermodynamic model, and the necessary derivatives are calculated from this model with automatic differentiation (or similar). In that way, implementing a model is all that is required to enable its use in the calculation of critical curves or to trace the phase equilibria.  Determining the starting values, on the other hand, may require model-specific assistance, for instance with superancillary equations.
\ No newline at end of file
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 5449eb132bd300c674dc79b743197507053894f2..debde9be0add829404924fae9526a1aedb3dee0c 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -4,6 +4,8 @@
 # list see the documentation:
 # https://www.sphinx-doc.org/en/master/usage/configuration.html
 
+import os, subprocess
+
 # -- Path setup --------------------------------------------------------------
 
 # If extensions (or modules to document with autodoc) are in another directory,
@@ -25,6 +27,12 @@ author = 'Ian Bell'
 import teqp
 release = teqp.__version__
 
+# -- Exeucute all notebooks --------------------------------------------------
+
+for path, dirs, files in os.walk('.'):
+    for file in files:
+        if file.endswith('.ipynb'):
+            subprocess.check_output(f'jupyter nbconvert --ExecutePreprocessor.allow_errors=True --to notebook --output {file} --execute {file}', cwd=path)
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/doc/source/derivs/derivs.ipynb b/doc/source/derivs/derivs.ipynb
index 7af50d86aa1f0ffc4e236a1302625a49d34749bb..aa39bc31243defd2b4cff0dd1a92e88f426da8fa 100644
--- a/doc/source/derivs/derivs.ipynb
+++ b/doc/source/derivs/derivs.ipynb
@@ -10,9 +10,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 1,
    "id": "ca35cc29",
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:39:50.124207Z",
+     "iopub.status.busy": "2022-07-06T18:39:50.124207Z",
+     "iopub.status.idle": "2022-07-06T18:39:50.261696Z",
+     "shell.execute_reply": "2022-07-06T18:39:50.261181Z"
+    }
+   },
    "outputs": [],
    "source": [
     "import teqp\n",
@@ -21,9 +28,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 2,
    "id": "8d254c40",
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:39:50.264736Z",
+     "iopub.status.busy": "2022-07-06T18:39:50.264736Z",
+     "iopub.status.idle": "2022-07-06T18:39:50.277541Z",
+     "shell.execute_reply": "2022-07-06T18:39:50.276780Z"
+    }
+   },
    "outputs": [],
    "source": [
     "Tc_K = [300]\n",
@@ -34,9 +48,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 3,
    "id": "49ee3453",
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:39:50.280537Z",
+     "iopub.status.busy": "2022-07-06T18:39:50.280537Z",
+     "iopub.status.idle": "2022-07-06T18:39:50.292555Z",
+     "shell.execute_reply": "2022-07-06T18:39:50.292555Z"
+    }
+   },
    "outputs": [
     {
      "data": {
@@ -44,7 +65,7 @@
        "-0.06836660379313926"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
diff --git a/doc/source/getting_started/index.rst b/doc/source/getting_started/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e573440f81c239c85fb946d2a63d3006ee49ac91
--- /dev/null
+++ b/doc/source/getting_started/index.rst
@@ -0,0 +1,12 @@
+Getting Started
+===============
+
+teqp is a C++-based library with wrappers. The documentation is based on the Python wrapper because it can be readily integrated with the documentation tools and can be auto-generated at documentation build time.
+
+The library can (for Python) be installed with:
+
+.. code::
+
+   pip install teqp
+
+as the binary wheels for all major platforms are provided on pypi.
\ No newline at end of file
diff --git a/doc/source/index.rst b/doc/source/index.rst
index d321eadc912705e22ba4ce5dee5673f3c10c71e5..9cbb358068e31cf5165b5c0b4c04738ff1d4506a 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -5,8 +5,9 @@ Welcome to teqp's documentation!
    :maxdepth: 2
    :caption: Contents:
 
-   derivs/index
+   getting_started/index
    models/index
+   derivs/index
    algorithms/index
 
 Indices and tables
diff --git a/doc/source/models/CPA.ipynb b/doc/source/models/CPA.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2a597f271e2528c8ae0722657d3172c60345f268
--- /dev/null
+++ b/doc/source/models/CPA.ipynb
@@ -0,0 +1,43 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "df447d16",
+   "metadata": {},
+   "source": [
+    "# Cubic EOS\n",
+    "\n",
+    "Something about cubic EOS"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a884dab",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/models/PCSAFT.ipynb b/doc/source/models/PCSAFT.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2a597f271e2528c8ae0722657d3172c60345f268
--- /dev/null
+++ b/doc/source/models/PCSAFT.ipynb
@@ -0,0 +1,43 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "df447d16",
+   "metadata": {},
+   "source": [
+    "# Cubic EOS\n",
+    "\n",
+    "Something about cubic EOS"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a884dab",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/models/cubics.ipynb b/doc/source/models/cubics.ipynb
index 2a597f271e2528c8ae0722657d3172c60345f268..1b1ba08af73bfbb8d92bbffbd7af4dca1954e601 100644
--- a/doc/source/models/cubics.ipynb
+++ b/doc/source/models/cubics.ipynb
@@ -5,18 +5,90 @@
    "id": "df447d16",
    "metadata": {},
    "source": [
-    "# Cubic EOS\n",
+    "# General cubics\n",
     "\n",
-    "Something about cubic EOS"
+    "The reduced residual Helmholtz energy for the main cubic EOS (van der Waals, Peng-Robinson, and Soave-Redlich-Kwong) can be written in a common form:\n",
+    "\n",
+    "$$ \\begin{equation}\n",
+    "\\label{eq:alphar_from_psi}\n",
+    "\\alpha^r = \\psi^{(-)} - \\dfrac{\\tau a_m}{RT_r } \\psi^{(+)}.\n",
+    "\\end{equation}$$\n",
+    "\n",
+    "$$ \\begin{eqnarray}\n",
+    "\\psi^{(-)} &=& \\int_0^\\delta\\dfrac{b_m\\rho_r }{1-b_m\\delta\\rho_r }{\\rm d}\\delta \\label{eq:psiminusintegral}\\\\\n",
+    "           &=&-\\ln(1-b_m\\rho ). \\label{eq:psiminusresult}\n",
+    "\\end{eqnarray} $$\n",
+    "\n",
+    "$$ \\begin{eqnarray}\n",
+    "\\psi^{(+)} &=& \\int_0^\\delta \\dfrac{\\rho_r}{\\left(1+\\Delta_1 b_m\\delta\\rho_r \\right)\\left(1+\\Delta_2 b_m\\delta\\rho_r \\right)} {\\rm d}\\delta \\label{eq:psiplusintegral}\\\\\n",
+    "           &=& \\dfrac{\\ln\\left(\\dfrac{\\Delta_1 b_m\\rho+1}{\\Delta_2b_m\\rho+1}\\right)}{b_m(\\Delta_1-\\Delta_2)}\\label{eq:psiplusresult}\n",
+    "\\end{eqnarray} $$\n",
+    "\n",
+    "with the constants given by:\n",
+    "\n",
+    "* vdW: $\\Delta_1=0$, $\\Delta_2=0$\n",
+    "* SRK: $\\Delta_1=1$, $\\Delta_2=0$\n",
+    "* PR: $\\Delta_1=1+\\sqrt{2}$, $\\Delta_2=1-\\sqrt{2}$\n",
+    "\n",
+    "The quantities $a_m$ and $b_m$ are described (with exact solutions for the numerical coefficients) for each of these EOS in https://pubs.acs.org/doi/abs/10.1021/acs.iecr.1c00847.\n",
+    "\n",
+    "The models in teqp are instantiated based on knowledge of the critical temperature, pressure, and acentric factor.  Thereafter all quantities are obtained from derivatives of $\\alpha^r$."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "id": "9a884dab",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "import teqp\n",
+    "\n",
+    "# Values taken from http://dx.doi.org/10.6028/jres.121.011\n",
+    "Tc_K = [ 190.564, 154.581, 150.687 ]\n",
+    "pc_Pa = [ 4599200, 5042800, 4863000 ]\n",
+    "acentric = [ 0.011, 0.022, -0.002 ]\n",
+    "\n",
+    "# Instantiate Peng-Robinson model\n",
+    "modelPR = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
+    "\n",
+    "# Instantiate Soave-Redlich-Kwong model\n",
+    "modelSRK = teqp.canonical_SRK(Tc_K, pc_Pa, acentric)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "615c7830",
+   "metadata": {},
+   "source": [
+    "## Adjusting k_ij\n",
+    "\n",
+    "Fine-tuned values of $k_{ij}$ can be provided when instantiating the model, for Peng-Robinson and SRK.  A complete matrix of all the $k_{ij}$ values must be provided. This allows for asymmetric mixing models in which $k_{ij}\\neq k_{ji}$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "a8c87890",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<teqp.teqp.GenericCubic at 0x20b222799b0>"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "k_12 = 0.01\n",
+    "kmat = [[0,k_12,0],[k_12,0,0],[0,0,0]]\n",
+    "teqp.canonical_PR(Tc_K, pc_Pa, acentric, kmat)\n",
+    "teqp.canonical_SRK(Tc_K, pc_Pa, acentric, kmat)"
+   ]
   }
  ],
  "metadata": {
diff --git a/doc/source/models/index.rst b/doc/source/models/index.rst
index 355677785073736c8116c9d9e04c2942ef887de9..30b3ec5c524c165dd24c8bc1e008cf3837232c3f 100644
--- a/doc/source/models/index.rst
+++ b/doc/source/models/index.rst
@@ -7,3 +7,5 @@ Models
 
    cubics
    multifluid
+   CPA
+   PCSAFT
diff --git a/doc/source/models/multifluid.ipynb b/doc/source/models/multifluid.ipynb
index 82d80bca097c09fe06d8ad6f52e7ee59b7db5428..76037228990d9c25f842ae2e0824bafc4595ef77 100644
--- a/doc/source/models/multifluid.ipynb
+++ b/doc/source/models/multifluid.ipynb
@@ -12,12 +12,19 @@
   {
    "cell_type": "code",
    "execution_count": 1,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:06.606089Z",
+     "iopub.status.busy": "2022-07-06T18:40:06.605465Z",
+     "iopub.status.idle": "2022-07-06T18:40:07.147629Z",
+     "shell.execute_reply": "2022-07-06T18:40:07.146669Z"
+    }
+   },
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "'0.8.1'"
+       "'0.9.2'"
       ]
      },
      "execution_count": 1,
@@ -43,13 +50,20 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:07.163785Z",
+     "iopub.status.busy": "2022-07-06T18:40:07.163785Z",
+     "iopub.status.idle": "2022-07-06T18:40:17.280030Z",
+     "shell.execute_reply": "2022-07-06T18:40:17.278333Z"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "13.4 ms ± 158 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+      "12.9 ms ± 4.32 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
      ]
     }
    ],
@@ -62,13 +76,20 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:17.286031Z",
+     "iopub.status.busy": "2022-07-06T18:40:17.285350Z",
+     "iopub.status.idle": "2022-07-06T18:40:27.300085Z",
+     "shell.execute_reply": "2022-07-06T18:40:27.299050Z"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "121 ms ± 3.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+      "120 ms ± 19.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
      ]
     }
    ],
@@ -89,13 +110,20 @@
   {
    "cell_type": "code",
    "execution_count": 4,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:27.307115Z",
+     "iopub.status.busy": "2022-07-06T18:40:27.307115Z",
+     "iopub.status.idle": "2022-07-06T18:40:34.876155Z",
+     "shell.execute_reply": "2022-07-06T18:40:34.875146Z"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "105 ms ± 968 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+      "85.4 ms ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
      ]
     }
    ],
@@ -110,13 +138,20 @@
   {
    "cell_type": "code",
    "execution_count": 5,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:34.879157Z",
+     "iopub.status.busy": "2022-07-06T18:40:34.879157Z",
+     "iopub.status.idle": "2022-07-06T18:40:42.930387Z",
+     "shell.execute_reply": "2022-07-06T18:40:42.930079Z"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "14.1 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+      "9.9 ms ± 164 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
      ]
     }
    ],
@@ -138,7 +173,14 @@
   {
    "cell_type": "code",
    "execution_count": 6,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:42.932720Z",
+     "iopub.status.busy": "2022-07-06T18:40:42.932720Z",
+     "iopub.status.idle": "2022-07-06T18:40:42.946939Z",
+     "shell.execute_reply": "2022-07-06T18:40:42.945927Z"
+    }
+   },
    "outputs": [],
    "source": [
     "# And you can also load the JSON that teqp is loading for the pure fluids\n",
@@ -155,7 +197,14 @@
   {
    "cell_type": "code",
    "execution_count": 7,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:42.949975Z",
+     "iopub.status.busy": "2022-07-06T18:40:42.948972Z",
+     "iopub.status.idle": "2022-07-06T18:40:42.962308Z",
+     "shell.execute_reply": "2022-07-06T18:40:42.961547Z"
+    }
+   },
    "outputs": [],
    "source": [
     "# Load the default JSON for the binary interaction parameters\n",
@@ -165,7 +214,14 @@
   {
    "cell_type": "code",
    "execution_count": 8,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:42.965595Z",
+     "iopub.status.busy": "2022-07-06T18:40:42.964944Z",
+     "iopub.status.idle": "2022-07-06T18:40:42.978409Z",
+     "shell.execute_reply": "2022-07-06T18:40:42.977232Z"
+    }
+   },
    "outputs": [
     {
      "data": {
@@ -197,7 +253,14 @@
   {
    "cell_type": "code",
    "execution_count": 9,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:42.981288Z",
+     "iopub.status.busy": "2022-07-06T18:40:42.981288Z",
+     "iopub.status.idle": "2022-07-06T18:40:42.993879Z",
+     "shell.execute_reply": "2022-07-06T18:40:42.993281Z"
+    }
+   },
    "outputs": [
     {
      "data": {
@@ -229,7 +292,14 @@
   {
    "cell_type": "code",
    "execution_count": 10,
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2022-07-06T18:40:42.996338Z",
+     "iopub.status.busy": "2022-07-06T18:40:42.996338Z",
+     "iopub.status.idle": "2022-07-06T18:40:43.423368Z",
+     "shell.execute_reply": "2022-07-06T18:40:43.422895Z"
+    }
+   },
    "outputs": [
     {
      "ename": "ValueError",