docs/source/examples/notebooks/creating_models/2-a-pde-model.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a simple PDE model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [previous notebook](./1-an-ode-model.ipynb) we show how to create, discretise and solve an ODE model in pybamm. In this notebook we show how to create and solve a PDE problem, which will require meshing of the spatial domain.\n",
    "\n",
    "As an example, we consider the problem of linear diffusion on a unit sphere,\n",
    "$$\n",
    "  \\frac{\\partial c}{\\partial t} = \\nabla \\cdot (\\nabla c),\n",
    "$$\n",
    "with the following boundary and initial conditions:\n",
    "$$\n",
    "  \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=0} = 0, \\quad \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=1} = 2, \\quad \\left.c\\right\\vert_{t=0} = 1.\n",
    "$$\n",
    "\n",
    "As before, we begin by importing the PyBaMM library into this notebook, along with any other packages we require:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import pybamm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the previous example, we start with a `pybamm.BaseModel` object and define our model variables. Since we are now solving a PDE we need to tell pybamm the domain each variable belongs to so that it can be discretised in space in the correct way. This is done by passing the keyword argument `domain`, and in this example we choose the domain \"negative particle\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.BaseModel()\n",
    "\n",
    "c = pybamm.Variable(\"Concentration\", domain=\"negative particle\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have given our variable the (useful) name \"Concentration\", but the symbol representing this variable is simply `c`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then state out governing equations. Sometime it is useful to define intermediate quantities in order to express the governing equations more easily. In this example we define the flux, then define the rhs to be minus the divergence of the flux. The equation is then added to the dictionary `model.rhs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = -pybamm.grad(c)  # define the flux\n",
    "dcdt = -pybamm.div(N)  # define the rhs equation\n",
    "\n",
    "model.rhs = {c: dcdt}  # add the equation to rhs dictionary"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike ODE models, PDE models require both initial and boundary conditions. Similar to initial conditions, boundary conditions can be added using the dictionary `model.boundary_conditions`. Boundary conditions for each variable are provided as a dictionary of the form `{side: (value, type)`, where, in 1D, side can be \"left\" or \"right\", value is the value of the boundary conditions, and type is the type of boundary condition (at present, this can be \"Dirichlet\" or \"Neumann\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial conditions\n",
    "model.initial_conditions = {c: pybamm.Scalar(1)}\n",
    "\n",
    "# boundary conditions\n",
    "lbc = pybamm.Scalar(0)\n",
    "rbc = pybamm.Scalar(2)\n",
    "model.boundary_conditions = {c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in our example the boundary conditions take constant values, but the value can be any valid pybamm expression.\n",
    "\n",
    "Finally, we add any variables of interest to the dictionary `model.variables`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.variables = {\"Concentration\": c, \"Flux\": N}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the model is now completely defined all that remains is to discretise and solve. Since this model is a PDE we need to define the geometry on which it will be solved, and choose how to mesh the geometry and discretise in space."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining a geometry and mesh\n",
    "\n",
    "We can define spatial variables in a similar way to how we defined model variables, providing a domain and a coordinate system. The geometry on which we wish to solve the model is defined using a nested dictionary. The first key is the domain name (here \"negative particle\") and the entry is a dictionary giving the limits of the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define geometry\n",
    "r = pybamm.SpatialVariable(\n",
    "    \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n",
    ")\n",
    "geometry = {\n",
    "    \"negative particle\": {r: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then create a uniform one-dimensional mesh with 20 points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mesh and discretise\n",
    "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n",
    "var_pts = {r: 20}\n",
    "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example of meshes that do require parameters include the `pybamm.Exponential1DSubMesh` which clusters points close to one or both boundaries using an exponential rule. It takes a parameter which sets how closely the points are clustered together, and also lets the users select the side on which more points should be clustered. For example, to create a mesh with more nodes clustered to the right (i.e. the surface in the particle problem), using a stretch factor of 2, we pass an instance of the exponential submesh class and a dictionary of parameters into the `MeshGenerator` class as follows: `pybamm.MeshGenerator(pybamm.Exponential1DSubMesh, submesh_params={\"side\": \"right\", \"stretch\": 2})`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining a mesh we choose a spatial method. Here we choose the Finite Volume Method. We then set up a discretisation by passing the mesh and spatial methods to the class `pybamm.Discretisation`. The model is then processed, turning the variables into (slices of) a statevector, spatial variables into vector and spatial operators into matrix-vector multiplications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n",
    "disc = pybamm.Discretisation(mesh, spatial_methods)\n",
    "disc.process_model(model);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the model has been discretised we are ready to solve. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we choose a solver and times at which we want the solution returned. We then solve, extract the variables we are interested in, and plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2021-11-19 15:31:50,774 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for negative particle. Using default of 1 [m].\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABQiUlEQVR4nO3dd3yV9fn/8deVRQKEGXYIewjIDOCeKO7RuveurbXaVm3t/nVacVstRdzbumqtioiCoLIFBCFhhD3DCBlknuv3R45+KWWcQM65T5L38/HIg5z73Oect7eQT677/tzXx9wdERERERERkWhLCDqAiIiIiIiINAwqQEVERERERCQmVICKiIiIiIhITKgAFRERERERkZhQASoiIiIiIiIxkRR0gN1lZGR4165dg44hIiIN2Jw5c/LdvU3QOQ6WxlIREYkH+xpP46oA7dq1K7Nnzw46hoiINGBmtiroDIdCY6mIiMSDfY2nmoIrIiIiIiIiMaECVERERERERGJCBaiIiIiIiIjEhApQERERERERiQkVoCIiIiIiIhITKkBFREREREQkJlSAioiIBMTMEs3sSzN7dz/7DDezKjO7IJbZREREokEFqIiISHBuAxbv60kzSwT+CkyIWSIREZEoUgEqIiJ1VlllFU9/lsfv3lkUdJQaM7NM4Exg/H52uxV4A9gck1AiItIgbSks4+qnZrJmW0nUPysp6p8gIiJSy6pCzptz1/LQR0tZt2MXR/VoTXlliJSkOnVe9SHgLiB9b0+aWSfgfOAkYPj+3sjMbgJuAsjKyqrVkCIiUr8VlFRw5ZMzWLW1hC1FZXRu1Tiqn6cCVERE6gx354OFG7l/Yi7LNhcxMLM593z3cI7pmYGZBR0vYmZ2FrDZ3eeY2Qn72O0h4GfuXnWg/zZ3HweMA8jOzvbaSyoiIvVZcVkl1z4zkxVbinnymmyGZrWM+meqABURkbjn7kxdms+YCTl8ta6Anm2bMvaKoYzu375OFZ67ORo4x8zOAFKBZmb2grtfsds+2cAr4f++DOAMM6t097djnlZEROqd0ooqvvf8HOat2cHjlw/j2F5tYvK5KkBFRCSuzVm1nTETljB9xTY6tUjjvgsHcf6QTiQm1MnCEwB3vxu4GyB8BfSOPYpP3L3bN9+b2TPAuyo+RUSkNlRWhbjtlS+Ztiyf+y4cxGkD2sfss1WAiohIXFqycSf3Tcjho8WbyWiawu/O7selI7NolJQYdLSoMbObAdx9bNBZRESkfgqFnLveWMCERZv43dn9uGBYZkw/XwWoiIjElVVbi3lgYi7vzF9P00ZJ3Dm6D9ce3ZXGKfVzyHL3ycDk8Pd7LTzd/ZrYJRIRkfrK3fn9u1/z5tx1/PSU3lxzdLcDv6iW1c/RXERE6pxNO0t5ZNJSXp21hqRE4+bje3DzcT1o3jg56GgiIiL1woMTc3nm85XceGw3fnhSz0AyqAAVEZFAbS8uZ+yU5Tzz+UqqQs6lI7K49aSetG2WGnQ0ERGRemP81BU88vEyLs7uzC/OOCywJn4qQEVEJBDFZZU8OS2PJz5dQVF5JecP7sTto3qT1Tq664+JiIg0NK/MXM0f/7OYMw/vwJ+/c3igHeRVgIqISEyVVVbx4vTVPPbJMrYWl3Nqv3b89NQ+9GmfHnQ0ERGReufdBeu5+62vOL53Gx68eHDgXeSjVoCaWR/g1d02dQd+4+4PReszRUQkflVWhXhz7joe+iiX9QWlHNWjNXeO7sOQGCx6LSIi0hB9krOZH786j+wuLRl7xTBSkhKCjhS9AtTdc4DBAGaWCKwD3orW54mISHwKhZwPFm3kvg9zWLGlmEGZzbn3gkEc0ysj6GgiIiL11sy8bdz8/Bz6tE/nyWuGk5YSH8uYxWoK7snAcndfFaPPExGRgLk7ny7NZ8yEJSxct5NebZsy9ophjO7fLtB7T0REROq7r9YWcN0zs8hsmcaz146gWWr8dJSPVQF6CfDy3p4ws5uAmwCysrJiFEdERKJpzqpt3PtBDjPytpHZMo37LxzEeUM6BX7fiYiISH23bHMhVz89k+Zpybxww0haN20UdKT/EvUC1MxSgHOAu/f2vLuPA8YBZGdne7TziIhI9CzesJP7JuQwaclmMpo24v+d059LRnSmUVJ8TPsRERGpz9ZsK+GK8TNJMOPFG0bSoXla0JH+RyyugJ4OzHX3TTH4LBERCcCqrcU8MDGXd+avJ71REneO7sO1R3elcYqarYuIiMTC5p2lXPHkDHZVVPHq946ga0aToCPtVSx+M7iUfUy/FRGRum1jQSmPfLyU12atITkxge8f34PvHdeD5o3j514TERGR+m5HSTlXPjmTLYVlvHDDSPq2bxZ0pH2KagFqZo2BU4DvRfNzREQktrYVl/P3yct47otVhNy5fGQWt5zUk7bpqUFHExERaVCKyiq5+ulZ5OUX8/S1wxka58ubRbUAdfcSoHU0P0NERGKnqKySJ6fm8cTUFZSUV3L+kExuH9WLzq0aBx1NRESkwSmtqOKm52azcF0Bf798KEf3jP8lznRzjoiIHFBpRRUvTF/F45OXs624nNH92/HTU/vQu1160NFEREQapIqqED986Us+X76VBy8exKn92wcdKSIqQEVEZJ8qq0K8PmctD09ayoaCUo7pmcGdo/swqHOLoKOJiIg0WKGQc+c/5/PR4k38/tz+nD8kM+hIEVMBKiIi/yMUct5buIEHPsxlRX4xgzu34P4LB3FUHZjaIyIiUp+5O799ZxFvz1vPHaf25qojuwYdqUZUgIqIyLfcncm5W7hvQg6L1u+kd7umjLtyGKf0a4eZBR1PRESkwbvvwxyen76K7x3XnVtO7Bl0nBpTASoiIgDMXrmNez/IYebKbXRulcYDFw3i3MGdSExQ4SkiIhIPxk5ZzmOfLOfSEZ35+el96+TJYRWgIiIN3KL1Bdw3IYdPcrbQJr0Rfzi3PxcPzyIlKSHoaCIiIhL20ozV3PP+Es4a2IE/nnd4nSw+QQWoiEiDlZdfzAMTc/n3/PU0T0vmZ6f15eqjutA4RUODiIhIPHln/np++fZXnNS3LQ9ePLhOz07SbxkiIg3MhoJdPDJpKa/NXktKYgK3nNiDm47rQfO05KCjiYiIyB4mLd7ET16dx/CurXj88qEkJ9btGUoqQEVEGohtxeU8/skynpu+CnfnyiO6cMuJPWmT3ijoaCIiIrIXXyzfyg9enMthHZrx5NXZpCYnBh3pkKkAFRGp54rKKhk/dQXjp+ZRUl7J+UMyuX1ULzq3ahx0NBEREdmHeWt2cONzs+ncqjHPXjeC9NT6MVNJBaiISD1VWlHFC9NX8fjk5WwrLuf0Ae35ySm96dUuPehoIiIish8L1xVw1ZMzaNUkhReuH0mrJilBR6o1KkBFROqZyqoQ/5yzlkcmLWVDQSnH9srgztF9GJjZIuhoIiIicgA5Gwu58skZpKcm89KNI2nfPDXoSLVKBaiISD0RCjnvfrWBByfmkpdfzJCsFtx/0SCO6pERdDQRERGJwPItRVw+fgbJiQm8eMNIMlvWv9tlVICKiNRx7s4nOZsZMyGXxRt20rd9OuOvyubkw9rW2TXCREREGpqV+cVc/sQMwHnpxiPpmtEk6EhRoQJURKQOm5m3jTETljBr5XayWjXmoYsHc/agjnV6fTAREZGGZtXWYi59YjpllVW8dOMR9GzbNOhIUaMCVESkDlq4roAxE3KYkruFtumN+ON5A7h4eOc6vzaYiIhIQ7N6awmXjpvOrooqXrrhCA7r0CzoSFGlAlREpA5ZvqWIBybm8p8FG2ielszPT+/L1Ud2JS2l7q8LJiIi0tCs2VbCpU9Mp7i8ihdvGEm/jvW7+AQVoCIidcL6Hbt4+KOlvD53LY2SErj1pJ7ceFx3mtWTNcFEREQamrXbq4vPwtIKXrrxCAZ0ah50pJhQASoiEse2FpXx2CfLeWH6KgCuOrILPzihJ23SGwWcTERERA7W+h27uPSJ6RTsquDFG0Y2mOITVICKiMSlwtIKnpiax5NTV7CroooLhmXyo5N71ct27CIiIg3JhoJdXDJuOjuKK3jhhpENbp1uFaAiInGktKKK579YxeOTl7G9pIIzDm/PT07pTc+26UFHExERkUO0saCUS8dNZ1txOc9fP4JBnVsEHSnmVICKiMSBiqoQ/5y9lkcmLWXjzlKO692GO0/tw+GZDWdKTkNkZonAbGCdu5+1x3OXAz8LPywCvu/u82McUUREasmmnaVc9sR08ovKefa6EQzJahl0pECoABURCVAo5Px7wXoenJjLyq0lDM1qwYMXD+bIHq2DjiaxcRuwGNhb28M84Hh3325mpwPjgJGxDCciIrVj885SLn1iOpt2lvLc9SMY1qVhFp+gAlREJBDuzsdLNjNmQg5LNhbSt306T16dzUl922JmQceTGDCzTOBM4E/AT/Z83t0/3+3hdCAzRtFERKQWbS6sLj43FpTy7HUjGNalVdCRAqUCVEQkxqav2MqYCTnMWbWdLq0b8/Algzl7YEcSElR4NjAPAXcBkdzgez3w/r6eNLObgJsAsrKyaiObiIjUgi2FZVz+xAzW7yjlmWuHM7xrwy4+IcoFqJm1AMYDAwAHrnP3L6L5mSIi8eqrtQWM+TCHT3O30K5ZI/58/uFcmJ1JcmJC0NEkxszsLGCzu88xsxMOsO+JVBegx+xrH3cfR/UUXbKzs732koqIyMHKLyrj8vHTWbO9hKevGcHI7rq9BqJ/BfRh4AN3v8DMUgCtHyAiDc6yzUU8MDGH977aSMvGyfzijL5cdWRXUpMTg44mwTkaOMfMzgBSgWZm9oK7X7H7TmY2kOoTuae7+9YAcoqIyEHYVlzOFeNnsHpbCU9dPVy9HXYTtQLUzJoBxwHXALh7OVAerc8TEYk363bs4qGJubwxdy1pyYn86ORe3HhsN9JTk4OOJgFz97uBuwHCV0Dv2EvxmQW8CVzp7rmxzigiIgdne3E5lz0xnbz8Yp66ZjhH9cwIOlJcieYV0O7AFuBpMxsEzAFuc/fi3XfSfSsiUt/kF5Xx2CfLeHH6ajC45qhu3HJiD1o3bRR0NIlzZnYzgLuPBX4DtAYeDzemqnT37ADjiYjIAewoKefy8TPIyy/myauHc7SKz/9h7tG5VcTMsqnu2ne0u88ws4eBne7+6329Jjs722fPnh2VPCIi0baztIInPl3Bk9PyKK2o4sJhnbltVC86tkgLOprUgJnNqcuFnsZSEZFgbC0q44onZ7J8SxHjr8rmuN5tgo4UqH2Np9G8AroWWOvuM8KPXwd+HsXPExEJxK7yKp77YiV/n7KcHSUVnHl4B35yam96tGkadDQRERGJgc2FpVwxfgartpbw5NXZHNurYRef+xO1AtTdN5rZGjPr4+45wMnA19H6PBGRWKuoCvHqrDU8MmkpmwvLOL53G+4c3YcBnZoHHU1ERERiZGNBKZeNn86GHaU8fe1wjuqhabf7E+0uuLcCL4Y74K4Aro3y54mIRF0o5Lwzfz0PTMxl9bYSsru05NFLh6i9umBmrdx9W9A5REQkNtbt2MVlT0wnv7CM564foXU+IxDVAtTd5wF19j4aEZHduTuTFm/mvg9zWLKxkMM6NOOpa7I5sU9bwk1ipAExs1+5+x/D3/cD3gaSrfovw8W73YIiIiL1UF5+MVeMn8HO0gqev2EkQ7NaBh2pToj2FVARkXrhi+VbGTNhCXNX76Br68Y8fMlgzh7YkYQEFZ4N2HeAP4a/H0N1p/f3zWwE8BBwVFDBREQkunI3FXL5+BlUVoV4+cYjdPtNDRywADWzRsB3ga677+/uv49eLBGR+PDV2gLunbCEqUvzad8slb9853AuGJZJcmJC0NEkvnR09/cB3H2mman1sYhIPfXV2gKuemoGyYkJvPa9I+nVLj3oSHVKJFdA/wUUUL2OZ1l044iIxIdlmwu5/8Nc3l+4kZaNk/nlGYdx5ZFdSE1ODDqaxI/uZvYOYECmmTV295Lwc8kB5hIRkSiZs2ob1zw1i2Zpybx4w0i6ZjQJOlKdE0kBmunup0U9iYhIHFi7vYSHP1rKG3PXkpacyG0n9+KGY7uRnqp6Qv7HuXs8TgAws3bA32MfR0REoumzZfnc8Oxs2jdP5YUbRtJJ63wflEgK0M/N7HB3/yrqaUREArKlsIzHPlnGSzNWg8F1R3fjByf2pFWTlKCjSZxy9yn72L4JeCzGcUREJIomLd7E91+cS7fWTXj+hhG0TU8NOlKdFUkBegxwjZnlUT0F1wB394FRTSYiEgMFuyp44tMVPPVZHmWVIS7KzuRHJ/eiQ3Od1ZSDZ2Y3ufu4oHOIiMih+8+CDdz2ypcc1qEZz103gpY6OX1IIilAT496ChGRGNtVXsUzn69k7JTlFOyq4KyBHfjJKb3p3qZp0NGkflB7ZBGReuD1OWu56/X5DOvSkievGU4z3ZJzyA5YgLr7KjMbBBwb3jTV3edHN5aISHSUV4Z4ddZqHvl4GVsKyzixTxvuGN2H/h3VPl1qj7v/I+gMIiJyaJ7/YiW//tcijumZwbirhtE4RStY1oZIlmG5DbgReDO86QUzG+fuj0Y1mYhILaoKOf+at44HP8plzbZdDO/akscvH8rwrq2CjiZ1mJn1pboZUSfAgfXAO+6+ONBgIiJySP4xZTl/eX8Jow5ry98uG6ou+LUokjL+emCkuxcDmNlfgS8AFaAiEvfcnYlfb+L+D3PJ2VRIvw7NePraAZzQuw1mmiUpB8/MfgZcCrwCzAxvzgReNrNX3P2ewMKJiMhBcXce+mgpD09aylkDO/DgxYO19ncti6QANaBqt8dV6N4WEakDPl+Wz70Tcpi3ZgfdMprw6KVDOPPwDiQk6EeY1Irrgf7uXrH7RjN7AFgEqAAVEalD3J0//Wcx46flceGwTO757kAS9TtDrYukAH0amGFmb4Ufnwc8GbVEIiKHaP6aHYyZkMO0Zfl0aJ7KPd85nAuGZZKkM5hSu0JAR2DVHts7hJ8TEZE6IhRyfvWvhbw0YzVXH9mF357dXyesoySSJkQPmNlkqpdjMeBad/8y2sFERGpq6aZC7vswhwmLNtGqSQq/OvMwrjiii+7bkGi5HZhkZkuBNeFtWUBP4IdBhRIRkZqprApx1+sLePPLdXz/hB7cNbqPbtOJon0WoGbWzN13mlkrYGX465vnWrn7tujHExE5sDXbSnjoo6W89eVaGqck8eNRvbn+2G40baRudRI97v6BmfUGRlDdhMiAtcAsd6/a74tFRCQulFZUcdsrXzJh0SbuOLU3PzypV9CR6r39/Xb2EnAWMIfqzn7fsPDj7lHMJSJyQFsKy/jbx0t5aeZqzIzrj+nG90/oSSstEC0x4u4hYHrQOUREpOYKSyu46bk5fLFiK785qx/XHdMt6EgNwj4LUHc/K/yn/k+ISFwp2FXBuE+X89S0lZRXhbgouzM/OrknHZqnBR1NGigzm+bux3zzZ9B5RERk/7YUlnHN0zPJ2VjIQxcP5rwhnYKO1GBEsg7oJHc/+UDbRESiraS8kmc+X8nYycvZWVrJ2YM68pNTetMto0nQ0UQah//UX0YRkTi3ZlsJVz45g407S3ni6mxO7NM26EgNyv7uAU2lekDNMLOW/N/SK82o7vonIhIT5ZUhXpm1mkc/XsaWwjJO7tuWn57ah34dmwUdTUREROqQResLuObpWZRXhnjxhiMY1qVl0JEanP1dAf0e1R3+OlJ9H+g3BehO4LHoxhIRgaqQ8/aX63jwo1zWbt/FiG6t+PvlQ8nu2iroaCIiIlLHfL4sn5uen0Oz1CReuvlIerVLDzpSg7S/e0AfBh42s1vd/dEYZhKRBs7dmbBoE/d/mMPSzUUM6NSMP51/OMf1ylBbdBEREamxf89fz09em0f3jKY8c91w9Y0IUCTrgD5qZgOAfkDqbtufi2YwEWmYpi3NZ8yEJcxfW0D3Nk147LKhnD6gvRaDlninv6AiInHqyWl5/PE/XzO8ayueuCqb5mnJQUdq0CJpQvRb4ASqC9D3gNOBaYAKUBGpNV+u3s6YCTl8vnwrHZuncu93B/KdoZ1ISkwIOppIJH68x58iIhKwUMj56wdL+MenKzitf3seumQwqcmJQcdq8CJZpf0CYBDwpbtfa2btgPHRjSUiDUXOxkLu+zCHiV9vonWTFH5zVj8uG5mlAULqDDO70N3/CeDuk/fcJiIisVdeGeKu1+fz9rz1XHVkF357dn8SNZsqLkRSgO5y95CZVZpZM2Az0D3KuUSknluzrYQHJ+by1rx1NE1J4qen9ObaY7rRtFEkP5ZE4srdwJ7F5t62iYhIDBSVVfL9F+YwdWk+d47uww9O6KEeEnEkkt/0ZptZC+AJqrvhFgEzI3lzM1sJFAJVQKW7Zx9cTBGpLzbvLOVvnyzj5ZmrSTDjpmO7c/PxPWjZJCXoaCI1YmanA2cAnczskd2eagZUBpNKRKRh21JYxrXPzGTxhkLuvWAgF2V3DjqS7GG/BahVnyr4i7vvAMaa2QdAM3dfUIPPONHd8w8ho4jUAwUlFYz9dDlPf5ZHZZVz8fDO3HpSL9o3Tz3wi0Xi03pgNnAO1Sdov1GI7gUVEYm5lfnFXPXUTDYXlvLEVcM4qW+7oCPJXuy3AHV3N7O3gWHhxytjkElE6pGS8kqe/mwlY6csp6isknMGdeTHo3rTNaNJ0NFEDom7zwfmm9lL7l4RdB4RkYZswdodXPv0LELuvHzjEQzJahl0JNmHSKbgTjez4e4+6yDe34EPzcyBf7j7uD13MLObgJsAsrKyDuIjRCQelVVW8crMNTz68TLyi8oYdVhbfnpqHw7r0CzoaCK1SsWniEiwpuRu4fsvzKFVkxSevW4EPdo0DTqS7EckBeiJwPfMbBVQTPVaZ+7uAyN47dHuvt7M2gITzWyJu3+6+w7honQcQHZ2ttcsvojEm6qQ89aX63hwYi7rduxiZLdW/OPKoQzr0iroaCIiIlLPvDl3LXe9voBe7dJ59trhtG2mW3viXSQF6OkH++buvj7852YzewsYAXy6/1eJSF3k7kxYtIn7P8xh6eYiDu/UnL9853CO7ZWhznNSr+1tyZVIl2Exs0Sq7yNd5+5n7fGcAQ9T3eioBLjG3efWXnIRkbrL3fnHpyu45/0lHNWjNWOvHEaz1OSgY0kEIilA/+juV+6+wcyeB67cx/7f7NMESHD3wvD3pwK/P+ikIhKX3J1py/IZMyGHBWsL6N6mCY9fPpTTB7RX4SkNxaEsw3IbsJjqzrl7Oh3oFf4aCfw9/KeISINWFXL+8O7XPPP5Ss4e1JH7LhxIoyStH15XRFKA9t/9Qfhs7bAIXtcOeCv8C2gS8JK7f1DjhCISt+au3s6YD3L4YsVWOrVI494LBvKdIZ1ISkwIOppI1B3qMixmlgmcCfwJ+MledjkXeM7dnep+DC3MrIO7bzj09CIidVNJeSU/enkeHy3exPXHdOOXZxxGQoJOeNcl+yxAzexu4BdAmpnt/GYzUE74ns39cfcVwKDaCCki8WXJxp3cNyGXjxZvIqNpCr89ux+XjczS2UdpaA51GZaHgLuA9H083wlYs9vjteFtKkBFpEHatLOU65+dxdfrd/K7s/txzdHdgo4kB2GfBai7/wX4i5n9xd3vjmEmEYlTq7YW8+DEXP41fz1NGyVxx6m9ufbobjRpFMlkCpH65VCWYTGzs4DN7j7HzE7Y1257+9h9vJ86yotIvbZ4w06uf2YWO3ZVMP7qbK3xWYcd8LdGd7/bzDoBXXbff89utiJSf23aWcojk5by6qw1JCUa3zuuBzcf350WjVOCjiYSD7qa2V+AfsC37Rfdvft+XnM0cI6ZnRF+TTMze8Hdr9htn7VA590eZ1J91fV/qKO8iNRnHy/ZxK0vfUl6ajL/vPlI+ndsHnQkOQQHLEDN7B7gEuBroCq82VE3W5F6b0dJOX+fspxnP19JZZVzyYjO3HpSL9qpxbnI7p4Gfgs8SPXSZdey96uX3wrPLLobIHwF9I49ik+Ad4AfmtkrVDcfKtD9nyLSkLg7T07L48/vLaZ/x+Y8cVU27Zvrd5C6LpJ5c+cDfdy9LNphRCQ+FJdV8tS0PMZ9uoKi8krOG9yJH4/qTVbrxkFHE4lHae4+yczM3VcBvzOzqVQXpTViZjcDuPtY4D2qmxwto3oZlmtrMbOISFyrqArxm38t5OWZazh9QHseuGgwaSnqNVEfRFKArgCSARWgIvVcWWUVL81YzWOfLCO/qJxT+rXjp6f2pm/7va0QISJhpWaWACw1sx8C64C2kb7Y3ScDk8Pfj91tuwO31GpSEZE6YEdJOd9/YS5frNjKLSf24Ken9FGn23okkgK0BJhnZpPYrQh19x9FLZWIxFRlVYg3v1zHwx8tZd2OXRzZvTXjrurD0KyWQUcTqQtuBxoDPwL+QPU03KuDDCQiUlet2FLE9c/OZt32Xdx/4SC+Oywz6EhSyyIpQN8Jf4lIPePufLBwI/d9mMPyLcUMymzOX787kKN7tia8hq+IHIC7zwp/W4SmyYqIHLTPl+Xz/RfnkphgvHjjSIZ3bRV0JImCSLrgPmtmaUCWu+fEIJOIRJm7M3VpPmMm5PDVugJ6tm3K2CuGMrp/exWeIiIiEnMvz1zNr99eSLeMJjx1zXA6t1Lfifoqki64ZwP3ASlANzMbDPze3c+JcjYRiYI5q7Zz7wdLmJG3jU4t0rjvwkGcP6QTibq3QkRERGKsKuT8+b3FPDktj+N7t+HRy4bQLDU56FgSRZFMwf0dMIL/a5Awz8y6RTGTiETB4g07uf/DHD5avJmMpo34f+f055IRnWmUpI5yIiIiEntFZZX86OUv+XjJZq45qiu/OvMwkhITgo4lURZJAVrp7gV7TMvTItcidcSqrcU8MDGXd+avp2mjJO4c3YdrjupKk0aR/PMXkQMxszbAjUBXdhtX3f26oDKJiMS7tdtLuOHZ2SzdXMQfzu3PlUd2DTqSxEgkv4EuNLPLgEQz60V1l7/PoxtLRA7Vpp2lPDxpKa/NWkNSonHz8T343nHdadE4JehoIvXNv4CpwEdAVcBZRETi3tzV27npudmUVYZ45trhHNurTdCRJIYiKUBvBX5J9RIsLwETgD9GM5SIHLztxeWMnbKcZz5fSVXIuXREFree1JO2zVKDjiZSXzV2958FHUJEpC7417x13Pn6Ato3S+WVm7Lp2TY96EgSY5F0wS2hugD9ZfTjiMjBKi6r5MlpeTzx6QqKyis5f3Anbh/Vm6zW6iInEmXvmtkZ7v5e0EFEROKVu/PQR0t5eNJSRnRrxdgrhtGqiWZlNUSRdMGdCFzo7jvCj1sCr7j76ChnE5EIlFVW8eL01Tz2yTK2Fpdzar92/PTUPvRprzOKIjFyG/ALMysHKsLb3N2bBZhJRCRuFJVVcsdr8/lg0UYuHJbJn84/nJQkNRtqqCKZgpvxTfEJ4O7bzaxt9CKJSCQqq0K8OXcdD32Uy/qCUo7u2Zo7R/dlcOcWQUcTaVDcXWd7RET2YcWWIm56fg55+cX86szDuP6YblpzvIGLpAANmVmWu68GMLMuqAuuSGBCIef9hRu5f2IOK7YUM6hzC8ZcOIije2YEHU2kwTKzc4Djwg8nu/u7QeYREYkHkxZv4vZX5pGclMDz14/gqB76XUUiK0B/CUwzsynhx8cBN0UvkojsjbszJXcL932Yw8J1O+ndrin/uHIYp/ZrpzOJIgEys3uA4cCL4U23mdkx7v7zAGOJiAQmFHIenlR9v+eATs34x5XZdGqRFnQsiRORNCH6wMyGAkcABvzY3fOjnkxEvjV75TbunZDDzLxtZLZM4/4LB3HekE4kJqjwFIkDZwCD3T0EYGbPAl8CKkBFpMEpKKng9le/5JOcLXx3aCZ/On8AqcmJQceSOBLpSvSNgG3h/fuZGe7+afRiiQjA1+t3ct+HOXy8ZDNt0hvxh3P7c/HwLN24LxJ/WlA9TgI0DzCHiEhgvl6/k5tfmMOGgl388bwBXD4yS7O05H9E0gX3r8DFwCIgFN7sgApQkShZmV/MAxNzeWf+epqlJnHXaX245qiuNE6J9JyRiMTQX4AvzewTqmcKHQfcHWwkEZHYeuvLtdz95le0SEvh1e8dydCslkFHkjgVyW+z5wF93L0syllEGrwNBbt4ZNIyXpu9hpTEBH5wQg++d1wPmjdODjqaiOyDu79sZpOpvg/UgJ+5+8ZgU4mIxEZ5ZYg//edrnv1iFSO7teJvlw2lTXqjoGNJHIukAF0BJAMqQEWiZFtxOX+fvIxnv1iFu3PFyCxuOaknbdNTg44mIvtgZn3dfUm4TwLA2vCfHc2so7vPDSqbiEgsbNpZyg9enMucVdu54Zhu/Pz0viQl6jYh2b9ICtASYJ6ZTWK3ItTdfxS1VCINRFFZJU9OzeOJqSsoKa/k/CGZ3D6qF51bNQ46mogc2E+o7gp//16ec+Ck2MYREYmdmXnbuOWluRSXVfLopUM4e1DHoCNJHRFJAfpO+EtEaklpRRUvTF/F45OXs624nNH92/HTU/vQu53WsxepK9z9myXJTnf30t2fMzNNXxCResndefqzlfz5vcV0btWYF64fSZ/2+v1FIhfJMizPmlkK0Du8KcfdKyL9ADNLBGYD69z9rIOLKVI/VFaFeH3OWh6etJQNBaUc0zODO0f3YVDnFkFHE5GD9zkwNIJtIiJ1Wkl5JXe/+RX/mreeU/q14/6LBtEsVX0qpGYi6YJ7AvAssJLq5gqdzezqGizDchuwGGh2cBFF6r5QyHlv4QYe+DCXFfnFDOrcgvsvHMRRPTOCjiYiB8nM2gOdgDQzG0L1GAnV453m0YtIvbIyv5ibX5hDzqZC7hzdh+8f34MErUcuByGSKbj3A6e6ew6AmfUGXgaGHeiFZpYJnAn8iep7ZUQaFHdncu4W7puQw6L1O+ndrinjrhzGKf3aaV0skbpvNHANkAk8sNv2QuAXQQQSEYmGSYs3cfur80hMMJ65dgTH924TdCSpwyIpQJO/KT4B3D3XzCK91v4QcBewz4nhZnYT1U0cyMrKivBtReLf7JXbuPeDHGau3EbnVmk8cNEgzh3ciUSdLRSpF9z9WeBZM/uuu78RdB4RkdoWCjkPTVrKI5OW0r9jM8ZeMUyNEuWQRVKAzjazJ4Hnw48vB+Yc6EVmdhaw2d3nhKfx7pW7jwPGAWRnZ3sEeUTi2qL1Bdw3IYdPcrbQJr0Rfzi3PxcPzyIlSW3JReojd3/DzM4E+gOpu23/fXCpREQOzZbCMn7y2jymLs3ngmGZ/PG8AaQmJwYdS+qBSArQ7wO3AD+i+v6WT4HHI3jd0cA5ZnYG1QNyMzN7wd2vONiwIvEsL7+YBybm8u/562melsxdp/XhmqO60jglkn9mIlJXmdlYqu/5PBEYD1wAzAw0lIjIIfhsWT63vzqPnbsq+PP5h3PpiM66dUhqTSS/GScBD7v7A/BtV9tGB3qRu98N3B1+zQnAHSo+pT7aULCLRyYt5bXZa0lJTOCHJ/bkxuO60zxNXeFEGoij3H2gmS1w9/9nZvcDbwYdSkSkpiqrQjz00VIem7yMHm2a8vz1I+jbXn1EpXZFUoBOAkYBReHHacCHwFHRCiVSF2wrLufxT5bx3PRVuDtXHtGFW07sSZv0A56fEZH65Zs1QEvMrCOwFegWYB4RkRpbu72E21+Zx+xV27koO5PfndNfs7gkKiL5W5Xq7t8Un7h7kZnV6O5jd58MTK5ZNJH4VFRWyfipKxg/NY+S8krOH5LJ7aN66aZ8kYbr32bWAhgDzAUceCLQRCIiNfD+Vxv42RsLCDk8fMlgzh3cKehIUo9FUoAWm9lQd58LYGbDgF3RjSUSf0orqnhh+ioen7ycbcXlnNa/PT89tTe92u2zybOI1HNmlgBMcvcdwBtm9i7VJ24Lgk0mInJgpRVV/P7dr3lpxmoGZTbnkUuH0KV1k6BjST0XSQF6O/BPM1sfftwBuDhqiUTiTGVViNfnrOXhSUvZUFDKsb0yuOPUPgzq3CLoaCISMHcPhe/5PDL8uAwoCzaViMiB5Wws5NaX55K7qYjvHdedn57aRx37JSYOWIC6+ywz6wv0oboL7hJ3r4h6MpGAhULOf77awAMTc8nLL2ZIVgvuv2gQR/XICDqaiMSXD83su8Cb7q7lxEQkrrk7L85YzR/e/Zr01CSeu24Ex/VuE3QsaUAiurM4XHAujHIWkbjg7kzO2cKYCTl8vWEnfdunM/6qbE4+rK1akIvI3vwEaAJUmlkp1Sdr3d3VOlJE4kpBSQU/f3MB7y/cyLG9MnjgosFqnigxp9ZWIruZtXIb936whFkrt5PVqjEPXTyYcwZ1JCFBhaeI7J2760ZwEYl7s1Zu47aXv2RzYRl3n96XG4/trt9vJBAqQEWAResLuG9CDp/kbKFteiP+eN4ALh7emeRE3QshIvtnZpPc/eQDbdvL61KBT6leWzsJeN3df7vHPs2BF4Cs8D73ufvTtZlfROq3qpDz2CfLeOijXDJbNub17x/FYPWxkAAdsAC16jmHlwPd3f33ZpYFtHf3mVFPJxJlK7YU8cDEXN5dsIHmacn8/PS+XH1kV9JSEoOOJiJxLlxANgYyzKwl1VNvAZoBHSN4izLgpPDyZsnANDN7392n77bPLcDX7n62mbUBcszsRXcvr8X/FBGppzYWlHL7q18yfcU2zh3ckT+eN4D01OSgY0kDF8kV0MeBEHAS8HugEHgDGB7FXCJRtX7HLh79eCmvzV5Lo6QEbj2pJzce151m+qEsIpH7HtWd4jsCc/i/AnQn8NiBXhxuWPTNOtvJ4a89mxg5kB4+GdwU2AZUHmpwEan/Jn69iTtfn095ZYj7LhzEd4d2Ui8LiQuRFKAj3X2omX0J4O7bzSwlyrlEomJrURmPT17O89NXgcOVR3ThlhN76gZ8Eakxd38YeNjMbnX3Rw/mPcwskeritSfwmLvP2GOXvwHvAOuBdOBidw/t5X1uAm4CyMrKOpgoIlJPlFZUcc/7S3jm85X079iMRy4dQo82TYOOJfKtSArQivAA6QDhKUD/M/iJxLPC0grGT83jyWl5lJRX8t2hmdw2qheZLRsHHU1E6jh3f9TMjgK6stu46u7PRfDaKmCwmbUA3jKzAe6+e9f50cA8qmch9QAmmtlUd9+5x/uMA8YBZGdnaykYkQbq6/U7+clr81iysZBrj+7Kz0/vS6Mk3VYk8SWSAvQR4C2grZn9CbgA+FVUU4nUktKKKp7/YhWPT17G9pIKTuvfnjtG96ZnWzWtFJHaYWbPU10czgOqwpsdOGAB+g1332Fmk4HT+O9lz64F7glP111mZnlAX0B9GETkW5VVIcZOWc7Dk5bSonEKT12TzUl92wUdS2SvDliAuvuLZjYHOJnq+1vOc/fFUU8mcggqqkK8Pmctj0xayoaCUo7tlcGdo/swMLNF0NFEpP7JBvqFi8SIhWcUVYSLzzRgFPDXPXZbTfX4O9XM2gF9gBW1kFlE6ollm4v46WvzmL+2gLMGduAP5w6gZRPdLSfxK5IuuEcAi9z9sfDjdDMbuZf7VEQCFwo57361gQcn5pKXX8yQrBbcf9EgjuqREXQ0Eam/FgLtgQ01fF0H4NnwbS4JwGvu/q6Z3Qzg7mOBPwDPmNlXVJ8E/pm759dedBGpq6pCzlPT8hjzYQ6NUxJ59NIhnD0okgbcIsGKZAru34Ghuz0u3ss2kUC5O5NztjBmQg5fb9hJn3bpPHFVNqMOa6uObyISbRnA12Y2k+qlVQBw93P29yJ3XwAM2cv2sbt9vx44tfaiikh9sDK/mDv+OZ/Zq7Yz6rB2/Pk7A2ibnhp0LJGIRFKA2u7Titw9ZGaRvE4kJmbmbWPMhCXMWrmdrFaNeejiwZw9qCOJCSo8RSQmfhd0ABFpGEIh57kvVnLPB0tISUzggYsGcf4QLa8idUskheQKM/sR1Vc9AX6A7j+ROLBofQH3Tcjhk5wttE1vxB/PG8DFwzuTnJgQdDQRaUDcfYqZdQF6uftHZtYYUNtJEalVq7eWcOfr85mRt40T+7ThL98ZSPvmuuopdU8kBejNVHfC/RXVXf0mEV5rTCQIK7YU8cDEXN5dsIHmacn8/PS+XH1kV9JS9PueiMSemd1I9bjYiupuuJ2AsVQ3DxIROSTuzoszVvPn9xaTaMa93x3IhdmZuuopdVYkXXA3A5fEIIvIfm0o2MXDHy3ln3PW0igpgR+e2JMbj+tO87TkoKOJSMN2CzACmAHg7kvNrG2wkUSkPli3Yxc/e30B05blc2yvDO757kA6tUgLOpbIIYmkC24qcD3QH/j2Or+7XxfFXCLf2lZczuOfLOO56avA4cojunDLiT1pk94o6GgiIgBl7l7+zdWIcJ+EGi3JIiKyO3fntdlr+MO7i3F3/nT+AC4bkaWrnlIvRDIF93lgCTAa+D1wOaB1QCXqisoqGT91BeOn5lFSXsl3hmZy+6heZLZsHHQ0EZHdTTGzXwBpZnYK1b0S/h1wJhGpozYWlPLzNxcwOWcLR3Zvzb0XDKRzK/3uI/VHJAVoT3e/0MzOdfdnzewlYEK0g0nDVVpRxQvTV/H45OVsKy7ntP7tuWN0b3q2TQ86mojI3vyc6plCXwHfA94DxgeaSETqHHfnzbnr+N2/F1FZ5fz+3P5cMbILCerqL/VMJAVoRfjPHWY2ANgIdI1aImmwKqtCvD5nLQ9PWsqGglKO7ZXBnaP7MDCzRdDRRET2Jw14yt2fADCzxPC2kkBTiUidsXlnKb946ys+WryZEV1bMebCgXRp3SToWCJREUkBOs7MWlLdBfcdoCnw66imkgYlFHLeW7iBBz7MZUV+MYM7t+D+CwdxVM+MoKOJiERiEjAKKAo/TgM+BI4KLJGI1AmhUPW9nn95fwmlFVX8+qx+XHtUV131lHptnwWomd3m7g8Di919O/Ap0D1myaTec3em5G5hzIQcFq3fSe92TRl35TBO6ddON9mLSF2S6u7fFJ+4e1F4LVARkX3K3VTIL978itmrtjOiWyv+8p3D6dGmadCxRKJuf1dArwUeBh4Fhtb0jcPdcz8FGoU/53V3/+3BhJT6Z86qbfz1gxxm5m0js2UaD1w0iHMHdyJRZ/xEpO4pNrOh7j4XwMyGAbsCziQicWpXeRWPfLyUJz5dQXpqEmMuGMgFw7SupzQc+ytAF5vZSqCtmS3YbbsB7u4DD/DeZcBJ4TPBycA0M3vf3acfWmSpy5Zs3Ml9E3L4aPFmMpo24vfn9ueS4VmkJCUEHU1E5GDdDvzTzNaHH3cALg4ujojEq4+XbOI3/1rE2u27uGBYJr844zBaNUkJOpZITO2zAHX3S82sPdUdb8+p6Ru7u/N/98Mkh7+0LloDtWprMQ9MzOWd+etJb5TEnaP7cO3RXWmcEsltyCIi8cvdZ5lZX6AP1Sdpl7h7xQFeJiINyIaCXfy/d77mg0Ub6dW2Ka/edAQju7cOOpZIIA702/8W4Ct3X3Uwbx7uBDgH6Ak85u4z9rLPTcBNAFlZWQfzMRLHNu0s5ZFJS3l11hqSEo2bj+/Bzcf1oHnj5KCjiYjUpuFUd4hPAoaYGe7+XLCRRCRolVUhnvl8JQ9OzKXKnbtO68MNx3TXzC9p0PZbgLp7lZllmFmKu5fX9M3dvQoYbGYtgLfMbIC7L9xjn3HAOIDs7GxdIa0nCkoq+PuU5TzzeR6VVc4lIzrzo5N60bZZatDRRERqlZk9D/QA5gFV4c0OqAAVacDmrt7OL99ayOINOzmxTxt+f+4AOrdSfzKRSOY/rgI+M7N3gOJvNrr7A5F+iLvvMLPJwGnAwgPsLnVYSXklT3+2krFTllNUVsm5gzry41N6ay0rEanPsoF+4VtPRKSBKyip4K8TlvDyzNW0S09l7BVDGd2/vZoMiYRFUoCuD38lAOmRvrGZtQEqwsVnGtVrpP31oFJK3CuvDPHyzNU8+vEy8ovKGHVYW+4Y3Ye+7ZsFHU1EJNoWAu2BDUEHEZHguDtvz1vHn/6zmG3F5Vx7VDd+cmpvmjZSvwuR3R3wX4S7/7+DfO8OwLPh+0ATgNfc/d2DfC+JU1Uh51/z1vHgR7ms2baLEd1a8Y8rhzKsS6ugo4mIxEoG8LWZzaS6AzwA7l7jBn4iUjct31LEr99eyOfLtzKocwueuXYEAzo1DzqWSFw6YAFqZp+wl+617n7S/l7n7guAIQcfTeKZu/PR4s3cNyGHnE2F9OvQjGeuHcDxvdtoiomINDS/CzqAiASjtKKKxycvZ+zk5TRKTuCP5w3g0hFZWtdcZD8imRNwx27fpwLfBSqjE0fqgukrtnLvB0uYu3oH3TKa8OilQzjz8A4k6IetiDRA7j7FzNpR3QkXYKa7bw4yk4hE36e5W/j1vxayamsJ5w7uyK/O7Eeb9EZBxxKJe5FMwZ2zx6bPzGxKlPJIHFu4roAxE3KYkruF9s1S+ct3DueCYZkkJ6qVuIg0XGZ2ETAGmEz1OqCPmtmd7v56oMFEJCry8ov5y3uL+fDrTXTLaMKLN4zk6J4ZQccSqTMimYK7+818CcAwqpstSAORl1/MAxNz+ff89bRonMwvzujLVUd2JTU5MehoIiLx4JfA8G+ueoab8H0EqAAVqUd2lJTzyKRlPPfFSholJXDn6D5cf0w3/T4kUkORTMGdQ/U9oEb11Ns84PpohpL4sGlnKQ9PWsqrs9aQkpjArSf15MbjutMsNTnoaCIi8SRhjym3W6k+YSsi9UB5ZYjnp6/ikUlLKSyt4OLhnfnxKb1pm661zUUORiRTcLvFIojEj4KSCv4+ZTlPf5ZHyJ0rRmbxw5N66b4GEZG9+8DMJgAvhx9fDLwfYB4RqQXuzoRFm7jn/cWs3FrCMT0z+OWZh3FYBy0xJ3Io9lmAmtlwYI27bww/vorqBkSrgN+5+7bYRJRYKSmv5OnPVjJ2ynKKyio5f3AnfnxKbzq3ahx0NBGRuOXud5rZd4BjqJ4tNM7d3wo4logcgq/WFvCH/3zNzLxt9GzblKevHc4J6vQvUiv2dwX0H8AoADM7DrgHuBUYDIwDLoh2OImNiqoQr8xawyOTlrKlsIxRh7XljtF96NteZ/hERPbFzHoC7dz9M3d/E3gzvP04M+vh7suDTSgiNbWhYBdjPsjhzS/X0bpJCn84bwCXDu9MkhouitSa/RWgibtd5byY6jO6bwBvmNm8qCeTqAuFnHe/2sD9H+awamsJw7u25O+XDyW7a6sDv1hERB4CfrGX7SXh586OZRgROXjFZZX8Y8pyxk1dQcjh+yf04Psn9FDfC5Eo2G8BamZJ7l4JnAzcFOHrJM65O58uzefeD5awaP1O+rZP56lrsjmxT1tNLRERiVxXd1+w50Z3n21mXQPIIyI1VBVyXp+zhvs+zGVLYRnnDOrInaP76PYjkSjaXyH5MjDFzPKBXcBU+HbKUUEMskkUfLl6O/d+kMMXK7aS2TKNBy8exLmDOpGQoMJTRKSG9tcCMy1mKUTkoExbms8f//M1SzYWMjSrBf+4chhDs1oGHUuk3ttnAerufzKzSUAH4EN39/BTCVTfCyp1yLLNRdw3IYcPFm2kdZMUfnd2Py4b2YWUJN3TICJykGaZ2Y3u/sTuG83seqqXMBOROLR0UyF/fm8xn+RsIbNlGn+7bAhnHt5Bs8BEYmS/U2ndffpetuVGL47Uto0FpTz0US6vzV5DWnIit4/qxQ3HdqdpI82iFhE5RLcDb5nZ5fxfwZkNpADnBxVKRPZuS2EZD0/K5eWZa2icksgvzujL1Ud1pVFSYtDRRBoUVSH11J5reV51ZFd+eFJPMppqLU8Rkdrg7puAo8zsRGBAePN/3P3jAGOJyB62FpUx7tMVPPfFKsqrQlwxMovbRvWmVZOUoKOJNEgqQOuZ0ooqnv18JY99sozCskrOG9yJn2gtTxGRqHH3T4BPgs4hIv9tW3F5uPBcSWlFFecM6siPTu5F9zZNg44m0qCpAK0nqkLOG3PX8uDEXDYUlHJCnzbcNbov/TpqLU8RERFpOHaUlPPE1BU889lKSiqqOHtgdeHZs60KT5F4oAK0jnN3Pl6ymb9+sITcTUUMymzO/RcN4qgeGUFHExGR/TCzVOBToBHV4/Hr7v7bvex3AtXriiYD+e5+fOxSitQdBSUVPDltBU99tpLi8krOOLwDt5/ci17t0oOOJiK7UQFah81dvZ173l/CzLxtdMtowmOXDeWMw9uri5uISN1QBpzk7kVmlgxMM7P3d28AaGYtgMeB09x9tZm1DSirSNwq2FXBU9PyeGpaHoVllZxxeHtuO7k3fdqr8BSJRypA66AVW4oYMyGH9xduJKNpI/5w3gAuGd6Z5EQtqSIiUleElzcrCj9MDn/5HrtdBrzp7qvDr9kcu4Qi8W1naQVPT1vJk9NWsLO0ktH923Hbyb11+5FInFMBWodsKSzjkUlLeXnmalKSErh9VC9uPLY7TbSkiohInWRmiVQv4dITeMzdZ+yxS28g2cwmA+nAw+7+3F7e5ybgJoCsrKyoZhYJWlFZJc98lscTU/Mo2FXBKf3acfuoXvTv2DzoaCISAVUudUBJeSXjp+bxjynLKasMcemILH50ci/apGtJFRGRuszdq4DB4am2b5nZAHdfuNsuScAw4GQgDfjCzKbvuSa3u48DxgFkZ2fveRVVpF4oKqvk2c9X8sTUFewoqeDkvm25fVRvDs9U4SlSl6gAjWOVVSFen7OWBybmsrmwjNP6t+eu0/qofbiISD3j7jvCVzlPA3YvQNdS3XioGCg2s0+BQUDu/76LSP1UVFbJc1+s5IlPV7C9pIIT+rThx6N6M6hzi6CjichBUAEah9ydyTlb+PN7i1m6uYihWS34+xVDGdalVdDRRESklphZG6AiXHymAaOAv+6x27+Av5lZEpACjAQejG1SkWDsKCnnxRmrGT/1/wrP207uxZCslkFHE5FDoAI0zixaX8Cf31vMZ8u20rV1Y/5++VBOG6DOtiIi9VAH4NnwfaAJwGvu/q6Z3Qzg7mPdfbGZfQAsAELA+D2m6IrUO6u2FvPUtDxem72WXRVVKjxF6hkVoHFiY0EpYybk8OaXa2mRlszvzu7HZSO7kJKkzrYiIvWRuy8Ahuxl+9g9Ho8BxsQql0hQ5qzaxhOf5jHh640kJRjnDu7E9cd047AO6morUp9ErQA1s87Ac0B7qs/ajnP3h6P1eXVVcVkl/5iynHFTVxAKwU3HducHJ/akeVpy0NFEREREoqqyKsSERZsYP20FX67eQfO0ZH5wQg+uPrIrbZulBh1PRKIgmldAK4GfuvtcM0sH5pjZRHf/OoqfWWdUhZw35qxlzIc5bCks46yBHfjZaX3p3Kpx0NFEREREoqqorJLXZq3hqc/yWLt9F11aN+b35/bngmGZNE7RBD2R+ixq/8LdfQOwIfx9oZktBjoBDb4Anb5iK39492sWrd/J4M4tGHvFMIZ10X0NIiIiUr9tKNjFM5+v5KUZqyksrSS7S0t+deZhnNKvPYkJ6nch0hDE5BSTmXWl+j6XPRfYblBWby3hz+8t5oNFG+nYPJWHLxnMOYM6qsGQiIiI1GsL1xUwfuoK3l2wgZA7pw/owA3HdlNjIZEGKOoFqJk1Bd4Abnf3nXt5/ibgJoCsrKxoxwlEcVklj32yjPFT80hKNH56Sm9uPK47qcmJQUcTERERiYpQyPkkZzNPTsvj8+VbaZKSyFVHduXao7vqliORBiyqBaiZJVNdfL7o7m/ubR93HweMA8jOzvZo5om1UMh5e9467nl/CZsLy/jOkE7cdVpf2jfXTfUiIiJSP5WUV/LG3HU8PS2PFfnFdGieyt2n9+WSEVlqsigiUe2Ca8CTwGJ3fyBanxOvvlpbwG/eWciXq3cwKLM5Y68cxlBNMxEREZF6asnGnbwycw1vfbmOgl0VDMpszsOXDOaMwzuQnKhl5USkWjSvgB4NXAl8ZWbzwtt+4e7vRfEzA7etuJwxE3J4ZdZqWjdJ4d4LBnLB0EwSdGO9iIiI1DNFZZX8e/56Xpm1hvlrdpCSmMDoAe25+sguDOvSUn0uROR/RLML7jSgwfzUqQo5r8xazb0f5FBUVsl1R3fjtlG9aJaqqSYiIiJSf7g789bs4JWZa/j3gvWUlFfRq21Tfn1WP74zpBMtm6QEHVFE4pgWWqoFX60t4Fdvf8X8tQUc0b0Vvz93AL3bpQcdS0RERKTW7Cgp5+0v1/HKrDUs2VhIWnIiZw/qwMXDsxia1UJXO0UkIipAD8HO0grum5DD89NX0bpJIy2rIiIiIvWKuzN9xTZenbWa9xZupLwyxMDM5vzp/AGcM6gj6ZrpJSI1pAL0ILg77y/cyO/eWUR+URlXHdGFn47uo+m2IiIiUi9sLizljTnreHXWalZuLSE9NYlLhnfm4uGd6d+xedDxRKQOUwFaQ2u3l/Cbfy3i4yWb6d+xGeOvzmZgZougY4mIiIgcksqqEJNztvDa7DV8vGQzlSFnRNdW/OjkXpw+oANpKVq/XEQOnQrQCIVCzgszVnHP+0twh1+deRjXHNWVJLUVFxERkTrK3VmysZC3563jzbnr2FJYRkbTFK47phsXD+9MjzZNg44oIvWMCtAIrMwv5q43FjAzbxvH9srgL985nMyWjYOOJSIiIlJj7k7upiL+s2A97361gRVbiklMME7s05aLsjM5sW9brdspIlGjAnQ/QiHnqc/yuO/DHJITE7j3goFcOCxTTYZERESkzlm2uZB3F2zg3QUbWLa5iASDI7q35rqju3HagPZkNG0UdEQRaQBUgO7Duh27uOO1+XyxYisn923Ln84/nPbNU4OOJSIiIhKx5VuK+M+CDfxnwQZyNhViBiO6tuLqc/tz2oAOtElX0SkisaUCdA/uztvz1vGbtxcRcufe7w7kwmxd9RQREZG6YWV+Mf/5qvpK5+INOzGD7C4t+X/n9Of0Ae1p20wn1EUkOCpAd1NUVskv3vyKd+avJ7tLSx64aDBZrXWvp4iIiMS31VtLwkXnehat3wnA0KwW/OasfpxxeAfN4hKRuKECNOzr9Tu55aW5rNpazB2n9ub7J/QkMUFXPUVERCQ+Ld1UyPsLN/LBwo18vaG66ByS1YJfnXkYZxzegY4t0gJOKCLyvxp8AeruvDJrDb97ZxHN05J5+cYjGNm9ddCxRERERP5LSXkl01dsZUrOFqbkbmHl1hLMYFhWS3515mGcNqC9uvSLSNxr0AVoeWWIX7+9kFdnr+HYXhk8ePFgdYATERGRuPDNcilTcjczJXcLs/K2U14VIi05kSN7tOb6Y7oxur/u6RSRuqXBFqA7Ssr5/gtz+WLFVm49qSe3j+qtKbciIiISqIKSCj5bns+UnC18unQLGwpKAejdrinXHN2V43u3IbtrSxolJQacVETk4DTIAjQvv5jrn5nF2u27eOjiwZw3pFPQkURERKQBCoWcr9YVMCW3elrtl6u3E3JIT03imJ4Z3D6qDcf1bkOH5rqfU0TqhwZXgH65ejvXPjOLBDNevHEkw7u2CjqSiIiINCDrduzis6X5TFtW/bWtuBwzGNipObec2JPje7dhcOcWJCUmBB1VRKTWNagCdOG6Aq56aiatmqTw3HUj6NK6SdCRREREpB5zd1ZuLWHuqu3MXrWd6Su2kpdfDEBG00ac0KcNx/duwzE9M2itPhQi0gA0mAJ06aZCrnpqJumNknjxhpHqEiciIiK1bld5FfPX7mDu6u3MXbWduat3sK24HID0RkkM79aKK4/owjG9MujVtilm6j8hIg1LgyhAV+YXc/n4GSQlGC/deISKTxERETlk7s7a7bv+q9j8esNOqkIOQPc2TTipb1uGdWnJsC4t6dmmKQlqeCgiDVy9L0Dzi8q4fPwMKkPOqzcdQdcMTbsVERGRmiutqGLhugLmrt7OnHDBuaWwDIDGKYkMymzBzcd3Z1iXlgzp3JKWTVICTiwiEn/qfQE67tMVbCjYxb9uOYZe7dKDjiMiIiJ1QMGuCpZuKmTJxkIWritg/toCcjcVfnt1M6tVY47u0ZqhXVoyNKslfdunq2mQiEgE6nUBur24nBemr+LsQR05PLN50HFEREQkzhSWVrB0cxFLNxWSu6mI3E2FLN1UxMadpd/u0zwtmYGZzTm5bw8Oz2zO0KyWtElXwyARkYNRrwvQZz5fSUl5FT84oWfQUURERCRAxWWVLNtcXWDmhovNpZsKWV/wf4VmanICPds25agerenVLp3e7ZrSu106mS3T1CxIRKSW1NsCtKiskmc+X8kp/drRp72m3oqIiNR3u8qrWF+wi3Xbd7FyazErthSTl1/M8i1FrN2+69v9UpIS6NGmKcO7taJ3u3R6tW1Kn/bpZLZsTKKaBImIRFW9LUBfmL6Kgl0V/PBEXf0UERGp60orqthYUMr6gl1s2FHKhoJdbCgoZUNBKet37GLjzlJ2lFT812uapCTSNaMJQ7JacnF252+vama1aqz7NUVEAhK1AtTMngLOAja7+4Bofc7elFZUMX5qHsf2ymBQ5xax/GgREZGImFkq8CnQiOrx+HV3/+0+9h0OTAcudvfXY5cyusorQxSWVlBUVsn2kgq2FpWRX1RGflE5G8PF5TeF5jdrae6uZeNkOjRPo1OLNLK7tqRD8zQ6NE+lY4s0umc0oU16I02dFRGJM9G8AvoM8DfguSh+xl69NnsN+UVl3HLikFh/tIiISKTKgJPcvcjMkoFpZva+u0/ffSczSwT+CkwIIuTeVFaFKCqrpLC0+qv6+4pvt337OPx84Z7Ph7eVV4b2+RnN05Lp0DyVDs1TGdS5BR2apdKhRRodm1f/2b5ZKmkpiTH8rxYRkdoQtQLU3T81s67Rev99qagK8Y8pKxjWpSUju7WK9ceLiIhExN0dKAo/TA5/+V52vRV4Axgeo2gAPPHpCuas2k5hWcW3BeM3xeOuiqoDvj4xwUhPTaJpo+qv9NQk2qan0j0jiaapSaSHt6WnJtO0URLN05LJSG9ERtMUMpo2IjVZxaWISH0U+D2gZnYTcBNAVlbWIb/fuu27SElK4Icn9tS0GxERiWvhq5tzgJ7AY+4+Y4/nOwHnAyexnwK0tsdSgFXbqhv4NE1NokXjFDq3avxtIZmemkyT8PfNUpNo2ii5uqj85qtRMqnJCRqHRUTkfwRegLr7OGAcQHZ29t7O/NZI14wmfPST41ETOxERiXfuXgUMNrMWwFtmNsDdF+62y0PAz9y9an/FXG2PpQB/PO/w2ngbERGR/xJ4ARoNaqEuIiJ1ibvvMLPJwGnA7gVoNvBKuPjMAM4ws0p3fzvmIUVERGpBvSxARURE4p2ZtQEqwsVnGjCK6mZD33L3brvt/wzwropPERGpy6K2CJaZvQx8AfQxs7Vmdn20PktERKQO6gB8YmYLgFnARHd/18xuNrObA84mIiISFdHsgntptN5bRESkrnP3BcD/rBfm7mP3sf810c4kIiISbVG7AioiIiIiIiKyOxWgIiIiIiIiEhMqQEVERERERCQmVICKiIiIiIhITKgAFRERERERkZgwdw86w7fMbAuwqpbeLgPIr6X3agh0vGpGx6vmdMxqRserZmrzeHVx9za19F4xp7E0UDpeNadjVjM6XjWj41VzUR9P46oArU1mNtvds4POUVfoeNWMjlfN6ZjVjI5Xzeh4RYeOa83oeNWcjlnN6HjVjI5XzcXimGkKroiIiIiIiMSEClARERERERGJifpcgI4LOkAdo+NVMzpeNadjVjM6XjWj4xUdOq41o+NVczpmNaPjVTM6XjUX9WNWb+8BFRERERERkfhSn6+AioiIiIiISBxRASoiIiIiIiIxUacLUDM7zcxyzGyZmf18L8+bmT0Sfn6BmQ0NImc8ieCYXR4+VgvM7HMzGxREznhxoOO1237DzazKzC6IZb54E8nxMrMTzGyemS0ysymxzhhPIvj32NzM/m1m88PH69ogcsYLM3vKzDab2cJ9PK+f+QdJ42nNaCytGY2lNafxtGY0ntZM4OOpu9fJLyARWA50B1KA+UC/PfY5A3gfMOAIYEbQuevAMTsKaBn+/vSGfMwiOV677fcx8B5wQdC54/l4AS2Ar4Gs8OO2QeeO8+P1C+Cv4e/bANuAlKCzB3jMjgOGAgv38bx+5h/ccdV4WvvHS2NpDY7Xbvs1+LE00mOm8bTGx0vj6X8fj0DH07p8BXQEsMzdV7h7OfAKcO4e+5wLPOfVpgMtzKxDrIPGkQMeM3f/3N23hx9OBzJjnDGeRPJ3DOBW4A1gcyzDxaFIjtdlwJvuvhrA3RvyMYvkeDmQbmYGNKV6wKyMbcz44e6fUn0M9kU/8w+OxtOa0VhaMxpLa07jac1oPK2hoMfTulyAdgLW7PZ4bXhbTfdpSGp6PK6n+uxHQ3XA42VmnYDzgbExzBWvIvn71RtoaWaTzWyOmV0Vs3TxJ5Lj9TfgMGA98BVwm7uHYhOvTtLP/IOj8bRmNJbWjMbSmtN4WjMaT2tfVH/mJ9XWGwXA9rJtzzVlItmnIYn4eJjZiVQPmsdENVF8i+R4PQT8zN2rqk+qNWiRHK8kYBhwMpAGfGFm0909N9rh4lAkx2s0MA84CegBTDSzqe6+M8rZ6ir9zD84Gk9rRmNpzWgsrTmNpzWj8bT2RfVnfl0uQNcCnXd7nEn1WY2a7tOQRHQ8zGwgMB443d23xihbPIrkeGUDr4QHzAzgDDOrdPe3Y5IwvkT6bzLf3YuBYjP7FBgENMQBM5LjdS1wj1ffkLHMzPKAvsDM2ESsc/Qz/+BoPK0ZjaU1o7G05jSe1ozG09oX1Z/5dXkK7iygl5l1M7MU4BLgnT32eQe4KtzJ6QigwN03xDpoHDngMTOzLOBN4MoGehZtdwc8Xu7ezd27untX4HXgBw14wIzk3+S/gGPNLMnMGgMjgcUxzhkvIjleq6k+u42ZtQP6ACtimrJu0c/8g6PxtGY0ltaMxtKa03haMxpPa19Uf+bX2Sug7l5pZj8EJlDd/eopd19kZjeHnx9LdSe1M4BlQAnVZz8arAiP2W+A1sDj4TORle6eHVTmIEV4vCQskuPl7ovN7ANgARACxrv7XluA13cR/v36A/CMmX1F9XSYn7l7fmChA2ZmLwMnABlmthb4LZAM+pl/KDSe1ozG0prRWFpzGk9rRuNpzQU9nlr1lWgRERERERGR6KrLU3BFRERERESkDlEBKiIiIiIiIjGhAlRERERERERiQgWoiIiIiIiIxIQKUBEREREREYkJFaAi9YSZtTCzHwSdQ0RERERkX1SAitQfLQAVoCIiIrXEqun3ZZFapH9QIvXHPUAPM5tnZmOCDiMiIlIXmVlXM1tsZo8Dc4HOQWcSqU/M3YPOICK1wMy6Au+6+4Cgs4iIiNRV4fF0BXCUu08POI5IvaMroCIiIiIi/22Vik+R6FABKiIiIiLy34qDDiBSX6kAFak/CoH0oEOIiIiIiOyLClCResLdtwKfmdlCNSESERERkXikJkQiIiIiIiISE7oCKiIiIiIiIjGhAlRERERERERiQgWoiIiIiIiIxIQKUBEREREREYkJFaAiIiIiIiISEypARUREREREJCZUgIqIiIiIiEhM/H8my20HSaJRvwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 936x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# solve\n",
    "solver = pybamm.ScipySolver()\n",
    "t = np.linspace(0, 1, 100)\n",
    "solution = solver.solve(model, t)\n",
    "\n",
    "# post-process, so that the solution can be called at any time t or space r\n",
    "# (using interpolation)\n",
    "c = solution[\"Concentration\"]\n",
    "\n",
    "# plot\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n",
    "\n",
    "ax1.plot(solution.t, c(solution.t, r=1))\n",
    "ax1.set_xlabel(\"t\")\n",
    "ax1.set_ylabel(\"Surface concentration\")\n",
    "r = np.linspace(0, 1, 100)\n",
    "ax2.plot(r, c(t=0.5, r=r))\n",
    "ax2.set_xlabel(\"r\")\n",
    "ax2.set_ylabel(\"Concentration at t=0.5\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [next notebook](./3-negative-particle-problem.ipynb) we build on the example here to to solve the problem of diffusion in the negative electrode particle within the single particle model. In doing so we will also cover how to include parameters in a model. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
      "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "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.8.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}