4 ODE-integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# 4. Integrating ODEs using `odeint`\n",
    "\n",
    "First, import some stuff we'll need. We give some of them nicknames to save typing later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `odeint`, which is found in the `integrate` submodule of the `scipy` package, integrates ordinary differential equations of the form\n",
    "$$\n",
    "\\newcommand{\\yy}{\\mathbf{y}}\n",
    "\\newcommand{\\ff}{\\mathbf{f}}\n",
    "\\frac{\\mathrm{d} \\yy}{\\mathrm{d}t} = \\ff(\\yy, t),\n",
    "$$where $\\yy$ is a vector in $\\mathbb{R}^n$, or\n",
    "$$\n",
    "\\dot{\\yy} =  \\ff(\\yy, t).\n",
    "$$\n",
    "That is, it provides a numerical solution of the *initial value problem* given by this differential equation, together with an initial condition\n",
    "$$\n",
    "\\yy(t = t_0) = \\yy_0,\n",
    "$$\n",
    "returning the solution curve (or trajectory) $\\yy(t)$ as a function of time $t$, within a certain (finite!) range of values of $t$.\n",
    "\n",
    "Writing $\\yy = (y_1, \\ldots, y_n)$, we have a system of the form\n",
    "\\begin{align}\n",
    "\\dot{y_1} &= f_1(y_1, \\ldots, y_n, t) \\\\\n",
    "\\dot{y_2} &= f_2(y_1, \\ldots, y_n, t) \\\\\n",
    "&\\ldots \\\\\n",
    "\\dot{y_n} &= f_n(y_1, \\ldots, y_n, t)\n",
    "\\end{align}\n",
    "\n",
    "Note that in fact, *any* ordinary (system of) ordinary differential equation(s) can be written in this form by defining new variables to represent higher-order derivatives,  \n",
    "and incorporating them into the vector $\\yy$ of variables.\n",
    "\n",
    "To use `odeint`, we must define at least the Python function `f` which implements the mathematical function $\\ff$. It must take a vector `y` and a time `t`,\n",
    "and return the new vector `f(y,t)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest case is when the spatial dimension is $n=1$. For example, we could wish to solve the ODE\n",
    "$\n",
    "\\dot{y} = ay,\n",
    "$\n",
    "whose exact solution is, of course, $y(t) = y_0 \\exp(a t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXhx2EQAABDQRQQEGsioqoqKN8BRUrfFUs\nbYWgVH4VrRtVwbYQqq2V2m8pWrW0iLiVIrYVFIEihGoLsriAbFKRLWyGJYo0QJLP74+5iUPInkkm\nGd7Px2Me3nzuOXfOTXA+c8659x5zd0RERIpSK9YNEBGR6k2JQkREiqVEISIixVKiEBGRYilRiIhI\nsZQoRESkWCUmCjObYma7zWxVROwcM1tiZh+a2TIzuyBi3xgz22hm68ysb0S8h5mtMrNPzWxiRLye\nmU0P6iwxs+SIfSlB+Q1mNjQ6pywiImVRmh7FVKBfgdgEYJy7nweMA34NYGbdgFuArsC1wDNmZkGd\nZ4Hh7t4F6GJmecccDuxz987AxODYmFkiMBa4ELgIGGdmTct1liIiUm4lJgp3fw/YXyCcC+R9aDcD\n0oPtG4Dp7p7t7puBjUBPM2sDNHH35UG5F4GBwfYAYFqwPRO4KtjuB8x390x3PwDMB64pw7mJiEgU\n1ClnvfuBeWb2G8CAS4J4ErAkolx6EMsGtkfEtwfxvDrbANw9x8wyzax5ZLzAsUREpAqVdzL7TuBe\nd08mnDSej16TsJKLiIhIVSlvjyLF3e8FcPeZZvanIJ4OtIso1zaIFRWPrLPDzGoDCe6+z8zSgVCB\nOosKa4yZ6YFVIiLl4O4lfjkvbY/COPabfrqZXQFgZn0Iz0UAzAIGB1cydQQ6AcvcfReQaWY9g8nt\nocAbEXVSgu1BwMJgex5wtZk1DSa2rw5ihXL3uH2NGzcu5m3Q+en8TsTzi+dzcy/99+sSexRm9irh\nb/YtzGwr4auc7gAmBT2ALGBE8GG91sxmAGuBo8BI/6Y1dwEvAA2AOe4+N4hPAV4ys43AXmBwcKz9\nZvYosAJwYLyHJ7VFRKQKlZgo3P17Rey6oLCguz8OPF5IfCVwdiHxw4QvqS3sWC8QTi4iIhIjujO7\nBgiFQrFuQqXS+dVs8Xx+8XxuZWFlGaeqrszM4+E8RESqkpnhpZjMLu9VTyIidOjQgS1btsS6GVKC\n9u3bs3nz5nLXV49CRMot+EYa62ZICYr6O5W2R6E5ChERKZYShYiIFEuJQkREiqVEISJxq3v37vzz\nn/+MdTMq3eLFi2nXrl3JBctJiUJE4tYnn3zC5ZdfXmK5jh07snDhwhLLVZbbbruNsWPHVugY3yz9\nE31KFCIilSgnJyfWTagwJQoRib70dOjfP/xKTy+5fCUdI6+nMH78eL7zne+QkpJCQkICZ599Nh98\n8AEAQ4cOZevWrXz7298mISGBJ598EoClS5dy6aWXkpiYyHnnncfixYvzj7t582auuOIKmjZtSt++\nfbn77rsZMmQIAFu2bKFWrVo8//zztG/fnj59+gBwyy23cMopp5CYmEgoFGLdunUA/PGPf+SVV15h\nwoQJJCQkMGDAAAB27tzJzTffTKtWrTj99NN56qmn8t8/KyuLYcOG0bx5c7p3787y5cupVLF+emGU\nnoDoIlL1ivx/77rr3CH8uu668h08Csfo0KGDv/POO56amuoNGzb0uXPnem5uro8ZM8Z79ep1TLmF\nCxfm/5yenu4tWrTwuXPnurv7ggULvEWLFp6RkeHu7hdffLE/9NBDfvToUX/vvfc8ISHBhwwZ4u7u\nmzdvdjPzlJQUP3TokGdlZbm7+9SpU/3rr7/2I0eO+P333+/nnntu/vsNGzbMf/azn+X/nJub6+ef\nf74/9thjnp2d7Z9//rmffvrpPn/+fHd3f/jhh/3yyy/3AwcO+Pbt27179+7erl27In8PRf2dgniJ\nn7HqUYjICaF3797069cPM2PIkCGsWrXqmP0ecUPayy+/TP/+/enXrx8Affr04YILLmDOnDls27aN\nFStWMH78eOrUqcOll17KDTfccMyxzIzx48fTsGFD6tevD8CwYcNo1KgRdevWZezYsXz88cd89dVX\nhbZ1+fLlZGRk8JOf/ITatWvToUMHfvCDHzB9+nQAXnvtNX7605/StGlTkpKSuOeee6L2eyqMHuEh\nItE3eTKMGPHNdqyOEaFNmzb5240aNSIrK4vc3Fxq1Tr++/KWLVuYMWMGs2fPBsJJJDs7m6uuuood\nO3bQvHlzGjRokF++Xbt2bN++/ZhjtG3bNn87NzeXRx55hJkzZ5KRkYGZYWZkZGTQpEmTQt8/PT2d\n5s2b579/bm5u/sT8jh07jjl++/bty/MrKTUlChGJvqQkeOut2B+jlApeMdSuXTuGDh3KH/7wh+PK\nbt26lX379pGVlZWfLLZt23bcMSJ/fvXVV5k9ezYLFy4kOTmZzMxMEhMT83sxhb3/aaedxoYNGwpt\n76mnnsq2bdvo2rUrQKU/b0tDTyJyQoocamrTpg2bNm3K//nWW29l9uzZzJ8/n9zcXLKysli8eDE7\nduwgOTmZCy64gNTUVI4ePcqSJUvyex6FHRvgq6++on79+iQmJvL1118zZsyYY5JD69atj3n/nj17\n0qRJEyZMmEBWVhY5OTmsWbOGFStWADBo0CAef/xxDhw4wPbt23n66afLdvJ5FwqUUomJwsymmNlu\nM1tVIP4jM1tnZqvN7FcR8TFmtjHY1zci3sPMVpnZp2Y2MSJez8ymB3WWmFlyxL6UoPwGMxta6rMS\nEaH4ewsi940ePZpHH32U5s2b83//93+0bduWN954g1/+8pecfPLJtG/fnieffJLc3FwAXnnlFf79\n73/TsmVLxo4dy+DBg/PnIgp736FDh5KcnExSUhLdu3fnkksuOWb/8OHDWbNmDc2bN+fGG2+kVq1a\nvPnmm3z00Ud07NiRVq1acccdd/Dll18CMG7cOJKTk+nYsSPXXHMNQ4eW8eNxxAg2z1lT6uIlPj3W\nzHoDB4EX3f1bQSwEPAJc5+7ZZtbS3TPMrCvwKnAh0BZYAHR2dzez94G73X25mc0Bfufu88zsTuBs\ndx9pZt8B/tfdBwfrZK8AehBer3sl0MPdMwtpo5d0HiISfXp6bNjgwYPp2rUr48aNi3VTCnXc36l/\nf/4ypzGDmRGdp8e6+3vA/gLhO4FfuXt2UCYjiA8Aprt7trtvBjYCPc2sDdDE3fMu9n0RGBhRZ1qw\nPRO4KtjuB8x390wPr5U9H7impPaKiFS2FStWsGnTJtyduXPnMmvWLAYOHFhyxepi8mTWd/52qYuX\nd46iC3C5mS01s0Vmdn4QTwK2RZRLD2JJQOQlAduD2DF13D0HyDSz5sUcS0Qkpnbt2kUoFKJJkybc\nd999PPfcc5xzzjmxblbpJSWx/vxbS128vFc91QES3b2XmV0IvAacVs5jFVSuB5akpqbmb4dCIa11\nKyKV5vrrr+f666+PdTPKLC0tjbS0NAAWLSp9vfImim3AXwGCOYccM2tB+Ft/ckS5tkEsHWhXSJyI\nfTvMrDaQ4O77zCwdCBWoU+SpRSYKERE5Xt6X6Nxc+PWvAcaXql5ph56MY7/p/51gLsHMugD13H0v\nMAv4TnAlU0egE7DM3XcRHlLqaeHLAYYCbwTHmgWkBNuDgLxHOM4DrjazpsHE9tVBTEREKmD7dmjW\nrPTlS+xRmNmrhL/ZtzCzrcA44HlgqpmtBg4T/uDH3dea2QxgLXAUGBlxOdJdwAtAA2COu88N4lOA\nl8xsI7AXGBwca7+ZPUr4yicHxgeT2iIiUgHr18OZZ8KOHaUrX+LlsTWBLo8ViQ1dHlszFPw7TZoE\nGzbAM89YdC6PFRGR+JLXoygtJQoRkTKo6LKjTZo0YfPmzRVqQ0VXxFOiEBGpZKVddvTKK6/k+eef\nPyb21Vdf0aFDh0poVemtXw9nnFH68koUIiInkMzM8CviKeUlUqIQkbj1xBNP0KlTJxISEujevTt/\n//vfAZg2bRqXXXYZDz74IM2bN+f0009n7ty5+fVeeOEFunXrRkJCAp06dWJyEethPPnkk9x8883H\nxO69917uv/9+fvrTn/Luu+9y9913k5CQkL+4UK1atfKfFJuVlcWoUaPo0KEDiYmJXH755Rw+fBg4\nfunUtWvXRuV3smFDuDdRyDIcRVKiEJG41alTJ/71r3/x5ZdfMm7cOIYMGcLu3bsBeP/99+natSt7\n9+7lwQcfZPjw4fn1WrduzZw5c/jyyy+ZOnUq999/Px999NFxx7/11luZN29e/lNdc3JymD59Oikp\nKTz22GNcdtllPP3003z55ZdMmjQJOHbYatSoUXz44YcsXbqUffv2MWHChPyFlK677jo+++wz9uzZ\nQ48ePfj+978fld/Jhg1lm58AJQoRqURm0XmV10033UTr1q2B8BoOnTp1YtmyZQB06NCB22+/HTMj\nJSWFXbt2sWfPHgCuvfba/HmEyy67jL59+/Luu+8ed/w2bdpw+eWX89prrwHw9ttvc/LJJ3PuuecW\n2aa8y1TdnalTpzJp0iTatGmDmdGrVy/q1q0LlG3p1LIo60Q2KFGISCVyj86rvF588UXOO+88EhMT\nSUxMZM2aNWRkhB92Hbk0asOGDXF3Dh48CIQ/8C+++GJatGhBYmIib7/9dn69goYOHcrLL78MhNep\nGDJkSKnalpGRweHDhznttOMfk5ebm8vo0aPp1KkTzZo1o2PHjvlLp1aUEoWISGDr1q2MGDGCZ555\nhv3797N//37OOuusEm8QPHLkCDfffDMPPfQQX3zxBfv37+faa68tst7AgQNZtWoVa9as4c033zxm\niKi4q6NatmxJgwYN+Oyzz47bF7l06oEDB9i8eTPuHpWbG5UoREQCX3/9NbVq1aJly5bk5uYydepU\nPvnkkxLrHTlyhCNHjtCyZUtq1arF22+/zfz584ssX79+fW666Sa+973vcdFFF9E24nKigkucRjIz\nbrvtNh544AF27txJbm4uS5cu5ciRIyUunVpe2dmwaRN07ly2ekoUIhKXunbtyqhRo+jVqxdt2rRh\nzZo19O7du8jyeR/EjRs3ZtKkSQwaNIjmzZszffp0BgwYUOx7paSksHr16uOWJL333nt57bXXaNGi\nBffdd98x7wPhq6bOPvtsLrzwQlq0aMHo0aNx9xKXTi2vzz+HU06Bhg3LVk/PehKRctOznsK2bdtG\n165d2bVrF40bN451c46T93eaPRuefRbmzDkmrmc9iYhUptzcXH7zm98wePDgapkkIpX1juw85V24\nSETkhHfo0CFat25Nx44defvtt2PdnBKtXw89e5a9nhKFiEg5NWrUKCr3NlSVDRugwDRKqWjoSUTk\nBFGeS2OhFInCzKaY2W4zW1XIvlFmlmtmzSNiY8xso5mtM7O+EfEeZrbKzD41s4kR8XpmNj2os8TM\nkiP2pQTlN5hZOfKgiIgAZGRATg60alX2uqXpUUwF+hUMmllbwutYb4mIdQVuAboC1wLP2DfXgj0L\nDHf3LkAXM8s75nBgn7t3BiYCE4JjJQJjgQuBi4BxZta0zGcoIiL5vYny3I5RYqJw9/eA/YXs+i3w\nYIHYAGC6u2e7+2ZgI9DTzNoATdx9eVDuRWBgRJ1pwfZM4Kpgux8w390zg7Wy5wPXlOqsRKRKtG/f\nHjPTq5q/2rdvX+5hJyjnZLaZ3QBsc/fV33QYAEgClkT8nB7EsoHtEfHtQTyvzjYAd88xs0wLD2Xl\nxwscS0SqiYqu1CZV58c/rsJEYWYNgUcIDztVhnLdp56ampq/HQqFCIVCUWqOiEjNt349NGuWRmpq\nWpnrlqdHcTrQAfg4mH9oC3xgZj0Jf+tPjijbNoilA+0KiROxb4eZ1QYS3H2fmaUDoQJ1FhXVqMhE\nISIix1q/Hn7zmxBnnBHKj40fP75UdUt7eawFL9z9E3dv4+6nuXtHwsNI57n7HmAW8B0LX8nUEegE\nLHP3XUCmmfUMkstQ4I3g2LOAlGB7ELAw2J4HXG1mTYOJ7auDmIiIlEFWFmzfDoU80bxUSuxRmNmr\nhL/ZtzCzrcA4d58aUcT5JomsNbMZwFrgKDAy4iFMdwEvAA2AOe6et+7gFOAlM9sI7AUGB8fab2aP\nAiuC9xgfTGqLiEgZ/Oc/0KEDBGsilZkeCigiEudefx1eegmCJcPzmR4KKCIiUP47svMoUYiIxDkl\nChERKVZFE4XmKERE4pg7JCTA1q2QmHjsPs1RiIgIO3bASScdnyTKQolCRCSOVXTYCZQoRETiWnmX\nP42kRCEiEsfUoxARkWJt2KBEISIihUlPh/79Wf/uHs5stqtCh1KiEBGJRyNGcHDOYjKyGpP86B0V\nOpQShYhInPqYcziLNdS23AodRzfciYjEo/R0JvWdzdqvknluyTmQdPwCoaW94a5cS6GKiEg1l5TE\nygt+SO/eVHgRaQ09iYjEqZUr4fzzK34cDT2JiMShr7+Gk0+GAwegXr3Cy0TtWU9mNsXMdpvZqojY\nBDNbZ2YfmdnrZpYQsW+MmW0M9veNiPcws1Vm9qmZTYyI1zOz6UGdJWaWHLEvJSi/wcyGltRWEREJ\n+/hj6Nat6CRRFqUZepoK9CsQmw+c5e7nAhuBMQBm1g24BegKXAs8E6yRDfAsMNzduwBdzCzvmMOB\nfe7eGZgITAiOlQiMBS4ELgLGmVnTcp2liMgJJlrDTlCKROHu7wH7C8QWuHve9VZLgbbB9g3AdHfP\ndvfNhJNITzNrAzRx9+VBuReBgcH2AGBasD0TuCrY7gfMd/fMYK3s+cA1ZTw/EZETUpUmilK4HZgT\nbCcB2yL2pQexJGB7RHw738zD59dx9xwg08yaF3MsEREpQTQTRYUujzWznwBH3f3P0WlO+LDlqZSa\nmpq/HQqFCIVCUWqOiEjNcugQfPYZdO9+bDwtLY20tLQyH6/cicLMhgHX8c1QEYS/9beL+LltECsq\nHllnh5nVBhLcfZ+ZpQOhAnUWFdWeyEQhInIi+/hj6NoV6tc/Nl7wS/T48eNLdbzSDj0ZEd/0zewa\n4EHgBnc/HFFuFjA4uJKpI9AJWObuuwgPKfUMJreHAm9E1EkJtgcBC4PtecDVZtY0mNi+OoiJiEgx\nojnsBKXoUZjZq4S/2bcws63AOOARoB7wj+CipqXuPtLd15rZDGAtcBQYGXGDw13AC0ADYI67zw3i\nU4CXzGwjsBcYDODu+83sUWAF4MD4YFJbRESKsXIl9OoVvePphjsRkTjzrW/B88/DBRcUX660N9wp\nUYiIxJH//hdatID9+4+foygoandmi4hIzfHxx+EV7UpKEmWhRCEiEkeiPZENShQiInFFiUJERIq1\nciX06BHdY2oyW0QkTmRlQfPmsG8fNGhQcnlNZouInGBWrYIuXUqXJMpCiUJEJE5UxvwEKFGIiMQN\nJQoRESlWZSUKTWaLiMSBvInsvXuhYcPS1dFktojICWT1aujcufRJoiyUKERE4kBlDTuBEoWISFxQ\nohARkWLFNFGY2RQz221mqyJiiWY238w2mNk8M2sasW+MmW00s3Vm1jci3sPMVpnZp2Y2MSJez8ym\nB3WWmFlyxL6UoPwGMxsanVMWEYkvhw/D+vVwzjmVc/zS9CimAv0KxEYDC9z9DMJLl44BMLNuwC1A\nV+Ba4Jlg6VOAZ4Hh7t4F6GJmecccDuxz987ARGBCcKxEYCxwIXARMC4yIYmISNjq1dCpU+VMZEMp\nEoW7vwfsLxAeAEwLtqcBA4PtG4Dp7p7t7puBjUBPM2sDNHH35UG5FyPqRB5rJnBVsN0PmO/umcES\nqPOBa8pwbiIiJ4TKHHaCUqyZXYRW7r4bwN13mVmrIJ4ELIkolx7EsoHtEfHtQTyvzrbgWDlmlmlm\nzSPjBY4lIiJ50tNZ+YvVnN/kc0i/AZKi/zEZrcnsaN7tVuLNHyIiEhgxgpXbTub8tS/CiBGV8hbl\n7VHsNrPW7r47GFbaE8TTgXYR5doGsaLikXV2mFltIMHd95lZOhAqUGdRUQ1KTU3N3w6FQoRCoaKK\niojEjYPZDdjAGZzLR3wzcl+4tLQ00tLSyvwepXqEh5l1AGa7+9nBz08QnoB+wsweBhLdfXQwmf0K\n4cnnJOAfQGd3dzNbCtwDLAfeAia5+1wzGwl0d/eRZjYYGOjug4PJ7BVAD8I9nxXA+cF8RcH26REe\nInJCmvtyBo/fu5PFvUbD5MllGnoq7SM8SuxRmNmrhL/ZtzCzrcA44FfAa2Z2O7CF8JVOuPtaM5sB\nrAWOAiMjPsHvAl4AGgBz3H1uEJ8CvGRmG4G9wODgWPvN7FHCCcKB8YUlCRGRE9mi1S258kctIfWt\nSnsPPRRQRKQG69kTfv1ruOKKstctbY9CiUJEpIbKzIS2bSEjA+rXL3t9PT1WRCTOvftuuEdRniRR\nFkoUIiI11KJFcOWVlf8+ShQiIjVUVSUKzVGIiNRA+/ZBhw7h+Yl69cp3DM1RiIjEscWL4ZJLyp8k\nykKJQkSkBqqqYSdQohARqZGqMlFojkJEpIb54gvo3Dk8P1GnvE/sQ3MUIiJxKy0NeveuWJIoCyUK\nEZEapiqHnUCJQkSkxlGiEBGRIu3cCbt3wznnVN17KlGIiNQgixaFnxRbu3bVvacShYhIDVLVw06g\nRCEiUqPUuERhZveb2SdmtsrMXjGzemaWaGbzzWyDmc0zs6YR5ceY2UYzW2dmfSPiPYJjfGpmEyPi\n9cxselBniZklV6S9IiI12bZt4TUozjqrat+33InCzE4FfgT0cPdvEV5W9bvAaGCBu58BLATGBOW7\nEV4ytStwLfCMmeXd6PEsMNzduwBdzKxfEB9OeG3uzsBEYEJ52ysiUtMtWgShENSq4rGgir5dbeAk\nM6sDNATSgQHAtGD/NGBgsH0DMN3ds919M7AR6GlmbYAm7r48KPdiRJ3IY80E+lSwvSIiNVYshp2g\nAonC3XcAvwG2Ek4Qme6+AGjt7ruDMruAVkGVJGBbxCHSg1gSsD0ivj2IHVPH3XOAA2bWvLxtFhGp\nyWpcojCzZoS/8bcHTiXcs/g+UPChS9F8CFOJzyQREYlHn38Ohw/DmWdW/XtX5Ekh/wNscvd9AGb2\nN+ASYLeZtXb33cGw0p6gfDrQLqJ+2yBWVDyyzg4zqw0k5L1fQampqfnboVCIUChUgVMTEalG0tNZ\nNOA1QnYetqMTJCWVXKcQaWlppKWllbleuZ8ea2Y9gSnAhcBhYCqwHEgmPAH9hJk9DCS6++hgMvsV\n4CLCQ0r/ADq7u5vZUuCeoP5bwCR3n2tmI4Hu7j7SzAYDA919cCFt0dNjRSR+9e/PrXO+yxUs5o7r\ndsBbb0XlsKV9emy5exTuvszMZgIfAkeD/04GmgAzzOx2YAvhK51w97VmNgNYG5QfGfHpfhfwAtAA\nmOPuc4P4FOAlM9sI7AWOSxIiIvEuO7cW8+jHL/gJ0L3K31/rUYiIVHPv/CWD0f9vH8svvR8mTy73\n0FNBpe1RKFGIiFRzd94JHTrAww9H97hKFCIicSAnJ9yBeO896NQpusfWCnciInHg3/+G1q2jnyTK\nQolCRKQae/11uOmm2LZBQ08iItVUbi60bw9z51bOgwA19CQiUsMtXw6NG0O3brFthxKFiEg1lTfs\nZDF+eJEShYhINeRePeYnQIlCRKRa+vjjcLI499xYt0SJQkSkWpo5s3oMO4EShYhItVRdhp1AiUJE\npNpZuxYOHoSePWPdkjAlChGRaub11+HGG6t+beyiVJNmiIhInrz5iepCiUJEpBr5z39g92649NJY\nt+QbShQiItXI66/DwIFQu3asW/KNCiUKM2tqZq+Z2TozW2NmF5lZopnNN7MNZjbPzJpGlB9jZhuD\n8n0j4j3MbJWZfWpmEyPi9cxselBniZklV6S9IiLV3euvw803x7oVx6poj+J3hJcu7QqcA6wHRgML\n3P0MYCEwBiBYM/sWoCtwLfCMWf4Vws8Cw929C9DFzPoF8eGE19/uDEwEJlSwvSIi1VN6OluuGsam\njzK5olN6rFtzjHInCjNLAC5z96kA7p7t7pnAAGBaUGwaMDDYvgGYHpTbDGwEeppZG6CJuy8Pyr0Y\nUSfyWDOBPuVtr4hItTZiBH9dlMgNR1+n7l0jYt2aY1SkR9ERyDCzqWb2gZlNNrNGQGt33w3g7ruA\nVkH5JGBbRP30IJYEbI+Ibw9ix9Rx9xzggJk1r0CbRUSqJXd4ntv5Hq/GuinHqVPBuj2Au9x9hZn9\nlvCwU8GFIaK5UESRN7Onpqbmb4dCIUKhUBTfVkSkcr33gxc4+s8s+lxeHyZPrpT3SEtLIy0trcz1\nyr1wkZm1Bpa4+2nBz70JJ4rTgZC77w6GlRa5e1czGw24uz8RlJ8LjAO25JUJ4oOBK9z9zrwy7v6+\nmdUGdrp7q0LaooWLRKRG++53oVcvuPfeqnvPSl+4KBhe2mZmXYJQH2ANMAsYFsRSgDeC7VnA4OBK\npo5AJ2BZMDyVaWY9g8ntoQXqpATbgwhPjouIxJVdu8Kr2KWklFw2Fioy9ARwD/CKmdUFNgG3AbWB\nGWZ2O+Hewi0A7r7WzGYAa4GjwMiIbsBdwAtAA8JXUc0N4lOAl8xsI7AXGFzB9oqIVDt/+hMMGgTN\nmsW6JYXTmtkiIjGUnQ0dO8Ls2VW/9oTWzBYRqQHefBPataseCxQVRYlCRCSGfv97uOuuWLeieBp6\nEhGJkQ0b4PLLYetWqF+/6t9fQ08iItXcc8/B7bfHJkmUhXoUIiIx8PXXkJwMH3wA7dvHpg3qUYiI\nVGN//nN4zYlYJYmyUKIQEali7uFJ7JEjY92S0lGiEBGpYu+/D199BX37lly2OlCiEBGpKunp0L8/\nv79pIXd+L5NaNeQTuIY0U0QkDowYwRdzljF7Rw9uW1K91pwojhKFiEgVeoaR3MTrNK93MNZNKTVd\nHisiUkUyVu/kjPNPYtkl93P6Kz+HpKSSK1Wi0l4eq0QhIlJFRo2CrKzwFU/VgRKFiEg1snUrnHce\nfPIJnHJKrFsTpkQhIlKNDB8ObdrAL34R65Z8o7SJoqILF4mISAnWrYNZs2Djxli3pHwqfNWTmdUy\nsw/MbFaYbJdTAAAPy0lEQVTwc6KZzTezDWY2z8yaRpQdY2YbzWydmfWNiPcws1Vm9qmZTYyI1zOz\n6UGdJWaWXNH2iohUtZ/9DB58sPquYFeSaFweey/h5U3zjAYWuPsZhNe4HgNgZt0IL4vaFbgWeCZY\nIxvgWWC4u3cBuphZvyA+HNjn7p2BicCEKLRXRKTKLF8OS5fC3XfHuiXlV6FEYWZtgeuAP0WEBwDT\ngu1pwMBg+wZgurtnu/tmYCPQ08zaAE3cfXlQ7sWIOpHHmgn0qUh7RUSq2pgxMHYsNGoU65aUX0V7\nFL8FHgQiZ5Jbu/tuAHffBbQK4knAtohy6UEsCdgeEd8exI6p4+45wAEza17BNouIVIkFC8JXO912\nW6xbUjHlnsw2s/7Abnf/yMxCxRSN5uVIRc7Op6am5m+HQiFCoVAU31ZEpGzcw72JRx+FunVj3Zqw\ntLQ00tLSylyv3JfHmtkvgVuBbKAh0AT4G3ABEHL33cGw0iJ372pmowF39yeC+nOBccCWvDJBfDBw\nhbvfmVfG3d83s9rATndvVaApujxWRKqP9HQYMYLXd17CL44+yIqP61Xbh/9V+sJF7v6Iuye7+2nA\nYGChuw8BZgPDgmIpwBvB9ixgcHAlU0egE7AsGJ7KNLOeweT20AJ1UoLtQYQnx0VEqq8RI8ieM4+f\nfHgTjzd8tNomibKojPsofgXMMLPbCfcWbgFw97VmNoPwFVJHgZER3YC7gBeABsAcd58bxKcAL5nZ\nRmAv4YQkIlKtTeIekkinb8sPYt2UqNCd2SIiUbR+8W56X92A9y95oFo8+K84eoSHiEgVy84Or4Od\nklIzljmt9DkKERE51oQJkJAAP/xhrFsSXepRiIhEwapV0KcPrFwJyTXkYUPqUYiIVJEjR2Do0HCP\noqYkibJQohARqaDHHoN27WDYsFi3pHLoMeMiIhWwYgX84Q/w0UdgJQ7i1ExKFCIiZRXcfZ2VU5eh\nm15j4sS61WbVusqgyWwRkbLq3x/mzOEhnmBTm0t4bUfvGtmb0GS2iEglepP+vMytPHvW72tkkigL\n9ShERMpo5dt7uGZgfd68YDwXzRhVre++Lo7uzBYRqQRbt8LFF8NTT8GNN8a6NRWjoScRkSjLzAxP\nT4waVfOTRFmoRyEiUgpHj8J110GXLvD00/FxKayGnkREosQd7rgDdu2Cv/8d6sTJjQUaehIRqaj0\ndOjfn8fPnMYH7x9h+vT4SRJloUQhIlKUESP485wE/vBpiDdbDadx41g3KDbKnSjMrK2ZLTSzNWa2\n2szuCeKJZjbfzDaY2TwzaxpRZ4yZbTSzdWbWNyLew8xWmdmnZjYxIl7PzKYHdZaYWRw+bktEqqtX\n06/gPibyJtdzaoN9sW5OzFSkR5ENPODuZwEXA3eZ2ZnAaGCBu59BeI3rMQBm1o3wsqhdgWuBZ4I1\nsgGeBYa7exegi5n1C+LDgX3u3hmYCEyoQHtFRErtt7+Fh/c8wDuXjefs65Jh8uRYNylmyj3a5u67\ngF3B9kEzWwe0BQYAVwTFpgFphJPHDcB0d88GNgfrYPc0sy1AE3dfHtR5ERgIzAuONS6IzwSeLm97\nRURKwx1Gj4ZZs+BfS+uQnPxMrJsUc1GZljGzDsC5wFKgtbvvhnAyMbNWQbEkYElEtfQglg1sj4hv\nD+J5dbYFx8oxswNm1tzdT9w+oIhUmqNH4Qc/gE8/hffegxYtYt2i6qHCicLMGhP+tn9v0LMoeJ1q\nNK9bLfIyrtTU1PztUChEKBSK4tuKSLz7+msYNAhq1YJ33oFGjWLdouhLS0sjLS2tzPUqdB+FmdUB\n3gTedvffBbF1QMjdd5tZG2CRu3c1s9GAu/sTQbm5hIeVtuSVCeKDgSvc/c68Mu7+vpnVBna6e6tC\n2qH7KESkfNLTyRj2Y/ovT6Vbv7ZMfvkk6taNdaOqRlXdR/E8sDYvSQRmAcOC7RTgjYj44OBKpo5A\nJ2BZMNeRaWY9g8ntoQXqpATbgwhPjouIRM3SQb/hogWPcVXmX3n+q1tOmCRRFuUeejKzS4HvA6vN\n7EPCQ0yPAE8AM8zsdsK9hVsA3H2tmc0A1gJHgZER3YC7gBeABsAcd58bxKcALwUT33uBweVtr4hI\npOxs+OUv4fcrf8Yz/ICb+CvYdbFuVrWkR3iIyAln0yYYMgQaNoRpv9pJ0rgfhHdMnlxjHxleHnrW\nk4hInmDpUnd4qd9LjHqsOWPGwH33hSevT1RKFCIiefr3Z/+cf/NDnmNN44t45b0OnHNOrBsVe3oo\noIgI4Xsjnt58PWeynjbsYvml9ytJlJF6FCISl9zhjTfg4YehfZssfp37Y85J+PyEm4cojoaeROTE\nEsxDACy7cyo/ntCK/fvhySehX78S6p6glChE5MTSvz+fztlIKqksrt+Xn/++JcOGQe3asW5Y9aU5\nChE5IbjD/PnQf3kqvXmPM1nPhiv+H8OHK0lEywm4VpOI1FgRw0uHfvdHXnrnVCZNCieEex/qxMwF\nt9Ow9pET+pHglUFDTyJSc/Tvz/o5nzGV23i+3g+55Jqm3HcfhEJgJQ6gSEEaehKRmilYp5r+/cPb\nwGefhR+3cc67T9GHd8imDksvfoA33oArr1SSqGzqUYhI9dK/P8yZwxaSmXHmWP5y0nC2bYObb4bv\nXPUFvafcRi1zXeYaBbrqSUSqt4j5BiZP5ssmSSxeDAvueYMFmzuxh1b8b7uVfGfqNVxxBdTRjGrU\nKVGISLX2db8bWTF/Lwu5igWJg1h1tBu9esH/XJjJ//xzLOcmbKL2H59Tr6ESKVGISOwFvYYjuXVY\nfc8fWba5FcuXw/LlsGldFt1zPuZKFnF1zy+5JO2XNGwY6wafWJQoRKRqRAwh5Tw7mc+PJLF2Laxb\nB2ufWsDa9KaspRunNf6CCwd1oGdPuPBCOLt5OvXu/mboST2HqhdXicLMrgEmEr5Ka0recqoR+5Uo\nRCqqwJzBMR/cEfuyJk1mS3YSn39O+PXkTDZvymUDZ7Cx1hm0ateAbt2ga1fo9o/f0XX1X/gWq2h8\n3RXw1lsxODEpStwkCjOrBXwK9AF2AMuBwe6+PqJMXCeKtLQ0QqFQrJtRaXR+laCoD/1iksHhawaw\nd95y9tCKnRfcwI4f/pydO2HHDtj513+zY3dtttGOfbVa0u60enToAB07Qu7bo+mzfQud2ciZfdvT\neN7rJbejhoj3f5ulTRQ14TqCnsBGd98CYGbTgQHA+mJrxZF4/8da486vjB/CaWlphDp3LtW39Yrs\ny26dxMGDhF/f/zlf/XMPmTTlwLenkXnXIxw4AJnP/ZMD/+nHAZqRcd4OMjokkZEBGRmQdXAmLcjg\nZL7g1M1HOeU9OPVU6NYN/uf9RZy6ezZJpHNqv3OpPWd2fjNSHzjCdzd8CbSGyZOO/V0lJdXoXkSN\n+7dZSWpCokgCtkX8vJ1w8jhWenql/M9X4X3RON4ZZxx/rtWtjRXZV5pzi3Ebc09JIjcXcnIgZ/iP\nyJn3T3KoTc6wUWS/PJ2cHMge8hNyFq0lmzoc/e6jHH3qObKzYft2ePfm35G99BBHqcuRAX/kyCOp\nHDkSXivhyK9ncmR1Moepz+E+88n67m0cPgxZWXD49Q/J2nojWTTg0PlbOHROEv/9Lxw6BIc25PDf\ng09zkMYcTE7gCNC4cfjVJPOnnMQemnGAZjsa0PRf0KwZNLMcTuczmnGAlu030fKpC2nZElq2hISv\n9mD/L/L3EfF3+d9hMOLfQAuY/Nyxf7OEhBqdDKRkNSFRlMr156VDz4h/2cvS4Ys7w9tVtM+/GBns\n2wEXBvuW74DC4qXc5xgb0/7Gssj+04od+Bd3B/V2wvkR9VbuxL/4UXj73F0F9u3CM+75Zl+PYN8H\nuyDjHhyDc3fDeUnkj+R9uAf23hfed84XcG7Evo+/gL0PBPsy8LPDx3MHVu/F9z0ULvetffhZEfvW\n7Mf3jwn/fHYmW0+CefPCRX39V/iBseFjnnUQ7xzU2XgI//LnOIZ3y8I7huPu4J8fxQ8+gWPknpGL\ntw3Hc3PBdxi5h54K7+tUl9yWQdwhd28jco9MJZda5HasT24j8pNBblZLcnP/Gk4GbcP/m9SqFX6m\nUO2cP1Obo9Qhm9qLa1HnW+Fr/Gvvm0AdDlKbHOp+WJ+6Q8Px3bth/aEU6nItdcim/tam1H0J6tUL\nXvu7Uo961OMIDY42oX4OnHQStGgB9Rtvoz4f0IhDNEzqRqP7L6FRo/Baz41GPU7Dd+fRmIM0vro3\nDd/+6zd3KKfXghFjw9uRH/qjroQRfw7ivzw2GTQt5tt/De8ZSMXUhDmKXkCqu18T/Dwa8MgJbTOr\n3ichIlJNxctkdm1gA+HJ7J3AMuC77r4upg0TETlBVPuhJ3fPMbO7gfl8c3mskoSISBWp9j0KERGJ\nrRr/mHEzu8bM1pvZp2b2cKzbE01mNsXMdpvZqli3pTKYWVszW2hma8xstZndE+s2RYuZ1Tez983s\nw+DcxsW6TZXBzGqZ2QdmNivWbYk2M9tsZh8Hf8NlsW5PtJlZUzN7zczWBf8PXlRk2ZrcoyjNzXg1\nmZn1Bg4CL7r7t2LdnmgzszZAG3f/yMwaAyuBAXH092vk7oeCebZ/Afe4e1x94JjZ/cD5QIK73xDr\n9kSTmW0Cznf3/bFuS2UwsxeAxe4+1czqAI3c/cvCytb0HkX+zXjufhTIuxkvLrj7e0Bc/iMFcPdd\n7v5RsH0QWMexF2zWaO5+KNisT3g+sOZ+KyuEmbUFrgP+FOu2VBKj5n9GFsrMEoDL3H0qgLtnF5Uk\noOb/Egq7GS9uPmhOJGbWATgXeD+2LYmeYFjmQ2AX8A93Xx7rNkXZb4EHibMEGMGBf5jZcjO7I9aN\nibKOQIaZTQ2GDiebWZHP7q3piULiQDDsNBO4N+hZxAV3z3X384C2wEVm1i3WbYoWM+sP7A56hBa8\n4s2l7t6DcK/prmAoOF7UAXoAvw/O8RAwuqjCNT1RpAPJET+3DWJSQwRjozOBl9z9jVi3pzIEXfpF\nwDWxbksUXQrcEIzj/xm40sxejHGbosrddwb//QL4G4U9Oqjm2g5sc/cVwc8zCSeOQtX0RLEc6GRm\n7c2sHjAYiLerL+L121qe54G17v67WDckmsyspZk1DbYbAlcTRw+ydPdH3D3Z3U8j/P/dQncfGut2\nRYuZNQp6upjZSUBf4JPYtip63H03sM3MugShPsDaospX+xvuihPvN+OZ2atACGhhZluBcXmTT/HA\nzC4Fvg+sDsbyHXjE3efGtmVRcQowLbgyrxbwF3efE+M2Sem1Bv4WPB6oDvCKu8+PcZui7R7gFTOr\nC2wCbiuqYI2+PFZERCpfTR96EhGRSqZEISIixVKiEBGRYilRiIhIsZQoRESkWEoUIiJSLCUKEREp\nlhKFiIgU6/8DeJSgEPP3/0gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105d8afd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 2.0\n",
    "def f(y, t):\n",
    "    return a*y\n",
    "times = np.linspace(0,6)\n",
    "y0 = 1.0\n",
    "y_output = odeint(f, y0, times)\n",
    "y_exact = y0*np.exp(a*times)\n",
    "plt.plot(times, y_output, 'r.', label='integrated')\n",
    "plt.plot(times, y_exact, label='analytical')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example E5-6, Calculating X in a reactor with pressure drop\n",
    "See page 183-188 in the textbook \"Essentials of Chemical Reaction Engineering\" by Fogler.\n",
    "\n",
    "This time we have two variables to integrate (conversion $X$ and pressure ratio $y$), so our function must accept and return a vector. We integrate over weight of catalyst $W$, not time $t$.\n",
    "The differential equations are:\n",
    "\n",
    "$$\\frac{\\mathrm{d} X}{\\mathrm{d} W} = \\frac{-r_A'}{F_{A0}}$$\n",
    "$$\\frac{\\mathrm{d} y}{\\mathrm{d} W} = -\\alpha \\frac{(1+\\epsilon X)}{2y}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dYvecdV(Yvec,W):\n",
    "    X, y = Yvec # unpack values from Y vector into three variables\n",
    "    # Explicit equations\n",
    "    epsilon = -0.15\n",
    "    kprime = 0.0074\n",
    "    Fa0 = 0.1362\n",
    "    alpha = 0.0367\n",
    "    raprime = -1*kprime*(1-X)/(1+epsilon*X) * y\n",
    "    f = (1+epsilon*X)/y # ratio of v / v0\n",
    "    # Differential equations\n",
    "    dXdW = -raprime / Fa0\n",
    "    dydW = -1*alpha * (1 + epsilon * X)/ (2*y) \n",
    "    return [dXdW, dydW]  # pack up three differential equations into the dYvecdV vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the steps you want to integrate over. In this case it's Weight not Time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "W_output = np.linspace(0, 27., num=50) # Range for the weight of catalyst of the reactor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the initial conditions, as a vector of the two values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Y0 = [0., 1.] # Initial values for the dependent variables i.e X and y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the integration, and extract the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Y_result = odeint(dYvecdV, Y0, W_output) # perform the integration\n",
    "X, y = Y_result.T # extract columns by transposing matrix and splitting into two variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEPCAYAAAC5sYRSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXh1XWkARIQgghhC0sAUR2gVCBCpVq3SpK\n0VZb2y9u9Su2X61Crd9a/Cq2Vv0pdaFYN6QtooiISKCALCKEPRACISSQkD0hJEByfn+cyWQhIZN1\nMpPP8/E4j3vnzpk75zLDe27OvfdcMcaglFLKu7RwdwOUUkrVPw13pZTyQhruSinlhTTclVLKC2m4\nK6WUF9JwV0opL1RtuItIWxHZLiK7RWSfiCyopM5kEckSke8c5XcN01yllFKuaFVdBWNMoYhMMcbk\ni0hLYIuIrDHG7KhQdZMx5ocN00yllFI14VK3jDEm3zHbFvuDUNmVT1JfjVJKKVU3LoW7iLQQkd3A\nGWCdMWZnJdXGicgeEVktIoPqtZVKKaVqxNU992JjzAigJzCmkvDeBfQyxgwHXgFW1m8zlVJK1YTU\ndGwZEXkKOGeMWXyFOseBkcaYjArLdSAbpZSqBWNMjbq+XTlbpquI+Djm2wHTgMMV6gSUmR+N/dEo\nF+xlGui1ZcGCBW5vg26fbl9z27bmsH21Ue3ZMkAQ8HcRaYH9MfjIGPO5iNxvs9osAW4VkV8BF4Hz\nwI9r1RqllFL1wpVTIfcBV1ey/I0y868Cr9Zv05RSStWWXqFaj6KiotzdhAal2+e5vHnbwPu3rzZq\nfEC1Tm8mYhrz/ZRSyhuICKaGB1Rd6XNvcL179yYhIcHdzVBuFhoayokTJ9zdDKW8QpPYc3f8KjVa\nO1TTpN8DpSpXmz137XNXSikvpOGulFJeSMNdKaWqkpvr7hbUmoZ7A3njjTd49NFH3d2MejFkyBA2\nbdpU5/U89thjvP766/XQIqUawblz0LMn5OdXX7cJ0gOq1Th37hxDhgzhj3/8I7NnzwYgLy+PwYMH\n89JLL3HzzTdf9pqLFy8SHh7Ojh07CAwMbOwmN1lnzpxh9OjRxMfH06rV5SdqNeXvgWqGVq+GF16A\nDRvc3RI9oNoQOnTowBtvvMHDDz9Meno6APPnz2f06NGVBjvAJ598QkREhMcEe1FRUaO8T2BgIBER\nEaxatapR3k+pOvniC/j+993dilrTcHfB9OnTueGGG3jwwQfZuHEjK1as4LXXXquy/po1a5g8eXK5\nZZs3b2bChAn4+voSGhrKsmXLAMjJyWHu3Ll0796dsLAw/vd//9f5mr///e9MnDiR+fPn4+fnR3h4\nOGvXrgVg+fLljBo1qtx7vPTSS9x0000AXLhwgccee4zQ0FCCgoL4r//6LwoLCwHYuHEjISEhPP/8\n8wQFBfGzn/2M9PR0Zs2aha+vL/7+/uXaHxYWxtdff+1c7yOPPEJwcDA9e/bk17/+NRcvXiy33sWL\nFxMQEEBwcDBLly4t18bJkyezevVql//tlXKbL76A6693dytqTcPdRYsXLyY6Oppbb72VF198kW7d\nulVZd9++fQwYMMD5+OTJk8ycOZOHH36YtLQ09uzZw/DhwwF44IEHyM3N5cSJE0RHR7Ns2TLeeecd\n52t37NhBREQE6enpzJ8/n5/97GcAzJo1iyNHjnDs2DFn3Q8++IC77roLgN/85jfExcWxd+9e4uLi\nSEpK4plnnnHWPXPmDFlZWZw8eZIlS5bw4osvEhISQnp6Oqmpqfzxj3+sdNueffZZduzYwd69e4mJ\niWHHjh08++yz5dabm5tLcnIyb775JvPmzSM7O9v5fEREBDExMS79myvlNseOQV4eDBvm7pbUXiMP\nW2kqU9Xy8nXqXupq6tSppkOHDiYnJ+eK9fr162fWrl3rfPzcc8+Zm2+++bJ6RUVFpk2bNubw4cPO\nZW+88YaZMmWKMcaYpUuXmn79+jmfy8/PNyJiUlJSjDHGzJkzx/zhD38wxhhz5MgR07lzZ1NQUGCM\nMaZDhw4mPj7e+dqtW7easLAwY4wx0dHRpm3btubChQvO559++mlz0003mbi4uMva2bt3b7N+/Xpj\njDHh4eHmiy++cD63du3acutt3769KSoqcj7fvXt3s337dufjdevWmfDw8Er/3Vz5HijVKF591Zi7\n73Z3K5wc/zdqlLces+deH/FeF//4xz9ISEhg6tSpPP7441es6+vrS26ZU6gSExMJDw+/rF5aWhqX\nLl2iV69ezmWhoaEkJSU5H5ftt2/Xrh1gD+gC3HnnnXzwwQcAvP/++9x00020bduWs2fPkp+fz8iR\nI/Hz88PPz48ZM2Y4jxkAdOvWjdatWzsfP/7444SHhzN9+nT69u3LokWLKt225OTky9qbnJzsfOzv\n70+LFqVfq/bt2zvbC5Cbm0uXLl0qXbdSTYaH97eDdsu4JDU1lUcffZQ333yT119/nY8//pgtW7ZU\nWT8yMpIjR444H4eEhBAXF3dZva5du9K6dety4+okJCQQHBzsUrumTZvG2bNniYmJ4cMPP+TOO+90\nrrd9+/YcOHCAjIwMMjIyyMrKKtc9IlL+wHuHDh144YUXOHbsGKtWrWLx4sVsqOQsgR49elzW3h49\nerjUXoBDhw4xzJP/1FXe78IF2LgRpk1zd0vqRMPdBQ888AA333wzkyZNIjAwkEWLFnHfffc5DyRW\nNHPmTKKjo52P77rrLtavX8+KFSsoKioiIyODmJgYWrRowe23386TTz5JXl4eCQkJvPTSS/zkJz9x\nqV2tWrXitttuY/78+WRmZjLN8WUUEX7+85/zyCOPcPbsWQCSkpL48ssvq1zX6tWrnf33nTp1olWr\nVrRs2fKyerNnz+bZZ58lLS2NtLQ0/vCHP7jcXrAHXWfMmOFyfaUa3ZYtMHAgdO3q7pbUiYZ7NT75\n5BO2bt3K888/71x27733EhwcXO4AZVmzZs0iNjaWM2fOAHbP/fPPP+eFF17Az8+PESNGsHfvXgBe\nfvll2rdvT58+fZg0aRJz5szhpz/9aZXtqbjHPXv2bNavX8/tt99erjtk0aJF9O3bl7Fjx9KlSxem\nT59e7q+Jio4ePcrUqVPp1KkTEyZMYN68eUyaNOmy9/zd737HNddcQ2RkJMOGDeOaa67hySefdKm9\np0+f5tChQ84zepRqkrygSwb0IqYG8+abb3Lw4EEWL67yPuLNzmOPPUbfvn355S9/Wenz3vg9UB5o\n+HB47TUYP97dLXGqzUVMGu6qydDvgXK75GQYMgRSU6GSq6jdRa9QVUqpuvjyS7juuiYV7LWl4a6U\nUiXWrvXoq1LL0m4Z1WTo90C5VVERBATAnj12NMgmRLtllFKqtr79FoKCmlyw11a14S4ibUVku4js\nFpF9IrKginovi8hREdkjIsPrv6lKKdWA1q71ilMgS1Qb7saYQmCKMWYEMByYISKjy9YRkRlAuDGm\nH3A/oHdkUEp5Fg8fBbIil7pljDEltyJpC7QCKnaM3ggsc9TdDviISEB9NVIppRpUZibs3w/XXuvu\nltQbl8JdRFqIyG7gDLDOGLOzQpVgILHM4yTHMqWUavq++gomToSrrnJ3S+qNSydzGmOKgREi0hlY\nKSKDjDEHa/OGCxcudM5HRUURFRVVm9UopVT9aWJdMtHR0eXGp6qNGp8KKSJPAeeMMYvLLHsd2GCM\n+cjx+DAw2RiTUuG1zeJUyCeeeILAwEAeeuihRn3fMWPGsHTpUiIiIhr1feuLt30PlIcwxp4hEx0N\n/fq5uzWVapBTIUWkq4j4OObbAdOAwxWqrQLmOuqMBbIqBntzkZaWxrvvvsv9999fq9fn5uYSGxtb\nq9fOnz+fp556qlavVarZ+u47aNcO+vZ1d0vqlSvdMkHA30WkBfbH4CNjzOcicj/27iBLHI9nikgc\ncA6oelhDL7d06VJmzpxJ27Ztyy3/9ttvWbBgAXl5edx9991cvHiRmJgYbr/99nJdU8uXLy83JK6r\nrwM7GuX9999Pamoq3bt3b8jNVMp7LF0Kc+aA1GjHuOmr6a2b6lKow2323On//u//zC233FJu2YMP\nPmgeeeSRy+p+73vfM++9916l67nlllvMsmXLnI/3799vunfvXq7OwoULa/W6EtOmTStX15M09e+B\n8kIFBcb4+xtT5paUTRHefJs9d5ozZw5r164lJycHgKKiIj766CPuvvvuy+pWvDl2WRs3bmTChAnO\nx/Hx8XTq1Mn5ODY2loEDB9b4dWXpDaiVqoFVqyAyEsLC3N2SeucxQ5/J7+v+J5NZULuDdYGBgUya\nNImPP/6Ye++9lzVr1tCtWzeGD7/8QtysrKxKg3fv3r20bt2aPn36AFBQUMCSJUt49dVXnXVWrlzJ\nr3/96xq/rqxOnTo5bxKilKrGO+/AFW6O48k8JtxrG8z1Ze7cubz++uvce++9vPfee1XeWq7izbFL\nbNiwgV69erF8+XIuXLhAXl4er7zyCqGhoQAUFxdz6dIl2rRp49LrOnfuzKJFiwgICGDo0KGMHDkS\n0BtQK+WypCT45htYscLdLWkYNe3HqUvBQ/vcjTGmoKDA+Pn5mf3795uOHTuaxMTESutNnTrVvP/+\n+5ctv/HGG82rr75a5fq/+OILs337dpdft3jxYrN9+3Zz8eJFc+eddzqXa5+7Ui567jlj7rvP3a1w\nCdrn3nDatm3LLbfcwp133smYMWPoWcXIcRVvjg12r3zTpk1cf4WLJLZt28bo0eWG7Lni6+Lj4wkK\nCqJVq1ZkZmYCUFhYyK5du5w3ylZKVcEYr+6SAR3yt0buvvtu9u3bx9y5c6usM3fuXNasWUNhYSEA\nMTExPPHEExQWFlZ5xVl2djZ+fn7lllX3uuLiYlq2bAmU3oR61apVTJkyhcDAwFpsnVLNyLZt9tTH\ncePc3ZIGozfrqIHExEQiIiI4c+YMHTt2rLLe7373O7p37+7yFap/+9vfmDVrVo1C+a9//SsTJ04k\nIiKCe+65hw8++IBx48bx1ltvMWjQIJfX05R4yvdAeYFf/AL69IHf/tbdLXGJ3iC7ARUXF/Poo4+S\nl5fHm2++Wa/rfuaZZ3j66adr9Jr09HTefvttunTpwpAhQxjnBXsgnvA9UF4gP98ON7BvHwR7xviG\ntQl3jzlbxp3y8/MJCAggLCyMNWvW1Ou64+PjiYyMrPHr/P39mT9/fr22Ralm4V//gjFjPCbYa0v3\n3FWTod8D1Siuuw7uvx9uv93dLXGZdssoj6bfA9XgTpyAa66BU6c8aux2vUG2UkpdybJlcMcdHhXs\ntaV97kqp5qG42I4A+fHH7m5Jo9A9d6VU8/D119CxI1x9tbtb0ig03JVSzcNzz8Gjj3rfuO1V0HBX\nSnm/rVshPh7uusvdLWk0Gu5KKe/37LP2atTWrd3dkkajB1SVUt5t1y7Yuxf+/W93t6RR6Z67Usq7\nPfsszJ8PFe5r7O003F105MgRRowYgY+PD6+88kqldZ544glefvnlBm/LmDFjOHToUIO/j1Ieb98+\ne0OOn//c3S1pdHqFqovuu+8+fHx8ePHFFyt9Pi0tjREjRhAXF0dbF/cQcnNzSU5OrvKeq1VZsWIF\nH374ISu87A4ynvA9UB7mjjvsqY+PP+7ultSJXqHagBISEhg8eHCVzy9dupSZM2e6HOwAy5cvr/JG\n11cya9YsNmzYQGpqao1fq1SzERsL69fDr37l7pa4hYa7C6677jo2bNjAvHnz6Ny5M3FxcZfVWbNm\nDZMnT67Rek+dOkWPHj1q3J62bdsycuRI1q5dW+PXKtVs/PGP8NBDUIsdKG9QbbiLSE8R+VpEDojI\nPhG57A4UIjJZRLJE5DtH+V3DNNc91q9fz8SJE3n11VfJycmhb9++l9XZt29fjbpXYmNjGThwYK3b\nFBERQUxMTK1fr5RXi4+Hzz6DBx90d0vcxpU990vAo8aYwcA4YJ6IVJZKm4wxVzvKs/XaSrBXldW1\nNKCsrCxnF8upU6dYuXIltzuGFH3ppZcu66tfuXIlP/rRj1yuX1GnTp3Iysqq781QyjssWmS7Y7p0\ncXdL3KbacDfGnDHG7HHM5wGHgMpGuW/Y9DSm7qUB+fr6kpubC8DRo0cZMWIEKSkpANxxxx20aFH6\nT11cXMylS5do06bNFetnZmayaNEili5dyq5du8q9X25uLl2a8RdXqSolJtrBwR55xN0tcasa9bmL\nSG9gOLC9kqfHicgeEVktIp55E886iIyM5MiRIwBMmTKFt99+m7vvvhuAjIwMpk6d6qy7bt06pk2b\n5nxcVf2lS5cyZcoU5syZw+LFi8u936FDhxg2bFhDb5ZSnufpp+09Urt2dXdL3MrlK1RFpCOwAnjY\nsQdf1i6glzEmX0RmACuB/pWtZ+HChc75qKgooqKiatjkpmnmzJlER0cze/ZsAL799lt+/OMfA7Bn\nzx7uKjOmxbZt21iwYEG511dWf8mSJdx66620atWKzMxMZ93CwkJ27drFsmXLGnqzlPIsW7bAunXg\n4deBREdHEx0dXad1uBTuItIKG+zvGmM+qfh82bA3xqwRkddExM8Yk1Gxbtlw9yRSTZ/93LlzGTFi\nBIWFhbRt25bZs2ezcuVKYmNjGT9+vLNednY2fn5+l72+svrFxcW0bNnysrqrVq1iypQpBAYG1nGr\nlPIily7BvHnwwgsef4ZMxR3f3//+9zVeh6t77m8DB40xf6nsSREJMMakOOZHYy+OuizYPdnXX399\nxef9/f2ZO3cub7zxBg899BBz5syptN7y5cu57bbbLlteWf2BAweSmpqKv78/Pj4+zuWLFy/mrbfe\nquEWKOXlXnsN/P3B8Rdwc1ftFaoiMgHYBOwDjKM8AYQCxhizRETmAb8CLgLngV8bYy7rl/fkK1Tr\nyzPPPMPTTz/tUt309HTefvttunTpwpAhQxg3blwDt869mtP3QNWzM2dg6FDYtAkiItzdmnqnN8hu\n4uLj49m7dy833XSTu5vSJDWX74FqAHPnQlCQPQXSC2m4K4+m3wNVK5s2wZw5cPCgvY2eF9KxZZRS\nzcvFi/Yg6uLFXhvstaXhrpTyXH/9q+2OueUWd7ekydFuGdVk6PdA1UhyMkRG2vuj9q/0shqvod0y\nSqnmwRjbHXP//V4f7LWl91BVSnme11+HhAT48EN3t6TJahLhHhoaWu0VoMr7hYaGursJyhPs22fH\nj9mypdndF7UmmkSfu1JKuSQ/H0aPhsceg3vucXdrGo3HnueulFIu+eUvITcX/vGPBr9HQ1NSm3Bv\nEt0ySilVrX/+0474uHt3swr22tI9d6VU03fyJIwaBZ9+artlmhk9FVIp5X0uXYI774T//u9mGey1\npeGulGrafv97aNfOHkRVLtM+d6VU0/XBB7BsGWzfDi10X7QmNNyVUk3T5s3w8MOwfj3oXcdqTH8K\nlVJNz9GjcOut8O679iYcqsY03JVSTUt6OvzgB/DMM/D977u7NR5LT4VUSjUdhYUwbRqMHQvPP+/u\n1jQZeoWqUspzGQM/+QkUFMDy5XoAtQy9QlUp5bkWLIC4ONiwQYO9Hmi4K6Xc709/ssP3/uc/9px2\nVWf686iUcq8//QnefhuioyEgwN2t8Rq6566Ucp/nnoN33rHB3qOHu1vjVardcxeRniLytYgcEJF9\nIvJQFfVeFpGjIrJHRIbXf1OVUl5Fg71BubLnfgl41BizR0Q6ArtE5EtjzOGSCiIyAwg3xvQTkTHA\n68DYhmmyUsrjabA3uGr33I0xZ4wxexzzecAhILhCtRuBZY462wEfEdHOM6XU5TTYG0WN+txFpDcw\nHNhe4algILHM4yTHspQ6tE0p5U2Ki+3IjmvWaLA3ApfD3dElswJ42LEHXysLFy50zkdFRREVFVXb\nVSmlPEV+vr1AKT0dtm4FX193t6hJi46OJjo6uk7rcOkKVRFpBXwGrDHG/KWS518HNhhjPnI8PgxM\nNsakVKinV6gq1dykpsIPfwjh4faUx7Zt3d0ij9OQd2J6GzhYWbA7rALmOhoxFsiqGOxKqWboyBEY\nPx6mTrU3tdZgbzTV7rmLyARgE7APMI7yBBAKGGPMEke9V4DrgXPAT40x31WyLt1zV6q52LzZDtv7\n7LNw333ubo1H04HDlFJNw9tvw29/a8dj12F760wHDlNKuVd+PjzwAGzbZs+IGTTI3S1qtnRsGaVU\n/Th6FMaNs0P27tihwe5mGu5Kqbr717/sgdP774f33oOOHd3domZPu2WUUrV38SL8z//AihWwejWM\nHu3uFikHDXelVO0cPAhz59phenftAn9/d7dIlaHdMkqpmikqghdfhEmT4Oc/h88+02BvgnTPXSnl\nuvh4uOceO79jB/Tp49bmqKrpnrtSqnrGwBtv2D71G2+09znVYG/SdM9dKXVlsbEwbx5kZcGmTXqK\no4fQPXelVOXOn4ennoIJE+CGG+yFSRrsHkP33JVSl/v8c3ul6TXXQEwMBFe8P49q6jTclVKlEhPh\nkUdg7174f/9Px4XxYNoto5SCnBx48kkYPhyGDoV9+zTYPZyGu1LN2cWL8Npr0L8/JCXBnj2wcCFc\ndZW7W6bqSLtllGqOjIFVq+A3v4GQEPjiC7vXrryGhrtSzc2GDfD00/bUxj//2Xa/SI2GClceQMNd\nqeYiOtp2uSQn21Mc77wTWrZ0d6tUA9FwV8rbbdxoQ/3UqdJQb6X/9b2dfsJKeSNjYO1aWLTInt74\n1FNw110a6s2IftJKeZMLF+DDD+GFF+zjxx7TPfVmSj9xpbxBTg787W/2AGn//vD883qgtJnTcFfK\nkx09as9TX7YMpk+HlSth5Eh3t0o1AXoRk1KepqjI3iDj+uvtoF5XXQXffQcffKDBrpyq3XMXkbeA\nG4AUY0xkJc9PBj4B4h2L/mWMebZeW6mUgrQ0WLrU7ql37WoH9lq5Uq8mVZVypVvmHeCvwLIr1Nlk\njPlh/TRJKeVUVATr1sFbb9npjTfaA6Z6I2pVjWrD3RizWURCq6mmR22Uqk8nTsA779gSEAD33gtv\nvgk+Pu5umfIQ9XVAdZyI7AGSgPnGmIP1tF6lmo/sbPjnP+Ef/7BD7t55J3z6KQwb5u6WKQ9UH+G+\nC+hljMkXkRnASqB/VZUXLlzonI+KiiIqKqoemqCUh7pwAdassYH+5Zdw3XW2L33mTO1Lb8aio6OJ\njo6u0zrEGFN9Jdst82llB1QrqXscGGmMyajkOePK+ynl1S5dsoN3ffwx/OtfMHgwzJkDt94Kvr7u\nbp1qgkQEY0yNur9d3XMXquhXF5EAY0yKY3409gfjsmBXqlm7eLE00FeuhLAwuO022LULQqs7pKVU\nzblyKuT7QBTgLyIngQVAG8AYY5YAt4rIr4CLwHngxw3XXKU8yPnz8NVXNsw/+QTCw22g79wJvXu7\nu3XKy7nULVNvb6bdMsrbpaXZC4w++QS+/hquvtqevvijH+keuqq12nTLaLgrVRfGwP798PnnsHo1\nxMTA1Kk20H/wA/D3d3cLlRfQcFeqMeTlwfr1NtA//xxat7ZBPmMGTJkC7dq5u4XKy2i4K9UQiors\ngc9162zZtQvGjrWnK86caUdh1NEXVQPScFeqPhgDcXG2z3zdOjvt0QOmTbNl0iTo2NHdrVTNiIa7\nUrV1/Lg9VbGkgO1imTbN9qH36OHe9qlmTcNdKVcYA7Gx8J//2LJpExQWQlSUDfQpU6BvX+1qUU2G\nhrtSlblwAfbsgS1bbJhv3gzt28PEibZMmgQDBmiYqyZLw10pgDNn4JtvSsvu3fYCogkTbJhfey2E\nhLi7lUq5TMNdNT/nztm7EO3YUVqys+3ZLOPGwfjxduzzTp3c3VKlak3DXXm3ggI7FO6uXbbs3GnP\nahk61Ab4qFF22q8ftNA7SCrvoeGuvEdeng3yPXtKw/zIEXtO+ciRtoweDZGR0KaNu1urVIPScFee\nxxhITi4N8pKSmGiHwh0+3I7PMnKkDXId41w1QxruqmnLzYWDB22Q79tXOm3Z0natjBhhw3z4cBg4\nEFrV143ClPJsGu6qaTh3Dg4fhgMH7KBaBw7YkppqQzsy0oZ5ZKQtAQHubrFSTZqGu2pcaWk2xA8d\nKl9SU+1BzcGDy5c+fexeulKqRjTcVf0rKIBjx+wVnUeO2GlJKSqye+IREeVLWJiGuFL1SMNd1U5+\nPsTH2xCPi4OjR0tLaqq9a1D//vYqzpLSvz90765XdSrVCDTcVeWMgZQUOzhWfLydHjtWGuYZGTbA\nw8Ntd0rfvnbarx/06qUHNpVyMw335soY2/994gQkJNhpSTl+3JYOHWx3SZ8+dhoeXlqCg7UbRakm\nTMPdWxUWwqlT9tzvkyfLl4QEO73qKnuPzt69S0toqA3ysDC9/F4pD6bh7ony8yEpyZZTp0qniYm2\nnDoFmZl2PPGQEBvYvXqVltBQu7xzZ3dviVKqgWi4NyUFBXZ0wjNn7BWYyclw+nTpfHKyDfLz521w\n9+xpu0eCg+18SEjpNCBAx0pRqhnTcG9oRUWQnm4DOyWlNLxL5k+fLp3m5dlQDgqy4V2xBAXZ8Pbz\n0zNOlFJX1CDhLiJvATcAKcaYyCrqvAzMAM4B9xhj9lRRr2mFuzF2eNjUVDh71pay8ykp9nFJycgA\nHx8IDLQlIKD8fFCQnQ8KsqGte9tKqXrQUOF+LZAHLKss3EVkBvCAMeYHIjIG+IsxZmwV62q4cC8u\nhpwcu2ddUtLSys9XLOnp0K4ddOtmS/fupfPdutnADgiwy7t3h65d9bRApVSjq024V5tUxpjNIhJ6\nhSo3AsscdbeLiI+IBBhjUmrSEMeb2QOMmZmQlWWnFUtGRum0bMnKsqf7+fuXL1272unQoXa+YtHh\nYpVSXqg+dkODgcQyj5McyyoN99i5P8CnUOiUf4l2+RdokZNrg7mktG4NXbqAr6/tAvHzs/MlpW9f\nO/Xzs8Xf3067dLGvVUopVS/hXiO/3fMtOXKBbCmkVWAhncd2op1/IB0ChtElIJQAv14Edw4muFOw\nc+rXzg/Rg45KqWYiOjqa6OjoOq3DpbNlHN0yn1bR5/46sMEY85Hj8WFgcmXdMhX73ItNMWfPnSUp\nN4mknKRbFJ8OAAAYSUlEQVTy09wkTuWcIiknicKiwnJh37NzT+fjkvmgTkG0aqH94Uop79Ngp0KK\nSG9suA+t5LmZwDzHAdWxwJ/r+4DquQvnygV/Seg753OTSD2XSrf23ZyB37NTz9L5Mj8G7Vu3r/H7\nK6WUOzXU2TLvA1GAP7YffQHQBjDGmCWOOq8A12NPhfypMea7KtbVYGfLXCq+xOnc084fgVM5p5zB\nXzJNykmiQ5sOzrAvCf6SEtI5hJ6de9KprV6qr5RqOvQipmoYY0jLT3MGf8UfgFM5p0jMSaRVi1al\nod/JEfw+IeV+CHza+uhxAKVUo9BwrwfGGLILs0nMTiz3A5CYk1huCpTb2y+ZD/EJcS7zucrHzVuj\nlPIGGu6NxBhDTmFOaeBnJzrnS34AErMTaSEtnGFfNvhL/goI6RxChzYd3L05SqkmTsO9CSn7F0BJ\n2Cfm2HIy+6TzL4MObToQ0jmEXj69nMFfMt/Lpxc9OvWgdUs9f1+p5kzD3cMYYzibf7bKH4CT2SdJ\nyUuhe4fuNvB9QujVuRe9fHqVPvbphX87f+3/V8qLabh7oUvFl0jOTSYxO5GE7ATnD0BJ+J/MPknB\npQJn4If6hDrne/n0IrRLKD0796RNSx1mQSlPpeHeTOUW5pYL/ISsBE7mOKbZJzmdd5qu7bte9gMQ\n6hNKaBc73+WqLu7eDKVUFTTcVaWKiotIzk0uDf/sBBKyEuzUMd+yRUtn2If62NK7S2/n4+4dumvX\nj1JuouGuasUYQ8b5DBKy7Z7+iawTl4V//sV8evn0soFfJvh7d+lN7y69CewYSAvR8euVagga7qrB\n5F3IcwZ/SfifyHZMs06QVZDlDP+SEtYlzE59wwjoEKB7/krVkoa7cpv8i/nlwv945nFOZDumWSfI\nvZBLqE8oYb5h9PaxgR/WJcw51ZE/laqahrtqss5dOGdDP+u4M/yPZzlK5nGKTbEz6Pv49nEGfx/f\nPvTu0lsHfFPNmoa78liZ5zOdQX886zjxmfHOaUJWAr7tfOnj28eWLn1K5337ENQpSPv7lVfTcFde\nqdgUk5ybTHxmvA39zOMcyzzG8azjHMs4RnZhtnOPv49vH8J9w+3Uz06vanWVuzdBqTrRcFfN0rkL\n55xBfyzzGPGZ8RzLPMaxjGMkZCfQrX03wv3C6evbl3C/cMJ9w+nrZ+f1/H7lCTTclaqgqLiIxJxE\nZ/AfyzhGXGacnWbEcVWrq+jr17dc6efXj75+ffUgr2oyNNyVqgFjDKnnUjmacdQZ9nGZccRlxHE0\n/Sgi4gz6fn796Offz/nYv72/u5uvGlhRERQUQPv24O7feA13peqJMYb08+kcTT9qwz7jqC3pdtqq\nRSv6+/enn1+/8lP/fnRs09Hdzfd6ly5BXp4t587ZUjJfMs3Pv3xatpw/X37+/Hkb5iXTS5egXTvI\nyoLWbh6YVcNdqUZQMprn0fSjHEk/wpH0IxzNsPNxGXH4tvNlgP8A+vv3p79/f+d8mG9Ys76J+4UL\nkJ0NOTlXLrm55UtOTmmQl5QLF6BDB+jUyU47dICOHUvny5b27UunFUu7dqXTdu3gqqtKp23auH+P\nvYSGu1JuVmyKOZVziti0WGfwx6bHEpsey+nc04T5hjHAf4AtXQcwsOtABnYdiF87P3c3vVoFBZCZ\nafdkMzNL5yuW7Ozy8yXl4kXw8YHOnUunZUunTpdPy5aOHW3p1MmGb1MJ3sag4a5UE3b+4nniMuJs\n2KfZwD+cdpjDaYdp17qdDXr/gc7Aj+gWQS+fXvV6Dr8xds83Pd2WjAxbys6XlJIAL5kvLgZfX1u6\ndCmdlp338Sld5uNTvrRr17wCuT5puCvlgYwxnM477Qz6w2mHOZR2iENnD5FZkMkA/wFEdIsgoqst\ng7oNoq9fX1q3bE1+Ppw9a0taWvlSsiw9vXSakWH7j/39bfHzK51WLL6+pVNfXw1nd9JwV8rDXbxo\nQzklBVJTIeFMDgdSDhOXdYiTBYdILT5IdptDXGh7CjLDkbRBdDg/CL9Lg+jRajAhHfoR2K0NXbtC\n1642uMvO+/vbLg3lWTTclWqCiorsnvOZM5WXlJTSkp1tg7h7dwgIsNOypVs3O+3ke56sVkc4ce4g\nh84e5MDZAxw4e4CT2Sfp49uHwd0G29J9MEO6D6GvX99mfTDX0zVYuIvI9cCfgRbAW8aYRRWenwx8\nAsQ7Fv3LGPNsJevRcFdeo6jIBnJycvly+nT5cvas7YMODISgIDstKQEB5Yu/P7RsWfs2FVwqIDYt\nlv2p+zlw9gD7U/ezP3U/Z/LOMKDrAIZ2H8qQ7kMY0n0IkQGRBHcK1gu1PECDhLuItACOANcBycBO\n4A5jzOEydSYD/22M+WE169JwVx7h3Dk4dcqWpCRbys4nJ9tuE39/6NGjtAQFlS89etjQdvd50nkX\n8jh49iD7UvaxP3U/+1L3sTdlL5eKLxEZEMnQ7kOJDIgkMiCSId2H0KFNB/c2WJXTUOE+FlhgjJnh\nePxbwJTde3eE+2PGmFnVrEvDXbldQQEkJtpy8qSdnjpVflpQAMHB0LOnLcHBl5fAQPeHdl2l5KU4\ng76kHE47TIhPCMMChjEsYBjDA4czLHCY7uW7UUOF+y3A940xv3A8ngOMNsY8VKbOZOCfwCkgCZhv\njDlYybo03FWDMsaeFZKQYIM7IaG0lAR5VpYN7JAQ6NXLTkNCSoM8JMSeJdJcc+xi0UVi02OJORPD\nnjN7iEmx02JTzIigEYwItGV44HD6+/enZYs69CMpl9Qm3OvrCMsuoJcxJl9EZgArgf6VVVy4cKFz\nPioqiqioqHpqgmoOjLGn8x0/DidOVF7atLGhHRpaWsaPL13WvTu00OHfq9S6ZWtnv/xdkXcBpadr\n7j69m91ndrPi0Aqe/PpJUs+lEhkQycigkYzsMZKrg65mULdBevC2jqKjo4mOjq7TOlztlllojLne\n8fiybplKXnMcGGmMyaiwXPfcVbUKCmxIHzsG8fG2HD9eWlq0gLAwW3r3tqVkPjTUXt2oGkdWQRa7\nT+/mu9Pfsev0Lnad3sWpnFMM7T6UkUEjGRU8ilE9RjGw60Ddw6+DhuqWaQnEYg+ongZ2ALONMYfK\n1AkwxqQ45kcDy40xvStZl4a7Aux4IXFxpeXYsdKSmmr3ssPDoU8fW0rCPCzMXlCjmq6cwhz2nNnD\nt8nfsjN5JzuTdpJyLoWrg65mVA8b9mN6jiHUJ1T78F3U0KdC/oXSUyH/JCL3Y/fgl4jIPOBXwEXg\nPPBrY8z2Staj4d6M5OXB0aO2HDlSOh8XZ89G6du3tISHl5aePaGV/lXvVTLOZ9iwT9rJzuSdbDu1\nDYAxPccwNngsY3qOYVSPUXRq28nNLW2a9CIm1eguXbJdJbGxNsBjY0vns7JsWPfrB/3722lJCQxs\nvgcsle3DP5l9ku1J29l+ajvbkrYRcyaGcL9wxvccz4ReExgfMp6wLmG6d4+Gu2pAOTk2tA8ftuXQ\nITuNj7fncw8YYAN8wIDS+Z499cClct2FogvsObOHrYlb2Zq4lS2JWygqLmJ8yHiu7XUtE3tNZETQ\niGZ5sFbDXdVZRgYcPHh5ycy0gR0RAQMH2hIRYbtU2rVzd6uVNzLGkJiTyJaTW9h8cjObTm4iISuB\nsT3HMrHXRCaFTmJMzzHN4gboGu7KZdnZcOAA7N9ffnr+PAwaZIN70KDS+dBQ3QtX7peen86WxC38\nJ+E/bDq5iQOpBxgVPIrv9f4eU8KmMDp4NG1atnF3M+udhru6TGGh7T7Zt8+WvXttkGdm2uAeMgQG\nDy6dBgdrX7jyHDmFOWw+uZmvj3/NhhMbOJJ+hPEh47ku7Dqmh08nMiCyXsfDdxcN92bMGDvCYExM\nadm7155aGBYGQ4eWL71765648j6Z5zPZmLCR9fHr+TL+S7IKspjWZxrTw6czrc80gjoFubuJtaLh\n3kwUFdmDm3v2wO7ddhoTY++UM2xYaYmMtF0qOn63aq6OZx5nXfw6vjz2JeuPr6eXTy9u6HcDN/S/\ngdHBoz3mwioNdy9UWGi7U777zpbdu223So8eMGIEDB9uy7Bhdpl2qShVuUvFl9h+ajufHfmMz45+\nRkpeCjP7zeSG/jcwPXw6nds23UubNdw9XEGB7UrZtau0xMbaM1JGjoSrr7aBHhmpl9grVVcnsk6w\n+shqPjv6GVtObmFS6CRuHXQrPxzwwyZ3w3INdw9y8aI9O2XnTvj2W1sOHbKnG15zTWmYR0bqqYZK\nNbScwhxWH1nNikMr+Cr+K8b1HMdtg27jxoE30rV9V3c3T8O9qTLGXsW5fTvs2GHLnj329MJRo2yY\njxplu1Y0yJVyr7wLeXx+9HNWHFxB9IloTjxygvat27u1TRruTUR2tg3wbdts2b7dHtQcPdqWMWPs\nnrl2rSjVtF0sukjrlu6/I4uGuxsUF9t+8a1b4ZtvbJifOGG7VMaOtWXMGHv+uFJK1YaGeyPIz7f9\n5Fu2lAa6j4+9GcS4cTbMIyM9//ZrSqmmQ8O9AaSnw+bN8J//2Om+ffZqzgkTbBk/3g6cpZRSDUXD\nvR4kJcHGjbZs3mxvmDx2LEycCNdea/vM27v32IpSqpnRcK+FEydKw3zjRju07aRJtkycaLtY9MYR\nSil30nB3QWIibNhQWgoKYPLk0hIRoWOuKKWaFg33SqSmwtdf27Jhg707UFQUTJliy8CBesm+Uqpp\n03DH3rdz0yZYvx6++goSEmwXy3XX2TAfMkT3zJVSnqVZhntRkR1Qa+1a+PJLO3/NNTB1qg30UaO0\nz1wp5dmaTbgnJtog//JLu3ceGAjf/z5Mm2b30jt0qIfGKqVUE+G14V5YaE9LXLMGvvjC3pRi2jSY\nPt1Oe/ZsgMYqpVQT0WDhLiLXA38GWgBvGWMWVVLnZWAGcA64xxizp5I6Lof7yZOwerUN9I0b7Vks\nM2bA9dfbbpeWnjHGvlJK1Vltwr3aQ4si0gJ4Bfg+MBiYLSIDK9SZAYQbY/oB9wOv16QRAJcu2b3z\n//kfexu4kSPt5f2zZ9tbxW3bBgsW2HFammqwR0dHu7sJDUq3z3N587aB929fbbhy3sho4KgxJsEY\ncxH4ELixQp0bgWUAxpjtgI+IBFS34uxs+OgjuOsu22/+4IM2uJcssV0v775rw72r+4dTdom3f8F0\n+zyXN28beP/21YYr55EEA4llHp/CBv6V6iQ5lqVUXNmJE/Dpp7BqlR0Kd9IkmDULFi3SvnOllKov\njX6S4OjRcMMNMG8erFypZ7YopVRDqPaAqoiMBRYaY653PP4tYMoeVBWR14ENxpiPHI8PA5ONMSkV\n1tW0BpZRSikPUdMDqq7sue8E+opIKHAauAOYXaHOKmAe8JHjxyCrYrDXpnFKKaVqp9pwN8YUicgD\nwJeUngp5SETut0+bJcaYz0VkpojEYU+F/GnDNlsppdSVNOpFTEoppRpHow2hJSLXi8hhETkiIr9p\nrPdtLCJyQkRiRGS3iOxwd3vqSkTeEpEUEdlbZpmviHwpIrEislZEfNzZxtqqYtsWiMgpEfnOUa53\nZxvrQkR6isjXInJARPaJyEOO5d7y+VXcvgcdyz3+MxSRtiKy3ZEj+0RkgWN5jT+7Rtlzd1wIdQS4\nDkjG9uPfYYw53OBv3khEJB4YaYzJdHdb6oOIXAvkAcuMMZGOZYuAdGPM844faF9jzG/d2c7aqGLb\nFgC5xpjFbm1cPRCRQCDQGLNHRDoCu7DXovwU7/j8qtq+H+MFn6GItDfG5ItIS2AL8BBwCzX87Bpr\nz92VC6E8ndCIfwk1NGPMZqDiD9WNwN8d838HbmrURtWTKrYN7Gfo8YwxZ0qG/zDG5AGHgJ54z+dX\n2fYFO572+M/QGJPvmG2LPS5qqMVn11hhVNmFUMFV1PVUBlgnIjtF5OfubkwD6V5yFpQx5gzQ3c3t\nqW8PiMgeEXnTU7ssKhKR3sBwYBsQ4G2fX5nt2+5Y5PGfoYi0EJHdwBlgnTFmJ7X47LxmT7MJmGCM\nuRqYCcxz/Onv7bzpaPxrQB9jzHDsfyqP/tMewNFlsQJ42LGHW/Hz8ujPr5Lt84rP0BhTbIwZgf1r\na7SIDKYWn11jhXsS0KvM456OZV7DGHPaMT0L/JvLh2jwBiklYwY5+j1T3dyeemOMOVtmyNK/AaPc\n2Z66EpFW2OB71xjziWOx13x+lW2ft32GxpgcIBq4nlp8do0V7s4LoUSkDfZCqFWN9N4NTkTaO/Yi\nEJEOwHRgv3tbVS+E8n2Yq4B7HPN3A59UfIEHKbdtjv8wJW7G8z+/t4GDxpi/lFnmTZ/fZdvnDZ+h\niHQt6U4SkXbANOwxhRp/do12nrvjtKS/UHoh1J8a5Y0bgYiEYffWDfYAyHuevn0i8j4QBfhjB4Bb\nAKwEPgZCgATgdmNMlrvaWFtVbNsUbN9tMXACuL+yq6w9gYhMADYB+7DfSQM8AewAluP5n19V23cn\nHv4ZishQ7AHTFo7ykTHmf0XEjxp+dnoRk1JKeSE9oKqUUl5Iw10ppbyQhrtSSnkhDXellPJCGu5K\nKeWFNNyVUsoLabirKxKRxSVDxjoefyEiS8o8fkFEHqlmHZtdeJ/jjnN5Ky6fLCLjatjmNiKyzjHs\n6201eW2F9fiIyK9crJtb1/WLSJCILK+ibqCIfOqYv1tE/lqD9+kqImtq2j7l2TTcVXW2AOMBRESA\nrsDgMs+PB7ZeaQXGGFfG2anqgouokvevgavt25qrjTEf1/C1ZfkC/+Vi3dpcMFJu/caY08aY26uo\n+yiwpMxjl9/PGJMGJNf0R1J5Ng13VZ2tlIbrYOwl3bmOvc42wEDgOwAReUxEdjhG5VtQsoKSvVqx\nXhORg44bDqwWkZtLqgEPicgusTc96S/2vr2/BB5x7IVPKNswxw0M/u2ov1VEhohIN+BdYJTjNWEV\nXhPu2KvfIyLfikiYiHQQka8cj2NEZJaj+nNAH8d6FlVS74cV/7FE5O9ll4vIP0RklogMEnsThu8c\n7x3uWH94mfWHisi+Kj6HW4AvKnm/H4jIFhHxE5E+IvKNo21/qPDXxCfAnCrWrbyRMUaLlisW4Bh2\nsLdfOMrvsYMZjQc2OupMA95wzAvwKXCt43GOY3or8JljPgDIAG52PD4O/Jdj/lfAEsf8AuDRKtr1\nMvCUY34KsNsxPxlYVcVrtgE/dMy3Aa7C7uR0dCzzx957ACAU2FvmtZXWq7CNk4B/O+Y7O/7tWjja\nOtuxvBV2rO6K6y/3uMzy3sDOMo/vdqzvJmAj0Nmx/FPsZekA95e0yfG4R2Xr1uK9RffclSu2AhOw\nYf4NNiBLHm9x1JkOTBOR77B78gOAfhXWMwE7Ng3GjvmxocLz/3ZMd2EDrTrXYvfSMcZsAPxKBnCr\njOO5HsaYVY7XXDDGFGDD9zkRiQG+AnqISGXjZV+pnjjWuQk7SJ4/MBv4pzGmGPvv9qSIPA70NsYU\nurB9JYKAsxWWXQc8DvzA2NEDAcZhR0oEeL9C/VTHelQz0crdDVAeoaRrZgi2W+YU8N9ANvCOo44A\nzxlj/laH9ykJvCJc+25W7Heu7V147sIeSxhhjCkWkePYPfqa1CvblmXAT7Cjn94DYIz5QES2ATcA\nn4vIL7B/rbjifCXtOQaEYX9Ed1XShor/Flc51qOaCd1zV67Yig2lDGNlAl2we4olB1PXAj9zDHmM\niPQQka6O50qCZgtwi6PvPQB7sLQ6udjujcr8B0c/sohEAWeNvWlDpRzPnRKRGx2vaSN2WFUfINUR\n2FOw3SMl792pzCqqqld2G8GO6veIfUt7n2ARCTPGHDfG/BXb/x1ZyforrqfEEWyQl3UC2w+/TEQi\nHMu2Ybu+wP6wlNUfDxwCV9WehrtyxT5sH/M3FZZlGWMyAIwx67BdAd+IyF5s90tJcJXsUf4Tu9d/\nALt3uwu791+2TkWfAj+q7IAqtu9/pKOb5I/Yvujq/AR74DYG+2MTALyHPQAbg/2xOOTYpgxgi4js\nFXtz8H9UVq9i+40xqY7n3inz/O0isl/s7dMGY2/OXXH9lf47GHtPzTgR6VNh+RHsXxMfOw4c/xp4\nVET2AOGU/tuCPSax2oV/H+UldMhf1ahEpIMx5pzYc9q3Y29P6LF3BKqMiLQHYoCrjTE1Pv+9inXe\nCIw0xjx9hTrtjDHnHfM/Bu4wxvzI8TgauNEYk13V65V30T531dg+E5EuQGvgGS8M9uuAt4AX6yvY\nAYwxnzgO0l7JSBF5Bdu1kwn8zNGmrsBiDfbmRffclVLKC2mfu1JKeSENd6WU8kIa7kop5YU03JVS\nygtpuCullBfScFdKKS/0/wGjSTya4nwEHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105d8add0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "epsilon = -0.15\n",
    "f = (1+epsilon*X)/y # volumetric flow rate ratio\n",
    "\n",
    "plt.plot(W_output, X, label='X (conversion)')\n",
    "plt.plot(W_output, y, label='y ($P/P_0$)')\n",
    "plt.plot(W_output, f, label='f ($v/v_0$)')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel('Weight of catalysti (kg)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integrating other variables\n",
    "Sometimes it's convenient to integrate things like $N_{tot}$ or Volume as well as the variables that you think you care most about. Just add them to your vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}