Repository: tinosulzer/PyBaMM

Issue:
Discretising "diffusivity times gradient" in Finite Volumes
Flux is of the form `N = -D(c) * grad(c)`

With the current implementation, after discretisation, `grad(c)` will evaluate to a vector of size (n-1,) (at cell edges) while `D(c)` (an arbitrary function of `c`, e.g. `c**2`) will evaluate to a vector of size (n,) (at cell centres), where (n,) is the size of c. So trying to multiply these two will give an error.

The Finite Volume discretisation needs to figure out that we have to take the average of `D(c)` to give `D_edges` which is a vector of size (n-1,) at cell edges. This could probably be done at the `BinaryOperator` stage of `disc.process_symbol`, i.e. if `right` is a `Gradient` object then take the average of `left`.

Discussion:

Patch:
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index c4cae7d6..2acc8bf5 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -194,7 +194,7 @@ $ python run-tests.py --books
 If notebooks fail because of changes to pybamm, it can be a bit of a hassle to debug. In these cases, you can create a temporary export of a notebook's Python content using
 
 ```
-$ python run-tests.py --debook examples/notebooks/notebook-name.ipynb script.py
+$ python run-tests.py --debook examples/notebook-name.ipynb script.py
 ```
 
 
@@ -225,9 +225,9 @@ Next, open a browser, and navigate to your local PyBaMM directory (by typing the
 
 ### Example notebooks
 
-Major PyBaMM features are showcased in [Jupyter notebooks](https://jupyter.org/) stored in the [examples directory](examples/notebooks). Which features are "major" is of course wholly subjective, so please discuss on GitHub first!
+Major PyBaMM features are showcased in [Jupyter notebooks](https://jupyter.org/) stored in the [examples directory](examples). Which features are "major" is of course wholly subjective, so please discuss on GitHub first!
 
-All example notebooks should be listed in [examples/README.md](https://github.com/tinosulzer/pybamm/examples/notebooks/README.md). Please follow the (naming and writing) style of existing notebooks where possible.
+All example notebooks should be listed in [examples/README.md](https://github.com/tinosulzer/pybamm/examples/README.md). Please follow the (naming and writing) style of existing notebooks where possible.
 
 Where possible, notebooks are tested daily. A list of slow notebooks (which time-out and fail tests) is maintained in `.slow-books`, these notebooks will be excluded from daily testing.
 
diff --git a/examples/notebooks/DFN.ipynb b/examples/notebooks/DFN.ipynb
deleted file mode 100644
index 22cb44eb..00000000
--- a/examples/notebooks/DFN.ipynb
+++ /dev/null
@@ -1,146 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Doyle-Fuller-Newman (DFN) Model"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Here we simply state the dimensionless form of the standard DFN model. For more information please see [1] or other classical texts. \n",
-    "\n",
-    "#### Electrolyte Concentration\n",
-    "\n",
-    "$$\n",
-    "\\delta \\epsilon_k \\frac{\\partial c_{\\text{e},k}}{\\partial t} = -\\frac{\\partial N_{\\text{e},k}}{\\partial x} + \\delta \\nu (1-t^+) \\frac{\\partial i_{\\text{e},k}}{ \\partial x}, \\\\ \n",
-    "N_{\\text{e},k} = -\\epsilon_k^b  D_e(c_{\\text{e},k})  \\frac{\\partial c_{\\text{e},k}}{\\partial x}\\\\\n",
-    "    N_{\\text{e,n}}\\big|_{x=0}= N_{\\text{e,p}}\\big|_{x=1}=0, \\\\ \n",
-    "    c_{\\text{e},k}\\big|_{t=0} = 1\n",
-    "$$ \n",
-    "\n",
-    "#### Electrolyte Current \n",
-    "\n",
-    "$$\n",
-    "\\frac{\\partial i_{\\text{e},k}}{\\partial x} = \\begin{cases}\n",
-    "                        G_{k}, \\quad &k=\\text{n, p},\\\\ \n",
-    "                        0, \\quad &k=\\text{s},\n",
-    "                        \\end{cases} \\\\\n",
-    "                        \\delta \\nu i_{\\text{e},k} = \\epsilon_k^b \\sigma_\\text{e}(c_{\\text{e},k}) \\left( - \\Lambda \\frac{\\partial\\phi_{\\text{e},k}}{\\partial x} + 2(1-t^+) \\frac{\\partial}{\\partial x}\\left(\\log(c_{\\text{e},k})\\right)\\right)\\\\\n",
-    "i_{\\text{e,n}}\\big|_{x=0}= i_{\\text{e,p}}\\big|_{x=1}=0,\n",
-    "$$\n",
-    "\n",
-    "#### Electrode Current \n",
-    "\n",
-    "\\begin{align} \n",
-    "\\frac{\\partial i_{k}}{\\partial x} &= -G_{k},\\\\\n",
-    "                        i_{k} &= -\\sigma_k\\frac{\\partial\\phi_{k}}{\\partial x}\n",
-    "\\end{align}\n",
-    "$$\n",
-    "i_{\\text{n}}\\big|_{x=0}= i_{\\text{p}}\\big|_{x=1}=\\mathcal{I},\n",
-    "$$\n",
-    "\n",
-    "\n",
-    "#### Particles \n",
-    "\n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = G_k \\\\ \n",
-    "c_k\\big|_{t=0} = c_{k,0}\n",
-    "$$ \n",
-    "\n",
-    "\n",
-    "#### Electrochemical Reactions\n",
-    "\n",
-    "$$\n",
-    "G_k = g_k \\sinh\\left(\\frac{\\Lambda \\eta_k}{2} \\right), \\quad \n",
-    "g_k = m_k \\hat{C}_k c_k^{1/2} (1-c_k)^{1/2}c_{ek}^{1/2}\\big|_{r_k=1}, \\quad\n",
-    "\\eta_k = \\phi_k - \\phi_{ek} - U_k(c_k)\\big|_{r_k=1},\n",
-    "$$\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We provide here the dimensionless parameters in the DFN in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities. We have\n",
-    "\\begin{gather}\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\quad\n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*}, \\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*}, \\\\\n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}, \\quad \n",
-    "\\delta = \\frac{L^* I^*}{F^* D_{\\text{e,typ}}^* c_{\\text{n,max}}^*}, \\quad \n",
-    "\\nu = \\frac{c_{\\text{n, max}}^*}{c_{\\text{e,typ}}^*} , \\quad \n",
-    "\\sigma_k = \\frac{\\sigma_k^* \\Phi^*}{I^* L_x^*}\n",
-    "\\end{gather} \n",
-    "\n",
-    "We note that the dimensionless parameters $\\epsilon_k$, $b$, and $t^+$ are already provided in the parameter file LCO.csv.  "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## References\n",
-    "[1] Doyle, Marc, Thomas F. Fuller, and John Newman. \"Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell.\" Journal of the Electrochemical Society 140.6 (1993): 1526-1533."
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/README.md b/examples/notebooks/README.md
index dcd3dd57..758c0dd6 100644
--- a/examples/notebooks/README.md
+++ b/examples/notebooks/README.md
@@ -1,67 +1,13 @@
 # Examples
 
-[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/tinosulzer/PyBaMM/master)
-
-This page contains a number of examples showing how to use PyBaMM.
+This page contains (or will eventually contain) a number of examples showing how to use PyBaMM.
 
 Each example was created as a _Jupyter notebook_ (http://jupyter.org/).
-These notebooks can be downloaded and used locally by running
-```
-$ jupyter notebook
-```
-from your local PyBaMM repository, or used online through [Binder](https://mybinder.org/v2/gh/tinosulzer/PyBaMM/master), or you can simply copy/paste the relevant code.
+These notebooks can be downloaded and used, or you can simply copy/paste the
+relevant code.
 
 ## Getting started
 
-The easiest way to start with PyBaMM is by running and comparing some of the inbuilt models:
-- [Run a pre-defined model](./run-model.ipynb)
-- [Compare models](./compare-models.ipynb)
-
-It is also easy to new models or change the setting that are used:
-- [Add a model](./add-model.ipynb)
-- [Change the settings](./change-settings.ipynb) (parameters, discretisation or solver)
-
-For more advanced usage, new sets of parameters, discretisations and solvers can be added:
-- [Add parameters](./add-parameters.ipynb)
-- [Add a discretisation](./add-discretisation.ipynb)
-- [Add a solver](./add-solver.ipynb)
-
-## Expression tree structure
-
-PyBaMM is built around an expression tree structure, similar to FEniCS/Firedrake's [UFL](https://fenics.readthedocs.io/projects/ufl/en/latest/). The following notebooks explains how this works, from model creation to solution.
-1. [Overview](./expression-tree/overview.ipynb)
-2. [Expression tree](./expression-tree/expression-tree.ipynb)
-3. [Setting parameters](./expression-tree/parameters.ipynb)
-4. [Discretising a model](./expression-tree/discretisation.ipynb)
-5. [Solving a model](./expression-tree/solver.ipynb)
-
-### Models
-
-The following models are implemented and can easily be [used](./run-model.ipynb) or [compared](./compare-models.ipynb). We always welcome [new models](./add-model.ipynb)!
-
-#### Lithium-ion models
-
-- [Single-Particle Model](./models/SPM.ipynb)
-- [Single-Particle Model with electrolyte](./models/SPMe.ipynb)
-- [Doyle-Fuller-Newman Model](./models/DFN.ipynb)
-
-#### Lead-acid models
-
-- [Full porous-electrode](./models/lead-acid-full.ipynb)
-- [Leading-Order Quasi-Static](./models/lead-acid-LOQS.ipynb)
-- [First-Order Quasi-Static](./models/lead-acid-FOQS.ipynb)
-- [Composite](./models/lead-acid-composite.ipynb)
-
-### Discretisations
-
-The following discretisation is implemented
-- [Finite Volumes](./discretisations/finite-volumes.ipynb)
-
-See [here](./add-discretisation.ipynb) for instructions on adding new discretisations.
-
-### Solvers
-
-The following solver is implemented
-- [Inbuilt SciPy solver](./solvers/scipy-integrate.ipynb)
-
-See [here](./add-solver.ipynb) for instructions on adding new solvers.
+- [Solve a simple diffusion equation](./simple-diffusion.ipynb)
+<!-- - [Compare models](./compare-models.ipynb)
+- [Add a model](./add-model.ipynb) -->
diff --git a/examples/notebooks/SPM.ipynb b/examples/notebooks/SPM.ipynb
deleted file mode 100644
index 02231806..00000000
--- a/examples/notebooks/SPM.ipynb
+++ /dev/null
@@ -1,97 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Single Particle Model (SPM) "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The SPM consists of two spherically symmetric diffusion equations: one within a representative negative particle ($k=\\text{n}$) and one within a representative positive particle ($k=\\text{p}$). In the centre of the particle the classical no-flux condition is imposed. Since the SPM assumes that all particles in an electrode behave in exactly the same way, the flux on the surface of a particle is simply the current divided by the thickness of the electrode. We shall use $r_k\\in[0,1]$ for the radial coordinate of the particle in electrode $k$. The concentration of lithium in electrode $k$ is denoted $c_k$ and the current is denoted by $\\mathcal{I}$. All parameters in the model stated here are dimensionless and are given in terms of dimensional parameters at the end of this notebook.The model equations for the SPM are then: \n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = \\begin{cases} \n",
-    "    \\frac{\\mathcal{I}}{L_n}, \\quad k=\\text{n}, \\\\ \n",
-    "    -\\frac{\\mathcal{I}}{L_p}, \\quad k=\\text{p},\n",
-    "\\end{cases} ,\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, p}\\}$.\n",
-    "\n",
-    "### Voltage Expression\n",
-    "The terminal voltage is obtained from the expression: \n",
-    "$$\n",
-    "V = U_p(c_p)\\big|_{r_p=1}-U_n(c_n)\\big|_{r_n=1}-\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{g_p L_p} \\right)  -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{g_n L_n} \\right)\n",
-    "$$\n",
-    "with the exchange current densities given by\n",
-    "$$\n",
-    "g_k = m_k \\hat{C}_k (c_k)^{1/2} (1-c_k)^{1/2}, \\quad \\text{for } k\\in\\{\\text{n, p}\\}.\n",
-    "$$"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters \n",
-    "We provide here the dimensionless parameters in the SPM in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities here. We have\n",
-    "$$\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\\\\n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*},\\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*} \\quad \n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}\n",
-    "$$"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/SPMe.ipynb b/examples/notebooks/SPMe.ipynb
deleted file mode 100644
index 253aca69..00000000
--- a/examples/notebooks/SPMe.ipynb
+++ /dev/null
@@ -1,138 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Single Particle Model with Electrolyte (SPMe)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The SPMe consists of two spherically symmetric diffusion equations: one within a representative negative particle ($k=\\text{n}$) and one within a representative positive particle ($k=\\text{p}$) and a diffusion equation within the electrolyte. In the centre of each particle the classical no-flux condition is imposed and the flux on the surface of the particle is simply the current divided by the thickness of the electrode. Since lithium is transferred between the electrolyte and particles, the flux through the particle surface also enters the electrolyte diffusion equation as a source/sink. The boundary conditions on the electrolyte model are no flux at each current collector to represent no transfer of lithium between the electrolyte and current collectors. We shall use $r_k\\in[0,1]$ for the radial coordinate of the particle in electrode $k$ and $x\\in[0,1]$ for the macroscopic through cell distance. The concentration of lithium in electrode $k$ is denoted $c_k$ and the concentration of lithium ions in the electrolyte within the $k$ electrode region is denoted $c_{\\text{e},k}$. We also denote the flux of lithium ions in the electrolyte within the $k$ electrode region by $N_{\\text{e},k}$ and the current by $\\mathcal{I}$. All parameters in the model stated here are dimensionless and are given in terms of dimensional parameters at the end of this notebook. The model equations of the SPMe are then: \n",
-    "\n",
-    "#### Particles: \n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = \\begin{cases} \n",
-    "    \\frac{\\mathcal{I}}{L_n}, \\quad k=\\text{n}, \\\\ \n",
-    "    -\\frac{\\mathcal{I}}{L_p}, \\quad k=\\text{p},\n",
-    "\\end{cases},\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, p}\\}$. \n",
-    "\n",
-    "#### Electrolyte: \n",
-    "$$ \n",
-    "\\delta \\epsilon_k \\frac{\\partial c_{\\text{e},k}}{\\partial t} = - \\frac{\\partial N_{\\text{e},k}}{\\partial c} + \n",
-    "\\begin{cases} \n",
-    "    \\nu (1-t^+) \\frac{\\mathcal{I}}{L_n}, \\quad &\\text{in negative electrode region}, \\\\\n",
-    "    0, \\quad &\\text{in separator region}, \\\\ \n",
-    "    -\\nu(1-t^+)\\frac{\\mathcal{I}}{L_p}, \\quad &\\text{in positive electrode region}\n",
-    "\\end{cases}, \\\\ \n",
-    "N_{\\text{e},k} = - \\epsilon_k^b D_e(1) \\frac{\\partial c_{\\text{e},k}}{\\partial x}, \\\\ \n",
-    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=1}=0\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, s, p}\\}$. \n",
-    "\n",
-    "### Voltage Expression\n",
-    "The terminal voltage is obtained from the expression: \n",
-    "$$ \n",
-    "V = U_{\\text{eq}} + \\eta_r + \\eta_c + \\Delta\\Phi_{\\text{Elec}} + \\Delta\\Phi_{\\text{Solid}}\n",
-    "$$ \n",
-    "where \n",
-    "\\begin{align}\n",
-    "& U_{\\text{eq}} = U_p(c_p)\\big|_{r_p=1} - U_n(c_n)\\big|_{r_n=1}, \\\\ \n",
-    "     &\\eta_{r} = -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{\\bar{g}_p L_n}\\right)\n",
-    "     -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{\\bar{g}_n L_n}\\right), \\\\\n",
-    "     &\\eta_c =  \\frac{2 \\delta (1-t^+)}{\\Lambda}\\left(\\bar{c}_{ep} - \\bar{c}_{en}\\right), \\\\\n",
-    "     \\label{eqn:SPMe:negative_exchange}\n",
-    "     &\\bar{g}_n =  \\frac{1}{L_n}\\int_0^{L_n} m_n \\hat{C}_n (c_n)^{1/2}(1-c_n)^{1/2} (1+\\delta c_{en})^{1/2} \\, \\text{d}x, \\\\\n",
-    "     \\label{eqn:SPMe:positive_exchange}\n",
-    "    &\\bar{g}_p =  \\frac{1}{L_p}\\int_{1-L_p}^{1} m_p \\hat{C}_p  (c_p)^{1/2}(1-c_p)^{1/2} (1+\\delta c_{ep} )^{1/2} \\, \\text{d}x, \\\\\n",
-    "    &\\Delta \\Phi_{\\text{Elec}}= -\\frac{\\delta \\nu \\mathcal{I}}{\\Lambda \\sigma_e(1)}\\left(\\frac{L_n}{3\\epsilon_n^b} + \\frac{L_s}{\\epsilon_{s}^b} + \\frac{L_p}{3\\epsilon_p^b} \\right),\n",
-    "     \\label{eqn:SPMe:electrolyte_ohmic_losses} \\\\\n",
-    "     &\\Delta \\Phi_{\\text{Solid}} =  -\\frac{\\mathcal{I}}{3}\\left(\\frac{L_p}{\\sigma_p} + \\frac{L_n}{\\sigma_n} \\right), \n",
-    "\\end{align} \n",
-    "with \n",
-    "\\begin{equation} \n",
-    "    \\bar{c}_{\\text{e},n} =  \\frac{1}{L_n}\\int_0^{L_n} c_{\\text{e},n} \\, \\text{d}x, \\quad\n",
-    "    \\bar{c}_{\\text{e},p} =  \\frac{1}{L_p}\\int_{1-L_p}^{1} c_{\\text{e},p} \\, \\text{d}x,\n",
-    "\\end{equation} \n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We provide here the dimensionless parameters in the SPMe in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities. We have\n",
-    "\\begin{gather}\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\quad \n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*}, \\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*} \\\\\n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}, \\quad \n",
-    "\\delta = \\frac{L^* I^*}{F^* D_{\\text{e,typ}}^* c_{\\text{n,max}}^*}, \\quad \n",
-    "\\nu = \\frac{c_{\\text{n, max}}^*}{c_{\\text{e,typ}}^*} , \\quad \n",
-    "\\sigma_k = \\frac{\\sigma_k^* \\Phi^*}{I^* L_x^*}\n",
-    "\\end{gather} \n",
-    "\n",
-    "We note that the dimensionless parameters $\\epsilon_k$, $b$, and $t^+$ are already provided in the parameter file LCO.csv.  "
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/add-discretisation.ipynb b/examples/notebooks/add-discretisation.ipynb
deleted file mode 100644
index 5cdaebcc..00000000
--- a/examples/notebooks/add-discretisation.ipynb
+++ /dev/null
@@ -1,41 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Adding a Discretisation\n",
-    "\n",
-    "Before adding a new discretisation, please read the [contribution guidelines](https://github.com/tinosulzer/PyBaMM/blob/master/CONTRIBUTING.md)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/add-model.ipynb b/examples/notebooks/add-model.ipynb
deleted file mode 100644
index b4ca6e16..00000000
--- a/examples/notebooks/add-model.ipynb
+++ /dev/null
@@ -1,41 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Adding a Model\n",
-    "\n",
-    "Before adding a new model, please read the [contribution guidelines](https://github.com/tinosulzer/PyBaMM/blob/master/CONTRIBUTING.md)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/add-parameters.ipynb b/examples/notebooks/add-parameters.ipynb
deleted file mode 100644
index 250e19c4..00000000
--- a/examples/notebooks/add-parameters.ipynb
+++ /dev/null
@@ -1,34 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Adding Parameters\n",
-    "\n",
-    "Before adding new parameters, please read the [contribution guidelines](https://github.com/tinosulzer/PyBaMM/blob/master/CONTRIBUTING.md)"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/add-solver.ipynb b/examples/notebooks/add-solver.ipynb
deleted file mode 100644
index e0cb0f7c..00000000
--- a/examples/notebooks/add-solver.ipynb
+++ /dev/null
@@ -1,41 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Adding a Solver\n",
-    "\n",
-    "Before adding a new solver, please read the [contribution guidelines](https://github.com/tinosulzer/PyBaMM/blob/master/CONTRIBUTING.md)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/change-settings.ipynb b/examples/notebooks/change-settings.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/change-settings.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/compare-models.ipynb b/examples/notebooks/compare-models.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/compare-models.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/discretisations/finite-volumes.ipynb b/examples/notebooks/discretisations/finite-volumes.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/discretisations/finite-volumes.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/expression-tree/discretisation.ipynb b/examples/notebooks/expression-tree/discretisation.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/expression-tree/discretisation.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/expression-tree/expression-tree.ipynb b/examples/notebooks/expression-tree/expression-tree.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/expression-tree/expression-tree.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/expression-tree/overview.ipynb b/examples/notebooks/expression-tree/overview.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/expression-tree/overview.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/expression-tree/parameters.ipynb b/examples/notebooks/expression-tree/parameters.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/expression-tree/parameters.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/expression-tree/solver.ipynb b/examples/notebooks/expression-tree/solver.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/expression-tree/solver.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/DFN.ipynb b/examples/notebooks/models/DFN.ipynb
deleted file mode 100644
index 22cb44eb..00000000
--- a/examples/notebooks/models/DFN.ipynb
+++ /dev/null
@@ -1,146 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Doyle-Fuller-Newman (DFN) Model"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Here we simply state the dimensionless form of the standard DFN model. For more information please see [1] or other classical texts. \n",
-    "\n",
-    "#### Electrolyte Concentration\n",
-    "\n",
-    "$$\n",
-    "\\delta \\epsilon_k \\frac{\\partial c_{\\text{e},k}}{\\partial t} = -\\frac{\\partial N_{\\text{e},k}}{\\partial x} + \\delta \\nu (1-t^+) \\frac{\\partial i_{\\text{e},k}}{ \\partial x}, \\\\ \n",
-    "N_{\\text{e},k} = -\\epsilon_k^b  D_e(c_{\\text{e},k})  \\frac{\\partial c_{\\text{e},k}}{\\partial x}\\\\\n",
-    "    N_{\\text{e,n}}\\big|_{x=0}= N_{\\text{e,p}}\\big|_{x=1}=0, \\\\ \n",
-    "    c_{\\text{e},k}\\big|_{t=0} = 1\n",
-    "$$ \n",
-    "\n",
-    "#### Electrolyte Current \n",
-    "\n",
-    "$$\n",
-    "\\frac{\\partial i_{\\text{e},k}}{\\partial x} = \\begin{cases}\n",
-    "                        G_{k}, \\quad &k=\\text{n, p},\\\\ \n",
-    "                        0, \\quad &k=\\text{s},\n",
-    "                        \\end{cases} \\\\\n",
-    "                        \\delta \\nu i_{\\text{e},k} = \\epsilon_k^b \\sigma_\\text{e}(c_{\\text{e},k}) \\left( - \\Lambda \\frac{\\partial\\phi_{\\text{e},k}}{\\partial x} + 2(1-t^+) \\frac{\\partial}{\\partial x}\\left(\\log(c_{\\text{e},k})\\right)\\right)\\\\\n",
-    "i_{\\text{e,n}}\\big|_{x=0}= i_{\\text{e,p}}\\big|_{x=1}=0,\n",
-    "$$\n",
-    "\n",
-    "#### Electrode Current \n",
-    "\n",
-    "\\begin{align} \n",
-    "\\frac{\\partial i_{k}}{\\partial x} &= -G_{k},\\\\\n",
-    "                        i_{k} &= -\\sigma_k\\frac{\\partial\\phi_{k}}{\\partial x}\n",
-    "\\end{align}\n",
-    "$$\n",
-    "i_{\\text{n}}\\big|_{x=0}= i_{\\text{p}}\\big|_{x=1}=\\mathcal{I},\n",
-    "$$\n",
-    "\n",
-    "\n",
-    "#### Particles \n",
-    "\n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = G_k \\\\ \n",
-    "c_k\\big|_{t=0} = c_{k,0}\n",
-    "$$ \n",
-    "\n",
-    "\n",
-    "#### Electrochemical Reactions\n",
-    "\n",
-    "$$\n",
-    "G_k = g_k \\sinh\\left(\\frac{\\Lambda \\eta_k}{2} \\right), \\quad \n",
-    "g_k = m_k \\hat{C}_k c_k^{1/2} (1-c_k)^{1/2}c_{ek}^{1/2}\\big|_{r_k=1}, \\quad\n",
-    "\\eta_k = \\phi_k - \\phi_{ek} - U_k(c_k)\\big|_{r_k=1},\n",
-    "$$\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We provide here the dimensionless parameters in the DFN in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities. We have\n",
-    "\\begin{gather}\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\quad\n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*}, \\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*}, \\\\\n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}, \\quad \n",
-    "\\delta = \\frac{L^* I^*}{F^* D_{\\text{e,typ}}^* c_{\\text{n,max}}^*}, \\quad \n",
-    "\\nu = \\frac{c_{\\text{n, max}}^*}{c_{\\text{e,typ}}^*} , \\quad \n",
-    "\\sigma_k = \\frac{\\sigma_k^* \\Phi^*}{I^* L_x^*}\n",
-    "\\end{gather} \n",
-    "\n",
-    "We note that the dimensionless parameters $\\epsilon_k$, $b$, and $t^+$ are already provided in the parameter file LCO.csv.  "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## References\n",
-    "[1] Doyle, Marc, Thomas F. Fuller, and John Newman. \"Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell.\" Journal of the Electrochemical Society 140.6 (1993): 1526-1533."
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/SPM.ipynb b/examples/notebooks/models/SPM.ipynb
deleted file mode 100644
index 02231806..00000000
--- a/examples/notebooks/models/SPM.ipynb
+++ /dev/null
@@ -1,97 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Single Particle Model (SPM) "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The SPM consists of two spherically symmetric diffusion equations: one within a representative negative particle ($k=\\text{n}$) and one within a representative positive particle ($k=\\text{p}$). In the centre of the particle the classical no-flux condition is imposed. Since the SPM assumes that all particles in an electrode behave in exactly the same way, the flux on the surface of a particle is simply the current divided by the thickness of the electrode. We shall use $r_k\\in[0,1]$ for the radial coordinate of the particle in electrode $k$. The concentration of lithium in electrode $k$ is denoted $c_k$ and the current is denoted by $\\mathcal{I}$. All parameters in the model stated here are dimensionless and are given in terms of dimensional parameters at the end of this notebook.The model equations for the SPM are then: \n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = \\begin{cases} \n",
-    "    \\frac{\\mathcal{I}}{L_n}, \\quad k=\\text{n}, \\\\ \n",
-    "    -\\frac{\\mathcal{I}}{L_p}, \\quad k=\\text{p},\n",
-    "\\end{cases} ,\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, p}\\}$.\n",
-    "\n",
-    "### Voltage Expression\n",
-    "The terminal voltage is obtained from the expression: \n",
-    "$$\n",
-    "V = U_p(c_p)\\big|_{r_p=1}-U_n(c_n)\\big|_{r_n=1}-\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{g_p L_p} \\right)  -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{g_n L_n} \\right)\n",
-    "$$\n",
-    "with the exchange current densities given by\n",
-    "$$\n",
-    "g_k = m_k \\hat{C}_k (c_k)^{1/2} (1-c_k)^{1/2}, \\quad \\text{for } k\\in\\{\\text{n, p}\\}.\n",
-    "$$"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters \n",
-    "We provide here the dimensionless parameters in the SPM in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities here. We have\n",
-    "$$\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\\\\n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*},\\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*} \\quad \n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}\n",
-    "$$"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/SPMe.ipynb b/examples/notebooks/models/SPMe.ipynb
deleted file mode 100644
index 253aca69..00000000
--- a/examples/notebooks/models/SPMe.ipynb
+++ /dev/null
@@ -1,138 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Single Particle Model with Electrolyte (SPMe)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Model Equations"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The SPMe consists of two spherically symmetric diffusion equations: one within a representative negative particle ($k=\\text{n}$) and one within a representative positive particle ($k=\\text{p}$) and a diffusion equation within the electrolyte. In the centre of each particle the classical no-flux condition is imposed and the flux on the surface of the particle is simply the current divided by the thickness of the electrode. Since lithium is transferred between the electrolyte and particles, the flux through the particle surface also enters the electrolyte diffusion equation as a source/sink. The boundary conditions on the electrolyte model are no flux at each current collector to represent no transfer of lithium between the electrolyte and current collectors. We shall use $r_k\\in[0,1]$ for the radial coordinate of the particle in electrode $k$ and $x\\in[0,1]$ for the macroscopic through cell distance. The concentration of lithium in electrode $k$ is denoted $c_k$ and the concentration of lithium ions in the electrolyte within the $k$ electrode region is denoted $c_{\\text{e},k}$. We also denote the flux of lithium ions in the electrolyte within the $k$ electrode region by $N_{\\text{e},k}$ and the current by $\\mathcal{I}$. All parameters in the model stated here are dimensionless and are given in terms of dimensional parameters at the end of this notebook. The model equations of the SPMe are then: \n",
-    "\n",
-    "#### Particles: \n",
-    "$$\n",
-    "\\frac{\\partial c_k}{\\partial t} = \\gamma_k \\frac{\\partial}{\\partial r_k}\\left( r_k^2 \\frac{\\partial c_k}{\\partial r_k} \\right) \\\\ \n",
-    "\\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=0} = 0, \\quad -\\beta_k \\hat{C}_k \\gamma_k \\frac{\\partial c_k}{\\partial r_k}\\bigg|_{r_k=1} = \\begin{cases} \n",
-    "    \\frac{\\mathcal{I}}{L_n}, \\quad k=\\text{n}, \\\\ \n",
-    "    -\\frac{\\mathcal{I}}{L_p}, \\quad k=\\text{p},\n",
-    "\\end{cases},\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, p}\\}$. \n",
-    "\n",
-    "#### Electrolyte: \n",
-    "$$ \n",
-    "\\delta \\epsilon_k \\frac{\\partial c_{\\text{e},k}}{\\partial t} = - \\frac{\\partial N_{\\text{e},k}}{\\partial c} + \n",
-    "\\begin{cases} \n",
-    "    \\nu (1-t^+) \\frac{\\mathcal{I}}{L_n}, \\quad &\\text{in negative electrode region}, \\\\\n",
-    "    0, \\quad &\\text{in separator region}, \\\\ \n",
-    "    -\\nu(1-t^+)\\frac{\\mathcal{I}}{L_p}, \\quad &\\text{in positive electrode region}\n",
-    "\\end{cases}, \\\\ \n",
-    "N_{\\text{e},k} = - \\epsilon_k^b D_e(1) \\frac{\\partial c_{\\text{e},k}}{\\partial x}, \\\\ \n",
-    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=1}=0\n",
-    "$$ \n",
-    "for $k\\in\\{\\text{n, s, p}\\}$. \n",
-    "\n",
-    "### Voltage Expression\n",
-    "The terminal voltage is obtained from the expression: \n",
-    "$$ \n",
-    "V = U_{\\text{eq}} + \\eta_r + \\eta_c + \\Delta\\Phi_{\\text{Elec}} + \\Delta\\Phi_{\\text{Solid}}\n",
-    "$$ \n",
-    "where \n",
-    "\\begin{align}\n",
-    "& U_{\\text{eq}} = U_p(c_p)\\big|_{r_p=1} - U_n(c_n)\\big|_{r_n=1}, \\\\ \n",
-    "     &\\eta_{r} = -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{\\bar{g}_p L_n}\\right)\n",
-    "     -\\frac{2}{\\Lambda}\\sinh^{-1}\\left(\\frac{\\mathcal{I}}{\\bar{g}_n L_n}\\right), \\\\\n",
-    "     &\\eta_c =  \\frac{2 \\delta (1-t^+)}{\\Lambda}\\left(\\bar{c}_{ep} - \\bar{c}_{en}\\right), \\\\\n",
-    "     \\label{eqn:SPMe:negative_exchange}\n",
-    "     &\\bar{g}_n =  \\frac{1}{L_n}\\int_0^{L_n} m_n \\hat{C}_n (c_n)^{1/2}(1-c_n)^{1/2} (1+\\delta c_{en})^{1/2} \\, \\text{d}x, \\\\\n",
-    "     \\label{eqn:SPMe:positive_exchange}\n",
-    "    &\\bar{g}_p =  \\frac{1}{L_p}\\int_{1-L_p}^{1} m_p \\hat{C}_p  (c_p)^{1/2}(1-c_p)^{1/2} (1+\\delta c_{ep} )^{1/2} \\, \\text{d}x, \\\\\n",
-    "    &\\Delta \\Phi_{\\text{Elec}}= -\\frac{\\delta \\nu \\mathcal{I}}{\\Lambda \\sigma_e(1)}\\left(\\frac{L_n}{3\\epsilon_n^b} + \\frac{L_s}{\\epsilon_{s}^b} + \\frac{L_p}{3\\epsilon_p^b} \\right),\n",
-    "     \\label{eqn:SPMe:electrolyte_ohmic_losses} \\\\\n",
-    "     &\\Delta \\Phi_{\\text{Solid}} =  -\\frac{\\mathcal{I}}{3}\\left(\\frac{L_p}{\\sigma_p} + \\frac{L_n}{\\sigma_n} \\right), \n",
-    "\\end{align} \n",
-    "with \n",
-    "\\begin{equation} \n",
-    "    \\bar{c}_{\\text{e},n} =  \\frac{1}{L_n}\\int_0^{L_n} c_{\\text{e},n} \\, \\text{d}x, \\quad\n",
-    "    \\bar{c}_{\\text{e},p} =  \\frac{1}{L_p}\\int_{1-L_p}^{1} c_{\\text{e},p} \\, \\text{d}x,\n",
-    "\\end{equation} \n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Example"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pybamm"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Dimensionless Parameters"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We provide here the dimensionless parameters in the SPMe in terms of the dimensional parameters in LCO.csv. We use a superscript * to indicate dimensional quanitities. We have\n",
-    "\\begin{gather}\n",
-    "\\gamma_k = \\frac{F^* c_{\\text{n,max}}^* L_x^* D_k^*}{I^* (R_k^*)^2}, \\quad \n",
-    "\\beta_k= a_k^*R_k^*, \\quad \n",
-    "\\hat{C}_k=\\frac{c_\\text{k,max}^*}{c_{\\text{n,max}}^*} \\quad \n",
-    "L_k = \\frac{L_k^*}{L_x}, \\quad \n",
-    "\\Lambda = \\frac{F^* \\Phi^*}{R^* T_{\\text{ref}}^*}, \\quad \n",
-    "c_{k,0} = \\frac{c_{k,0}^*}{c_{k,\\text{max}}^*} \\\\\n",
-    "m_k = \\frac{m_k^* c_{\\text{n,max}}^* (c_{\\text{e,typ}}^*)^{1/2} a_k^*}{I^*}, \\quad \n",
-    "\\delta = \\frac{L^* I^*}{F^* D_{\\text{e,typ}}^* c_{\\text{n,max}}^*}, \\quad \n",
-    "\\nu = \\frac{c_{\\text{n, max}}^*}{c_{\\text{e,typ}}^*} , \\quad \n",
-    "\\sigma_k = \\frac{\\sigma_k^* \\Phi^*}{I^* L_x^*}\n",
-    "\\end{gather} \n",
-    "\n",
-    "We note that the dimensionless parameters $\\epsilon_k$, $b$, and $t^+$ are already provided in the parameter file LCO.csv.  "
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/lead-acid-composite.ipynb b/examples/notebooks/models/lead-acid-composite.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/models/lead-acid-composite.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/lead-acid-foqs.ipynb b/examples/notebooks/models/lead-acid-foqs.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/models/lead-acid-foqs.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/lead-acid-full.ipynb b/examples/notebooks/models/lead-acid-full.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/models/lead-acid-full.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/models/lead-acid-loqs.ipynb b/examples/notebooks/models/lead-acid-loqs.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/models/lead-acid-loqs.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/run-model.ipynb b/examples/notebooks/run-model.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/run-model.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/notebooks/simple-diffusion.ipynb b/examples/notebooks/simple-diffusion.ipynb
index a246cf8f..53420fb9 100644
--- a/examples/notebooks/simple-diffusion.ipynb
+++ b/examples/notebooks/simple-diffusion.ipynb
@@ -32,8 +32,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "import os\n",
@@ -52,7 +54,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 13,
    "metadata": {
     "collapsed": true
    },
@@ -71,7 +73,9 @@
   {
    "cell_type": "code",
    "execution_count": 14,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "param = pybamm.ParameterValues()"
@@ -87,7 +91,9 @@
   {
    "cell_type": "code",
    "execution_count": 15,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "tsteps = 100\n",
@@ -106,7 +112,9 @@
   {
    "cell_type": "code",
    "execution_count": 16,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "sim = pybamm.Simulation(model, param, mesh, name=\"Electrolyte diffusion\")"
@@ -142,7 +150,9 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -513,7 +523,9 @@
   {
    "cell_type": "code",
    "execution_count": 8,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -549,18 +561,20 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "ename": "NameError",
      "evalue": "name 'sim' is not defined",
-     "output_type": "error",
      "traceback": [
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
       "\u001b[0;32m<ipython-input-11-f0e817a6d484>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msave\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
       "\u001b[0;31mNameError\u001b[0m: name 'sim' is not defined"
-     ]
+     ],
+     "output_type": "error"
     }
    ],
    "source": [
@@ -570,21 +584,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (anaconda)",
    "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.6.7"
+   "name": "myenv"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/solvers/scipy-integrate.ipynb b/examples/notebooks/solvers/scipy-integrate.ipynb
deleted file mode 100644
index 72c3edc8..00000000
--- a/examples/notebooks/solvers/scipy-integrate.ipynb
+++ /dev/null
@@ -1,32 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "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.6.7"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index 3cad9a12..355fa239 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -7,7 +7,6 @@ import pybamm
 
 import numpy as np
 import numbers
-import copy
 
 
 class BaseDiscretisation(object):
@@ -200,15 +199,22 @@ class BaseDiscretisation(object):
 
             return new_symbol
 
+        elif isinstance(symbol, pybamm.Scalar):
+            return pybamm.Scalar(symbol.value, domain=symbol.domain)
+
+        elif isinstance(symbol, pybamm.Array):
+            return symbol.__class__(symbol.entries, domain=symbol.domain)
+
         else:
+            raise NotImplementedError
             # hack to copy the symbol but without a parent
             # (building tree from bottom up)
             # simply setting new_symbol.parent = None, after copying, raises a TreeError
-            parent = symbol.parent
-            symbol.parent = None
-            new_symbol = copy.copy(symbol)
-            symbol.parent = parent
-            return new_symbol
+            # parent = symbol.parent
+            # symbol.parent = None
+            # new_symbol = copy.copy(symbol)
+            # symbol.parent = parent
+            # return new_symbol
 
     def process_binary_operators(self, bin_op, y_slices, boundary_conditions):
         """Discretise binary operators in model equations.
@@ -237,17 +243,25 @@ class BaseDiscretisation(object):
         left, right = bin_op.children
         new_left = self.process_symbol(left, y_slices, boundary_conditions)
         new_right = self.process_symbol(right, y_slices, boundary_conditions)
-
         # Post-processing to make sure discretised dimensions match
         # If neither child has gradients, or both children have gradients
         # no need to do any averaging
-        if left.has_gradient() == right.has_gradient():
+        if (
+            left.has_gradient_and_not_divergence()
+            == right.has_gradient_and_not_divergence()
+        ):
             pass
         # If only left child has gradient, compute diffusivity for right child
-        elif left.has_gradient() and not right.has_gradient():
+        elif (
+            left.has_gradient_and_not_divergence()
+            and not right.has_gradient_and_not_divergence()
+        ):
             new_right = self.compute_diffusivity(new_right)
         # If only right child has gradient, compute diffusivity for left child
-        elif right.has_gradient() and not left.has_gradient():
+        elif (
+            right.has_gradient_and_not_divergence()
+            and not left.has_gradient_and_not_divergence()
+        ):
             new_left = self.compute_diffusivity(new_left)
         # Return new binary operator with appropriate class
         return bin_op.__class__(new_left, new_right)
diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py
index 162d1095..3c247db5 100644
--- a/pybamm/expression_tree/array.py
+++ b/pybamm/expression_tree/array.py
@@ -29,6 +29,10 @@ class Array(pybamm.Symbol):
         super().__init__(name, domain=domain)
         self._entries = entries
 
+    @property
+    def entries(self):
+        return self._entries
+
     @property
     def ndim(self):
         """ returns the number of dimensions of the tensor"""
diff --git a/pybamm/expression_tree/matrix.py b/pybamm/expression_tree/matrix.py
index 3d6fca1b..7df8b706 100644
--- a/pybamm/expression_tree/matrix.py
+++ b/pybamm/expression_tree/matrix.py
@@ -21,5 +21,5 @@ class Matrix(pybamm.Array):
 
     """
 
-    def __init__(self, entries, name=None):
-        super().__init__(entries, name=name)
+    def __init__(self, entries, name=None, domain=[]):
+        super().__init__(entries, name=name, domain=domain)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 942713d3..5b9ae47c 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -308,12 +308,16 @@ class Symbol(anytree.NodeMixin):
         """Returns True if equation has spatial derivatives (grad or div)."""
         return self.has_gradient() or self.has_divergence()
 
+    def has_gradient_and_not_divergence(self):
+        """Returns True if equation has a Gradient term and not Divergence term."""
+        return self.has_gradient() and not self.has_divergence()
+
     def has_gradient(self):
-        """Returns True if equation has a Gradient."""
+        """Returns True if equation has a Gradient term."""
         return any([isinstance(symbol, pybamm.Gradient) for symbol in self.pre_order()])
 
     def has_divergence(self):
-        """Returns True if equation has a Divergence."""
+        """Returns True if equation has a Divergence term."""
         return any(
             [isinstance(symbol, pybamm.Divergence) for symbol in self.pre_order()]
         )


Test patch:
diff --git a/run-tests.py b/run-tests.py
index 525d1487..0f848fe3 100755
--- a/run-tests.py
+++ b/run-tests.py
@@ -277,10 +277,6 @@ if __name__ == "__main__":
         help="Export a Jupyter notebook to a Python file for manual testing.",
     )
     # Doctests
-    parser.add_argument(
-        "--flake8", action="store_true", help="Run flake8 to check for style issues"
-    )
-    # Doctests
     parser.add_argument(
         "--doctest",
         action="store_true",
@@ -312,10 +308,6 @@ if __name__ == "__main__":
     if args.nosub:
         has_run = True
         run_unit_tests()
-    # Flake8
-    if args.flake8:
-        has_run = True
-        run_flake8()
     # Doctests
     if args.doctest:
         has_run = True
diff --git a/tests/__init__.py b/tests/__init__.py
index 0dacf942..b5d1a622 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -5,4 +5,4 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
-from test_models.standard_model_tests import StandardModelTest
+from .test_models.standard_model_tests import StandardModelTest
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index 8cfe1d0a..30433059 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -77,12 +77,16 @@ class TestDiscretise(unittest.TestCase):
         scal_disc = disc.process_symbol(scal)
         self.assertIsInstance(scal_disc, pybamm.Scalar)
         self.assertEqual(scal_disc.value, scal.value)
-
-        # parameter
-        par = pybamm.Parameter("par")
-        par_disc = disc.process_symbol(par)
-        self.assertIsInstance(par_disc, pybamm.Parameter)
-        self.assertEqual(par_disc.name, par.name)
+        # vector
+        vec = pybamm.Vector(np.array([1, 2, 3, 4]))
+        vec_disc = disc.process_symbol(vec)
+        self.assertIsInstance(vec_disc, pybamm.Vector)
+        np.testing.assert_array_equal(vec_disc.entries, vec.entries)
+        # matrix
+        mat = pybamm.Matrix(np.array([[1, 2, 3, 4], [5, 6, 7, 8]]))
+        mat_disc = disc.process_symbol(mat)
+        self.assertIsInstance(mat_disc, pybamm.Matrix)
+        np.testing.assert_array_equal(mat_disc.entries, mat.entries)
 
         # binary operator
         bin = var + scal
@@ -91,6 +95,12 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(bin_disc.children[0], pybamm.StateVector)
         self.assertIsInstance(bin_disc.children[1], pybamm.Scalar)
 
+        bin2 = scal + var
+        bin2_disc = disc.process_symbol(bin2, y_slices)
+        self.assertIsInstance(bin2_disc, pybamm.Addition)
+        self.assertIsInstance(bin2_disc.children[0], pybamm.Scalar)
+        self.assertIsInstance(bin2_disc.children[1], pybamm.StateVector)
+
         # non-spatial unary operator
         un1 = -var
         un1_disc = disc.process_symbol(un1, y_slices)
@@ -102,14 +112,18 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(un2_disc, pybamm.AbsoluteValue)
         self.assertIsInstance(un2_disc.children[0], pybamm.Scalar)
 
+        # parameter should fail
+        with self.assertRaises(NotImplementedError):
+            disc.process_symbol(pybamm.Parameter("par"))
+
     def test_process_complex_expression(self):
         var1 = pybamm.Variable("var1")
         var2 = pybamm.Variable("var2")
-        par1 = pybamm.Parameter("par1")
-        par2 = pybamm.Parameter("par2")
         scal1 = pybamm.Scalar("scal1")
         scal2 = pybamm.Scalar("scal2")
-        expression = (scal1 * (par1 + var2)) / ((var1 - par2) + scal2)
+        scal3 = pybamm.Scalar("scal3")
+        scal4 = pybamm.Scalar("scal4")
+        expression = (scal1 * (scal3 + var2)) / ((var1 - scal4) + scal2)
 
         disc = pybamm.BaseDiscretisation(None)
         y_slices = {var1.id: slice(53), var2.id: slice(53, 59)}
@@ -120,7 +134,7 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(exp_disc.children[0].children[0], pybamm.Scalar)
         self.assertIsInstance(exp_disc.children[0].children[1], pybamm.Addition)
         self.assertTrue(
-            isinstance(exp_disc.children[0].children[1].children[0], pybamm.Parameter)
+            isinstance(exp_disc.children[0].children[1].children[0], pybamm.Scalar)
         )
         self.assertTrue(
             isinstance(exp_disc.children[0].children[1].children[1], pybamm.StateVector)
@@ -140,7 +154,7 @@ class TestDiscretise(unittest.TestCase):
             exp_disc.children[1].children[0].children[0].y_slice, y_slices[var1.id]
         )
         self.assertTrue(
-            isinstance(exp_disc.children[1].children[0].children[1], pybamm.Parameter)
+            isinstance(exp_disc.children[1].children[0].children[1], pybamm.Scalar)
         )
         self.assertIsInstance(exp_disc.children[1].children[1], pybamm.Scalar)
 
diff --git a/tests/test_discretisations/test_finite_volume_discretisations.py b/tests/test_discretisations/test_finite_volume_discretisations.py
index c443bf30..2db1dd9f 100644
--- a/tests/test_discretisations/test_finite_volume_discretisations.py
+++ b/tests/test_discretisations/test_finite_volume_discretisations.py
@@ -57,6 +57,33 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             # Check that the equation can be evaluated
             eqn_disc.evaluate(None, y_test)
 
+        # more testing, with boundary conditions
+        for flux, eqn in zip(
+            [
+                pybamm.grad(var),
+                pybamm.grad(var),
+                pybamm.grad(var),
+                2 * pybamm.grad(var),
+                2 * pybamm.grad(var),
+                var * pybamm.grad(var) + 2 * pybamm.grad(var),
+            ],
+            [
+                pybamm.div(pybamm.grad(var)),
+                pybamm.div(pybamm.grad(var)) + 2,
+                pybamm.div(pybamm.grad(var)) + var,
+                pybamm.div(2 * pybamm.grad(var)),
+                pybamm.div(2 * pybamm.grad(var)) + 3 * var,
+                -2 * pybamm.div(var * pybamm.grad(var) + 2 * pybamm.grad(var)),
+            ],
+        ):
+            eqn_disc = disc.process_symbol(
+                eqn,
+                y_slices,
+                {flux.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}},
+            )
+            # Check that the equation can be evaluated
+            eqn_disc.evaluate(None, y_test)
+
     def test_add_ghost_nodes(self):
         # Set up
         param = pybamm.ParameterValues(