{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerically solving differential equations with python\n",
    "\n",
    "*This is a brief description of what numerical integration is and a practical tutorial on how to do it in Python.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Software required\n",
    "\n",
    "*In order to run this notebook in your own computer, you need to install the following software:*\n",
    "\n",
    "* [python](http://python.org)\n",
    "* [numpy](http://numpy.org) and [scipy](http://scipy.org) - python scientific libraries\n",
    "* [matplotlib](http://matplotlib.org) - a library for plotting\n",
    "* the [ipython notebook](http://ipython.org/notebook.html) (now renamed to [*Jupyter*](https://jupyter.org/))\n",
    "\n",
    "On Windows and Mac, we recommend installing the [Anaconda distribution](https://store.continuum.io/cshop/anaconda/), which includes all of the above in a single package (among several other libraries), available at http://continuum.io/downloads.\n",
    "\n",
    "On Linux, you can install everything using your distribution's prefered way, e.g.:\n",
    "\n",
    "* Debian/Ubuntu: `sudo apt-get install python-numpy python-scipy python-matplotlib python-ipython-notebook`\n",
    "* Fedora: `sudo yum install python-numpy python-scipy python-matplotlib python-ipython-notebook`\n",
    "* Arch: `sudo pacman -S python-numpy python-scipy python-matplotlib jupyter\n",
    "\n",
    "Code snippets shown here can also be copied into a pure text file with .py extension and ran outside the notebook (e.g., in an python or ipython shell).\n",
    "\n",
    "### From the web\n",
    "Alternatively, you can use a service that runs notebooks on the cloud, e.g. [SageMathCloud](https://cloud.sagemath.com/) or [wakari](https://www.wakari.io/). It is possible to visualize publicly-available notebooks using http://nbviewer.ipython.org, but no computation can be performed (it just shows saved pre-calculated results)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How numerical integration works\n",
    "\n",
    "Let's say we have a differential equation that we don't know how (or don't want) to derive its (analytical) solution. We can still find out what the solutions are through **numerical integration**. So, how dows that work?\n",
    "\n",
    "The idea is to approximate the solution at successive small time intervals, extrapolating the value of the derivative over each interval. For example, let's take the differential equation\n",
    "\n",
    "$$ \\frac{dx}{dt} = f(x) = x (1 - x) $$\n",
    "\n",
    "with an initial value $x_0 = 0.1$ at an initial time $t=0$ (that is, $x(0) = 0.1$). At $t=0$, the derivative $\\frac{dx}{dt}$ values $f(0.1) = 0.1 \\times (1-0.1) = 0.09$. We pick a small interval step, say, $\\Delta t = 0.5$, and assume that that value of the derivative is a good approximation over the whole interval from $t=0$ up to $t=0.5$. This means that in this time $x$ is going to increase by $\\frac{dx}{dt} \\times \\Delta t = 0.09 \\times 0.5 = 0.045$. So our approximate solution for $x$ at $t=0.5$ is $x(0) + 0.045 = 0.145$. We can then use this value of $x(0.5)$ to calculate the next point in time, $t=1$. We calculate the derivative at each step, multiply by the time step and add to the previous value of the solution, as in the table below:\n",
    "\n",
    "| $t$ | $x$      | $\\frac{dx}{dt}$  |\n",
    "| ---:|---------:|----------:|\n",
    "| 0   | 0.1      |  0.09     |\n",
    "| 0.5 | 0.145    |  0.123975 |\n",
    "| 1.0 | 0.206987 |  0.164144 |\n",
    "| 1.5 | 0.289059 |  0.205504 |\n",
    "| 2.0 | 0.391811 |  0.238295 |\n",
    "\n",
    "Of course, this is terribly tedious to do by hand, so we can write a simple program to do it and plot the solution. Below we compare it to the known analytical solution of this differential equation (the *logistic equation*). **Don't worry about the code just yet**: there are better and simpler ways to do it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f41d70014a8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEFCAYAAAD+A2xwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4lFXexvHvmTQILRB6SSAgVUUI\nSsAFxIKIjcVXV8WCCIgFXBVxdS2gwuKqi2tZkUVQsVdEEQTpRUpCxxBKQkInhAQSQuo87x8T2AAB\nZiBPJjO5P9eVa9o5md+E5ObMmeecx1iWhYiI+C6HtwsQEZELoyAXEfFxCnIRER+nIBcR8XEKchER\nH6cgFxHxcQpyEREfpyAXEfFxCnIRER8X6O0Cihs0aJDVuHFjb5chIuJ1o0eP/tCyrEHutC1XQd64\ncWNGjRrl7TJERLxu9OjRu9xtq6kVEREfpyAXEfFx5WpqpSQZGRns3bvX22VUSA0aNCAsLMzbZYjI\nOZT7ID948CBNmzalcuXK3i6lQjl27Bi7d+9WkIv4gHI/tZKfn0+lSpW8XUaFU6lSJfLz871dhoi4\nodwHOYAxxtslVDj6mYv4Dp8IcnfEJafz3vxtxCWne7uUs7ryyiuJi4vzuF9CQgIAI0eO5O233y7t\nssTLvPX7643nrUivtazYMkdujKkPTATetCxroTGmHjAMCABesSwruzSfLy45nf6TlpNX4CQ40MFn\ng2KIjqxZmk9RaqZMmUJUVJRHfbZs2cKwYcOYPXs2TzzxBCEhITZVJ94wc8Nehn+5hoJCiwCH4dbL\nGlK3eiWcTgunZeG0oNBpYRVdd91n4XRy4vET9x2/XqzvSdeL9T18LI/N+zJxWuAw0KpeNapVCrLn\nRVoWgeSTm5NN8oEMHFYhX5hCLqpdmeohEGAVEkDRl1X8soAACnFYhQSeuCzAUdTGgRMHToxluS6x\nMDhxYGGs/13PKyjgUFYODpwsnw8NL2tIg+rBYDmLfVlFl4Wn3F/8cQsodln02k6679TLq1+A2hfZ\n83MtYkuQW5a1zxhTHTj+/rwvUAA0BroBvx5va4wZAgwBuOmmm876fZv+bcY5nzsn38lt7y87a5sd\n424s8f4VK1bw22+/sXjxYrp27cpLL73EuHHj+OKLL5g7dy5PP/00a9eupWPHjlStWpVx48YRHR1N\nTEwMv//+OzNnzuS5556jXbt21K1bl969e9OlSxeGDBlCVlYWjz76KF27dmXSpEls2LCBDz74gIED\nBzJv3jwefPBBxo4dy5o1a/jiiy9wOp1MmjSJ77//nkmTJhEfH8+MGTN46623aNeuHSNGjOC5556j\nYcOGdO/enUsuuYSIiAjGjBnDhx9+yPz584mIiDjnz0vKlmVZpBzKZnliGisSD7Ei6RDp2XnkF7pC\nodBpceBILlF1quIwhgAHOIzBGIPDQIDjf9cdxhBgDOb4dcf/rh/va0q47jAQYAw/rtvD5n2ZBFJA\nGNn0rOugd6tqBORnE1CQjaMgm4D8bBwFR4suj99/lICCYzgKcjDOPByFuZjCXByFRdeL7nPdX3Td\nmYfTEUwegWQHOyggkAICCMkJJpQQLEcQliMAyxEIjkAsE3DSfZYJdF03x+8LwjIBYBxYxlHsugEc\nRdcDgEAs4yBh/1EOZxXixJDndJB8pJAGtauAcZT85QgAY/53m+PXTdH14peUfN/x25Vr2f57VVZH\nrQQXu37S5KtlWRNxjd4ZNWrUWc8EfaYAPj4izy9wEnQBI/KsrCwaN25Mfn4+a9asAeChhx5izpw5\nfP3110RERHDo0CFGjRp1IjTr1KlDjRo1WLVqFVOnTiUjI4M777yTDh06MHDgQF5++WWefPJJkpKS\nCAsLo2rVqgA0adKE2rVr069fP1555RWWLFnCwIEDSUxMJC8vj6ZNm5KSksLq1atp27Yt4eHh3Hjj\njXzzzTcAvPvuu1SrVo1evXoxfPhwNm7cCECfPn2YOHEia9euVZCXA5ZlkXjwaFFou8LbaVnERIXT\nOaoWj/RsweHsPPp/uOLE7+9fr2t5fu8oLQtyM+Foqusr6wAcPQBH0yAnA3IOw7EM1/VjGXTKOsSL\nwYcIIZ8sQqm8pyYhGdUhOBSCq0Bw1aLLKhASCtWqQnCd/90XWKnoK9h1GRACgcW+TtyuBAHBOBwO\nNp36t3pP2bx73peczn+KPW/3njFQTt+1nw+7plbCgJpAhDFmLnAf8ASQASwq7eeLjqzJZ4NiWJ6Y\nRkxU+Hn/YjgcDhYuXEijRo04cuTIifuDgoJO+vDv+O3j90VERFClShWcTifGGAoLC8nLy8OyLHr1\n6kVQUBDTpk1jwIAB56zB6XSSk5PD9u3bqVy5MpZllfjB46nPVdLjUvYsy2LbgSyWJx1iRWIaK5IO\nEegwxESF0yUqnL9e25Km4aGn/Zue8/e3sAAy98KR3XB4FxzeCYd3u25n7oOjB12hbQKgah2oUheq\n1oUqtSG0NlRvCHXbQuUwqBQGlWoQXDmMTamwbGcuMc1rl0mgltbfqq88b1mxa2olA2hfdPOTossR\ndjzXcdGRNS/4H2fVqlVs3ryZNm3anPhDmzp1KkeOHOG2227jvffeY+fOnbz55psMHz4ch8NBamoq\nO3fuBODuu+9m4cKFTJs2jddee43U1FTeeusthgwZwksvvUSvXr3IysoiJSWFrKwsMjIySExMBGDl\nypUA7N27l3nz5lG3bl0aNWpEYmIiPXv25NixY8TFxbF7924yMzN54403eOGFF1iwYAFjx45l1y7X\ntgxbtmwhKyvrxG2xl9NpseVA5okR98qkQ1QKCqBzs3C6t6zDyOtb06RW5XMeBRQdWZPohpUhfQfE\nL4W0bXBoO6Rth/RkyNrvCuUajV1f1RtBeAuI6gFV60GVOq7gDq7iUf0dwqCDvdO3pymNv1Vfet6y\nYCzrrLMZZWrUqFHWqZtmxcfH06ZNmzKvZceOHTRr1oz09PQTi2JGjRrF2rVrmTZtWpnX4w3e+tmX\nZ06nRfy+IycFd/XKQXRuVovOzVzTJY1rhp79mxQWuIJ6/0bYv8n1lRoPmfshLALCm7tCulaU67Jm\nJFRr6JrCkArDGDPasqxR7rQt9ys7vSU5ORmAxMREOnbseOK+ffv2kZ2dTWjoOf5YxWfFJaefeAt+\nWZMw/thzhBVJaSxPPMSqHYcIrxJM56ha3HBxA0bd0o4GNc6y6tjphIMJsGuV62vvOkjdAtUbQL12\nUO8S6Hgf1G0DYZEQoD9J8Zx+a86gR48enPpuZcqUKV6qRspKXHI6d//XdSirMVApKICGYZXp3KwW\nt1zWkLF/vpi61c+y0jj/GOxcATuWuIJ792oIDYfGl7u+OtzrmqsOqVp2L0r8noJcpJjPVyaTW+D6\noNgAA7o2ZWTv1mfuUFgAe9ZA0gJIWuQK7rptoGk36PwwNO7kmtsWsZGCXKTI+l0ZzN60n+BAB4WF\nrsPUrmlT7/SGuZmwbS4kzISts6FaA4i6CmIehciuUKl6WZcuFZyCXARITjvKoI9jefP29oRXDTn9\nMLVj6fDHdIifDikroMkV0OoGuOYF11EkIl7kN3uteMPChQsxxpCRkVHi48f3R3F3f5UHHniAvn37\nuvXc2dnZtGvXjrfeeuuMbRISEnA6nURFRbFv3z63vm9FlJaVy/2TVzL8movo1a4+0ZE1ebRnC9fh\ngJumwZf94a1LYdtv0OEeePIPuPd7uGKwQlzKBY3IL0BkZOQZH/v+++9ZsmQJ//rXv9zeXyUyMpL0\ndPc29AkNDaVOnTpnfHz8+PEYY2jVqhU///wz9evXd+v7VjTZeQUM/DiWmy5tyD0xRf+e+zdB7GTY\n8C00uBQuuQNufc+1mEakHNKIvMiKFSsYM2YMvXv35pdffmHnzp0YYxg7dizNmzcnJSXltDbHFRQU\n0KRJE8aMGcOiRYv4+OOPGT9+POvWrWP16tV07dqVn3/+mfj4eJ5++ml69OhBUlISL7/8MlOmTKFL\nly5kZmaeVlNKSgqDBw+mb9++FBQUMGzYMCZOnMgjjzxy0hE1r776KpdddhmfffYZYWFhZGRkMGHC\nBJYtW8bChQtp164da9eu5dNPP2XUqFHce++9pKen88orr9CkSRNGjhzJnXfeWSY/5/KkoNDJsM/X\n0KJOVZ66OgLWfQUfXg+f3uZaDfnwUrj/J+h4r0JcyjXfGpGPqlFK3+fwaXcV32dlzpw59OnTBzh5\n75IqVaqc1KZt27YABAYG8thjjzF37lw6d+7Mvffey/z58wkLCzuxwRa4Ardbt2706tWL0FDXMu2I\niAiWL19+YoVncRkZGXz55ZeMGzeOhIQEJk+ezNGjR4mIiGDgwIEn2jVu7Hp736hRIwDCwsJo0KAB\nXbt2pUePHifajRgxgp9++ol33nmHDz/8kCZNmhAeHk6fPn3o169f6fxsfYRlWbzw40aCCo7wz3pr\nMP++3XVcd9dh0LK3jucWn+Jbv60lBHBpKb7PSkmrXZ1O51nbDBw4kFGjRnHPPffgcDhKXJKdk5ND\nbGwsQ4cOxel0kpaWdmJPl5Kes2HDhnzwwQc89NBDLF269MT3PHV/FYfj9DdWnu7PUtH2Zpn8y1Ku\n2PJf+rIQE34D3DsN6rX1dlki50VTK0WO77MSFBREcnLyif1Tiu9dcmqb46s/U1JSqFOnDn379qV9\ne9cWMzExMcTHx5OcnHxif5WRI0cye/Zsnn32WdLS0pg7dy7Lli2jZcuWbN++neTkZPbv309ubi4A\nGzduZN26ddxxxx1ERkbSv39/Jk+ezLXXXkv79u1JTU0lJSWF6OhosrOzmTNnDrm5uRw8eJArr7yS\nlStXnlTjuHHjTkwJPfjgg+zatYuMjAySkpJO9PN7WQfY8tEj/N+qv3B9u/qYh5fAn99XiItP014r\npWTv3r1s2LCBXr16ebuUUuMrP3u3HEuHpW+Tv3Iy3+RfSdcBY2ga2czbVYmckSd7rWhEXgqysrK4\n+eabqVGjlObwpfQUFsCKifBOJ9IO7OHWgn/QasB7CnHxK741R15OVa1aldjYWG+XIadKWgQzn4Eq\ntdnT92v6fpPOq7dd7LdbmUrF5RNBfqaTK4h9ytOUm8eyDrgCfFcsXP8qB5tcz90Tfufxa10LfkT8\nTbmfWgkKCiInJ8fbZVQ4OTk5BAXZdCJeu1gWrPsS3u/q2tf7sZVkt7iRB4sW/PTvfOYFXCK+rNyP\nyGvXrs2OHTu8XUaF1KBBA2+X4L7Du+Cnx10nZ+j/DTTsQEGhk8emxtGibjWe6tXS2xWK2KbcB3lY\nWNiJM/SIlGjTDzBjBHR+CP70BAQEYVkWz0/bSIHTYtxtl2hqTvxauQ9ykTPKOwozR8KOpdD/a2gU\nfeKhf8/dysY9h/lqSBeCAsr9DKLIBdFvuPimfRvgg+6uU6kNXXxSiH+1KoXvVu9i8oDLqRKisYr4\nP/2Wi+/Z8K1rJN57HFx6x0kPzdu8n9d/3cLXD8VQt9pZTskm4kcU5OI7Cgvgt5dg889w349Q/5KT\nHl63M4MR36xn0v2diKqjc2JKxaEgF9+QfQi+uR8cgTB4PoTWOunhHQePMuiTWF677VI6RmjBj1Qs\nmiOX8i99B3zYCxq0h/7fnhbiB7NyuX/KSv567UVc17aEc2yK+DmNyKV8270avrgLuo9wnVrtFNl5\nBTz40Spuaa8FP1JxKcil/EqYBT8+Are8A61vPO3hgkInj362movqVePJ67TgRyouBbmUT+u+gjkv\nwN1fQ+NOpz1sWRZ//2EjTgv+0U8LfqRiU5BL+RP3ESx4De6bDnVbl9jk33O38sfeI3w5JEYLfqTC\nU5BL+bJ8Avz+Hgz4GcKbl9jky5UpfL96N9893FULfkRQkEt5suQt12j8gRmu3QtLMG/zft6Y7Vrw\nU6daSNnWJ1JOKcilfPj9P7D6Y3jgF6jesMQma7XgR6RECnLxvriPYPn7Zw3xHQePMviTWP6pBT8i\np1GQi3et/wYWjIMBMyCsSYlNii/4uVYLfkROoyAX79k8A359zrVvSgkfbMYlp7N4Syo/rd/DrVrw\nI3JGCnLxjpQVMH2462w+9dqe9nBccjr9Jy0nJ99JgDH0aFnHC0WK+AYdgCtlL207fHUP/HkCNOpY\nYpPliWnk5juLblksTzpUdvWJ+BgFuZStrFT49Da4+u9w0XVnbNa6fjUswGEgKNBBTFR42dUo4mM0\ntSJlJy8bvrgTLr4Nogecten8hAPcdGkD2jSoTkxUONGROlJF5ExsCXJjTE+gJdDEsqznjTF3AfWB\nK4FBlmVl2PG8Uo5ZFkwb6vpQ8+rnz9p0e2oWv2zYx9wne1CzSnAZFSjiu+yaWhkCJAB3GWNCAQNc\nC+wFjhZvaIwZYoyJNcbExsXF2VSOeN3iN+Dwbrj5bTjHBlf/nLWZId2jFOIibrIryIv/BRqgEvA5\n0AdoULyhZVkTLcvqZFlWp+joaMQPJcyEVZPhL59C0NnPoxm74xAbdx9hQNemZVObiB+wK8gnAB2A\n9UXXqwJhwDwg06bnlPIoNQF+fAzu+ASqNzhrU8uyGPtLPE/1akmloIAyKlDE99kyR25Z1hxgDjDe\nju8vPuJYhuvsPteNhiaXn7P5r5v2cSzfSd/LGpVBcSL+Q4cfij0sC6Y9As2vhg73nLN5fqGT12Yl\n8Fyf1jgcOkmEiCd0+KHYY/l/IHMv3D7FreZfrEyhcc3KdLtIKzhFPKUgl9K3cxUs/hcMnguB594z\nPDMnn7fnbuPjgeeefhGR02lqRUpX9iH4diDc/G+o2dStLhMXJdK9ZW3aNaxhb20ifkpBLqXn+Lx4\nm5uhzU1uddl/JIepy5MZ0auVzcWJ+C9NrUjpWTkRjh5wHWropvFztnDn5RE0DKtsY2Ei/k0jcikd\nB+Jh4WvQ778Q6N6KzC37M/ktfj8PX1XySZZFxD0KcrlwBbnw/WC45qUSTxBxJuNmbubhq1pQo3KQ\njcWJ+D8FuVy4+WOgRgR0vM/tLsu2H2TrgUzuiYmwsTCRikFz5HJhkhbDuq/g4aXn3AzrOKfTYtzM\nzYy8vjUhgVqKL3KhNCKX83csA34YCre+C1Vqu93t5w17McCNl5x97xURcY+CXM7fr3+Hltef9Uw/\np8otKOT1XzfztxvaaCm+SCnR1Iqcn62/QdIieGSZR92m/p5My7rV6NJcp24TKS0akYvnco7AT4/D\nLf+GkGpudzucnc/7C7bzzA2tbSxOpOJRkIvn5rwIzXu6djb0wH8WbuO6tvVoWc/98BeRc9PUingm\ncSFsnQ2P/O5Rt90Zx/hq1U5+/Wt3mwoTqbg0Ihf35WbB9GFw03io5NkGV2/OTuC+mEjqVT/7qd5E\nxHMKcnHfvFchoovrSBUPbNpzmMVbDzKkh5bii9hBUyvint2rYeN38OgKj7uOm7mZ4Ve3oGqIft1E\n7KARuZybsxB+fsJ17s3QWh51XbQllV3px7jzCi3FF7GLglzObdUkCK4K7e/yqFuh0+IfMzfzTO9W\nBAXoV03ELnqvK2d3ZI9re9oHZrm9l8px09bsJjQ4gOvb1bepOBEBjcjlXGY9C50GQp2WHnXLyS/k\nzdkJPNenNcbD/wBExDMKcjmzrXNg71ro9pTHXacs3cGljcOIjvRsTl1EPKcgl5LlZcOMp+DGNyHI\ns9OwHTqax8RF2xnZW+fhFCkLCnIp2ZJ/QaOO0OJaj7u+O28bN13akKg6VW0oTEROpSCX0x1Kch2p\n0muMx11T0rL5Yc0uhl9zkQ2FiUhJFORyutnPQ5fHoEYjj7u+PjuBgVc2o061EBsKE5GSKMjlZNvm\nwv6NriD30LqdGaxMSuPBbs1sKExEzkRBLv9TmA+z/gbX/wOCPNvcyrIsxv4SzxPXtiQ0WMsTRMqS\nglz+Z+VEqNEEWt3gcdd5mw9w6Gge/xfd2IbCRORsNHQSl6wDsPjN81rBWVDoZNzMzfzthtYEaim+\nSJnTX524zB3t2kvFwxWcAN/G7SK8ajBXt65rQ2Eici4akQvsinOdTPmxVR53zc4rYPxvW5h4byct\nxRfxEo3IKzrLcn3Aec0LUKm6x90nLU7iimbhtG8SZkNxIuIOBXlFF/8T5B+D9nd73DU1M5fJS5N4\nupeW4ot4k6ZWKrLCfPhtFPR5HRye/5/+9tyt9OvQmIjw0NKvTUTcphF5RRb3EYRFQItrPO66PTWL\nGRv2MuzqFqVfl4h4xJYRuTGmJ9ASaGJZ1vNF9w0BNgOxlmVl2/G84oHcTFj0OvT/5ry6vz4rgSHd\no6hZJbiUCxMRT9k1Ih8CJAB3GWNCjTH3Az0Az1eaiD2WvQNRV0GD9h53jUs+xPpdGQzo2rS0qxKR\n82BXkBcfphkgGlgAtAB6F29ojBlijIk1xsTGxcXZVI6cJHOfaxVnz7973NWyLMbMiOepXq2oFBRg\nQ3Ei4im7gnwC0AFYX3T9I6ANsAdYXLyhZVkTLcvqZFlWp+joaJvKkZMsGAeX9YeakR53/XXTPo7l\nO+nbwfOdEUXEHrbMkVuWNQeYA4wvdvdqO55LPJSaAPHT4bFYj7vmFzp5bVYCo29pR4BDi39Eygsd\ntVLR/DYarnwcQj0/l+aXK1NoXLMy3VvWsaEwETlfCvKKJPl32LcernjI465Lth1k3MzN3HpZQxsK\nE5ELoSCvKCwL5rwAVz/v8V7jccnpDJi8kqN5hTw/bSNxyek2FSki50NBXlHET4f8HLjkDo+7zt60\njwKnBUB+gZPliWmlXZ2IXAAt0a8ICvNdc+PnuRR/7c4MAh0Gy7IICnQQExVuQ5Eicr4U5BXBBSzF\nX7Qllb2Hc5j64BWsTskgJiqc6MiapV+jiJw3Bbm/u4Cl+HkFTkb9tIkXb2pLl+a16dK8tg0FisiF\nOuf7bGNMp1NuD7avHCl1F7AUf/LSJJqGV+HatvVKvSwRKT3uTJg+YozpaYy5wRjzO66VmuILLmAp\n/r7DOXywcDsv3tTWhsJEpDS5E+RJwDvADKAQeNzWiqT0XMBS/LG/xNO/cyRNa1exoTARKU3uBPlo\nYDtwRdH1y2ytSErH8aX43Z7yuOvv29OIS07n0Z7aa1zEF7jzYefdlmV9efyGMSbfxnqktJznUvz8\nQiejpm/i+RvbUDlYuxuK+IJzjsiLh3jR7QW2VSOl4wKW4k/9PZk61ULofXF9GwoTETvo8EN/cwFL\n8VMzc3l3/ja+fqgLxmh3QxFfoSX6/uYCluKPm7mZ26Mb06JuVRsKExG7KMj9yfGl+NeN9ngpflzy\nIZZuO8iway6yqTgRsYuC3J+c51L8QqfFiz9u4tk+rakaotk2EV+jv1p/kZsJC/8J93zrcdfPV6ZQ\nJSSQW9prr3ERX6Qg9xfL3oHmPT1ein/oaB5vzdnCp4M66wNOER+lIPcHx5fiD1nocdfXf03g5vYN\nadOgug2FiUhZ0By5PzjPpfjrd2XwW/x+nriupU2FiUhZ0Ijc1x1fiv9YrEfdnE6LF37cxMjrW1Gj\ncpBNxYlIWdCI3Ned51L8b+N24TBwW8fGNhUmImVFQe7LznMp/uHsfP75awKv3HoxDoc+4BTxdQpy\nX3UBS/H/NSeB69vV4+JGNWwqTkTKkoLcV53nUvw/9hxhxoa9jOjVyqbCRKSs6cNOX3R8KX6f1z1a\nim9ZFi9N38iT17WiZpVgGwsUkbKkEbkvOs+l+NPW7iYn38lfLm9iT10i4hUakfua81yKn5mTz7iZ\nm3n/nmgC9AGniF/RiNzXnOdS/H//tpUeLevQMaKmTYWJiLdoRO5Ljux1LcV/aJFH3bbsz+T7NbuZ\n/UR3mwoTEW/SiNyXLBgLHe51zY+7ybIsRk3fxPCrW1C7aoiNxYmIt2hE7isOxMPmX2BYnEfdZmzY\ny6GjedwT49k+LCLiOxTkvmLOi9DtKagc5naXo7kFjJkRz7/v7EBggN58ifgr/XX7gsSFcHALXD7I\no27vzd9G52a1uKKZZ/uwiIhvUZCXd04nzH4ernkJAt1fxJOYmsUXK1N4tk8bG4sTkfJAQV7ebfgG\nAoKh3Z/d7mJZFqN/+oNHrmpBveqe7cMiIr5Hc+TlWX4OzHsF+v0XPDgN25w/9rM74xgDrmxqX20i\nUm5oRF6erZjgWvgT2cXtLjn5hbz88x+MvqUdQfqAU6RCKPURuTGmJ9ASaGJZ1vNF99UHFlmWpXOK\nuSv7ECx7Gwb+6lG39xdsp33jMK5sUdumwkSkvLFjyDYESADuMsaEGmMaAZHARTY8l/9a+E9o1w9q\nu/9jS0nL5uPfd/DcjfqAU6QisSPIix9aYYCBwMNwYrR+EmPMEGNMrDEmNi7Os8UufittO6z/Cno8\n41G3V2b8weBuUTQKq2xTYSJSHtkR5BOADsB6YIJlWa9YljUAwLKs+ac2tixromVZnSzL6hQdHW1D\nOT5o9vNw5XCoWsftLvMTDrB1fyaDujWzsTARKY9KfY7csqw5wBxg/Cn3a+9Ud2yfDwf+gNs/crtL\nbkEho6dv4qVb2hESGGBfbSJSLumwhvKksABmPQu9XoVA9ze4mrQ4iRZ1q9GzVV0bixOR8kpBXp7E\nTYEqtaH1TW532Z1xjEmLE3np5rY2FiYi5ZkWBJUXx9Jh4Wtw7zSPFv+MnRHPfV2a0qRWqI3FiUh5\nphF5ebHgNddIvP7FbndZuu0g63Zl8PBVzW0sTETKO43Iy4PULbDha3h0pdtd8gqcvDR9Ey/c1JZK\nQfqAU6Qi04i8PPj1OfjTk675cTd9vGwHjcIq06ttPRsLExFfoBG5t22dA+lJcMXnbnc5cCSH/yzY\nxncPd8V4MJ8uIv5JI3Jvys+BmSPh+n+4vdd4XHI6A6as5KpWdYmqU9XmAkXEFyjIvWnZO1C3LbTs\n5VbzuOR07pq4nD/2ZjJzw17iktNtLlBEfIGC3FvSd8Dy/0Dvf7jdZcaGPeQVOgHIL3SyPDHNpuJE\nxJcoyL1l1rPQ5VEIi3Cr+YHMHH5at4egAEOAgaBABzFR4TYXKSK+QB92ekPCLEhNcHs/ley8Ah78\nKJb+nSPpdlEdliemERMVTnRkTXvrFBGfoCAva/nHXB9w3jTerf1UCgqdDPt8Da3qV+Pxay7CGKMA\nF5GTaGqlrC0ZDw0vgxbXnLNfQv78AAAJ3UlEQVTp8ZMo5xY4GfvnS3SooYiUSCPyspS2HVZOhKFL\n3Gr+38WJrNpxiK+HdiE4UP/nikjJFORlxbLgp8eh2wio0ficzX9ev4cpS3fw3cNdqV4pqAwKFBFf\npWFeWVnzKeRmQueh52wau+MQL/24iUn3d6KhTtsmIuegEXlZyDoAv42Ce3+AgLP/yBNTsxj66Wr+\n9ZfLaNewRtnUJyI+TSPysjDzGehwDzS49KzNDmbl8sBHqxjRqyU9Wrp/vk4RqdgU5HZLmAV71sBV\nfztrs2N5hQz6OJabL23InVe4t0hIRAQ0tWKv3EyY8RT0fQ+CzjzXXei0+OtXa2gaHspTvVqWYYEi\n4g8U5Hb6bTRE9YCoq87abMyMeA4fy+ftuzroWHER8ZiC3C6JCyDhF3h46VmbTV6SxKKtqXw3tCsh\ngTrTj4h4TnPkdsg5DD8+Bje/DZXPvJx+1sZ9fLBoO1MGXE6NUB0rLiLnRyNyO/z6nGsJ/kXXnrHJ\nmpR0nvthAx89cDlNaoWWYXEi4m8U5KUtYRYkLYKHl52xSXLaUYZMjeP1/7uUSxuHlWFxIuKPFOSl\nKfsQ/PxX6PdfCKlWYpP0o3k8MGUVw69uwTVtdOJkEblwmiMvLZblOtSwbV9o1q3EJjn5hQz+JJbr\n2tbj3i5Ny7Y+EfFbGpGXlrWfwYF46PufEh92Oi2e+mYd9WpU4pnercu4OBHxZwry0nBwK8x5Ee7/\n+YwLf16btZn9h3P4dFBnHA4dKy4ipUdBfqEKcuHbgXDVs1CvbYlNpi5PZs4f+/nu4a5UCtKx4iJS\nuhTkF2ruy1CjCVw+qOSH4/fz9tytfDu0CzWrBJdxcSJSESjIL8TWObDpB9cZf0pYWr9+VwZPf7ue\nD+/vRGR4FS8UKCIVgYL8fGWkwLRH4PYpEFrrtId3Hspm8Cex/KPfJXSI0MmSRcQ+OvzwfBTkwtf3\nQ9dh0PRPpz18ODufBz5axUPdm3N9u/peKFBEKhIF+fmY+QzUaOQK8lPkFhTy0KexdL+oDgP/1MwL\nxYlIRaOpFU+t/Rx2LIbB80+bF7csi2e+XU+NykH8/cY2XipQRCoaBbkn9q6H2c/DgBlQqfppD785\news70rL5YnAMATpWXETKiILcXZn74cu7oc8bUPf00faXK1P4af0evnu4K5WDday4iJQdW+bIjTE9\njTEPGWNeLbp9mTFmnjFmuTHG9w6mzs+Br/q7TqB8cb/THl64JZU3Zm9hyoDLqV01xAsFikhFZteH\nnUOABOAuY0wo8AdwDVDz1Oc0xgwxxsQaY2Lj4uJsKucCWBZMf8y16KfHM6c9vGnPYZ78ai0T7ulI\nVJ2qXihQRCo6u4K8+KjbWJaVB1wN9Aes4g0ty5poWVYny7I6RUdH21TOBVj8BqRtd22GdcqHm79u\n2sedHyxnQNemdGp6+rHkIiJlwa4gnwB0ANYDE4wx1wATgSlAK5ues/Rt+BZiP4I7Pz9tM6xvYncy\ndGocWbkFvLdgG3HJ6d6pUUQqPFs+7LQsaw4wBxhf7O7mdjyXbbbPg1l/g/umQ/UGJ+7OK3Dy7ryt\nTFqcBLjeXuQXOFmemEZ0pFZwikjZ04KgkuxeDd8Nhjs+OWlHw017DnPLu0vYtOcIb9/VgZAgBwEG\nggIdxESFe7FgEanIdPjhqQ5ugy/uhFvehsiuAOQXOnlv/jam/p7Mc33a0K9jI4wxfDYohuWJacRE\nhWs0LiJeoyAv7vBu+LQf9Pw7tL4RgPi9R3jq63XUqx7CjOHdqF+j0onm0ZE1FeAi4nUK8uOO7IWP\nb3LtKx59P/mFTiYs2M5Hy3bwzA2tuT26MaaErWpFRLxNQQ6QuQ8+vhk63gdXDidhXyYjvllHzSrB\n/DTsTzQMK/n0bSIi5YGCPOuAK8Tb/4WCLo/zwfxtfLgkiZHXt+IvlzfRKFxEyr2KHeSHd8HUP8PF\nt7G11VBGvL+MapWC+GnYn2ikUbiI+IiKG+Rp2+GTvjivGMzEghuZOHE5T/Vqyd1XRGgULiI+pWIG\n+b4N8On/sb/TCIaua0vloFR+fPRKmtQK9XZlIiIeq3hBnrQY69sHmNvsaUYubsoT1zaif+dIHNo/\nXER8VMUK8rWfU/jr87wS8hTxaZcy7ZH2RIRrFC4ivq1iBLllYc17lczYLxiQ8zy3druaF2M0ChcR\n/+D/QZ6bxdFvhrI7eSv/DHuT8YO6ExlexdtViYiUGr8OcueBLRz++C8sPBrJ4Z4fMbFba43CRcTv\n+G2Qp678luBZT/Jllfu5/pGRRNWt5u2SRERs4XdB7szLIf7zp6m14xeWXv4eQ264SWe0FxG/5hdB\nHpeczvLENC4N3kOT+Y+TGVCf0IHz6RsR4e3SRERs5/NBHpecTv9Jv3OncxZtA79nUcSj3DLgGQIC\ndM4MEakYfD7tliemEeNcQ9+AJdyeP4o9UbcrxEWkQvH5EXlMVDjvODpwR/4lOAKDdMo1EalwfD7I\noyNr8tmgLjrlmohUWD4f5KBTrolIxabJZBERH6cgFxHxcQpyEREfpyAXEfFxCnIRER+nIBcR8XHG\nsixv13CCMWYSsOs8u0cDcaVYTnmm1+qf9Fr90/m+1saWZQ1yp2G5CvILYYyJtSyrk7frKAt6rf5J\nr9U/lcVr1dSKiIiPU5CLiPg4fwryid4uoAzptfonvVb/ZPtr9Zs5chGRisqfRuQiIhWSglxExMcp\nyEVEfJyCXETEx/lFkBtjehpjHjLGvOrtWuxW9FpXGmN+8HYtdjPGtDHGrPB2HWXBGBNojBlmjLnC\nGBPs7XrsZIwZYYy53xjzozHGeLseOxhj6htjphtjehhjXjfGPG2MaWPX8/lFkANDgATgLmNMqLeL\nsdlSoDfQzNuF2MkYcxHQHKjn7VrKyN+BlsCNQIGXa7FbINATSLT89LA5y7L2AdWB1sDlQCpwn13P\n5y9BXnwE45f/wx9nWVYecA3QzxgT4u16bPQYMACobYy5xMu1lIVoYDrQB7jMy7XYLQiYBvyfMaay\nt4uxWZlkk1+csxOYAHQAvrcs66i3i7GTMeYe4GUgF+hSdOl3LMt63BjTFOhkWdYGL5dTFt7H9e8Z\nC/zh5VrsVh0IAX4D8r1ciy2MMWFATeAIrn/PxsBntj2fn76zERGpMPxlakVEpMJSkIuI+DgFuYiI\nj1OQi4j4OAW5yCmMMfWMMRXl+HXxAwpykWKMMbWA74EG3q5FxF0KcpGTPQJ0BR41xjT3djEi7lCQ\ni5xsSdHle5ZlbfdqJSJuUpCLiPg4BbnIyZzHrxhjArxZiIi7FOQiJ1sPbAXuB/x9QyfxE9prRUTE\nx2lELiLi4xTkIiI+TkEuIuLjFOQiIj5OQS4i4uMU5CIiPk5BLiLi4/4f8v4TKmhvhu0AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f41d70e3128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from numpy import *\n",
    "from matplotlib.pyplot import *\n",
    "\n",
    "# time intervals\n",
    "dt = 1.0\n",
    "tt = arange(0, 10, dt)\n",
    "# initial condition\n",
    "xx = [0.1]\n",
    "\n",
    "def f(x):\n",
    "    return x * (1.-x)\n",
    "\n",
    "# loop over time\n",
    "for t in tt[1:]:\n",
    "    xx.append(xx[-1] + dt * f(xx[-1]))\n",
    "\n",
    "# plotting\n",
    "plot(tt, xx, '.-')\n",
    "ta = arange(0, 10, 0.01)\n",
    "plot(ta, 0.1 * exp(ta)/(1+0.1*(exp(ta)-1.)))\n",
    "xlabel('t')\n",
    "ylabel('x')\n",
    "legend(['approximation', 'analytical solution'], loc='best',)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why use scientific libraries?\n",
    "\n",
    "The method we just used above is called the *Euler method*, and is the simplest one available. The problem is that, although it works reasonably well for the differential equation above, in many cases it doesn't perform very well. There are many ways to improve it: in fact, there are many books entirely dedicated to this. Although many math or physics students do learn how to implement more sophisticated methods, the topic is really deep. Luckily, we can rely on the expertise of lots of people to come up with good algorithms that work well in most situations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Then, how... ?\n",
    "\n",
    "We are going to demonstrate how to use scientific libraries to integrate differential equations. Although the specific commands depend on the software, the general procedure is usually the same:\n",
    "\n",
    "* define the derivative function (the right hand side of the differential equation)\n",
    "* choose a time step or a sequence of times where you want the solution\n",
    "* provide the parameters and the initial condition\n",
    "* pass the function, time sequence, parameters and initial conditions to a computer routine that runs the integration.\n",
    "\n",
    "### A single equation\n",
    "\n",
    "So, let's start with the same equation as above, the logistic equation, now with any parameters for growth rate and carrying capacity:\n",
    "\n",
    "$$ \\frac{dx}{dt} = f(x) = r x \\left(1 - \\frac{x}{K} \\right) $$\n",
    "\n",
    "with $r=2$, $K=10$ and $x(0) = 0.1$. We show how to integrate it using python below, introducing key language syntax as necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f41cd808b00>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEFCAYAAAD69rxNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHQxJREFUeJzt3Xt4VOWh7/HvOyTchWC4BZIAgXJV\nuQQ1YC2l1ZSinlKsHupdxNDtrbs+itWebuMFDt1bj+yqtbBRtlUK1aoUEVTk5gUDJBCEGoISSABB\nQ0h2CJDrrPPHQIhcJAkz886a9fs8zzwzk5nM+1uT8GPlnXUxjuMgIiLu47MdQEREmkcFLiLiUipw\nERGXUoGLiLiUClxExKVU4CIiLqUCFxFxKRW4iIhLqcBFRFwqxnaA46ZMmeIkJibajiEiYt1jjz32\nouM4U872vIgp8MTERDIzM23HEBGx7rHHHtvTmOdpCkVExKVU4CIiLhUxUygnKysrY9++fbZjeFJC\nQgJxcXG2Y4jIWURsgR84cIDevXvTpk0b21E85ejRo+zdu1cFLuICETuFUlNTQ+vWrW3H8JzWrVtT\nU1NjO4aINELEFjiAMcZ2BM/Rey7iHhFd4I2RU1jK86u+JKew1HaU73TZZZeRk5PT5O/Lz88HYNq0\nafzxj38MdiyxzNbvr41xvbSs4RKSOXBjTHdgDvA0cDXwDbDEcZy8YI6TU1jKjXOzqK710zLGx/wp\naaT26hTMIYJm3rx5pKSkNOl7tm/fzr333sv777/Pb37zG1q1ahWidAKB36esghLSUuJJTY6D6sNQ\ndQiqK6C2Emqrjl2O3a5reL86cN9fC/66wMU5du2vPXbb3+B2LQfKj1KYv58eTh1frXLo0y+e89vG\nguMAzhmuOXG/4e0mXJdX1lKzr5zhjkPNakN5Qgc6tA7tx2GHKmuo3VfOCAdqV8OhhA6c1zq2wTIE\nw6mvc6iyBv/+Q4x0HPyrDYe6n3di3DO+TJDyjJkGfccG57XOICQ/Ncdx9htjOgADgYuB/wZuAR5u\n+DxjTAaQAXD11Vc3eZysghKqa/34Haip9ZNVUNKsAl+3bh0ffPABH330EaNHj+bRRx9l5syZLFiw\ngBUrVvDggw+Sm5vLiBEjaN++PTNnziQ1NZW0tDQ+/fRTli1bxiOPPMKQIUPo2rUr48aNY9SoUWRk\nZFBRUcHdd9/N6NGjmTt3Llu2bGH27NlMnjyZlStXcscddzBjxgw2bdrEggUL8Pv9zJ07lzfffJO5\nc+eSl5fHO++8w6xZsxgyZAgPPPAAjzzyCD169OAHP/gBF154IcnJyUyfPp0XX3yRVatWkZyc3OT3\nIKo5DhwpgbIiOHwADhc3uATuV5QfpMM3xfycI5y3+iiOrwoT0wZanQct20FMa4hpdey65Yn7LVqd\n+HqLluBrcewSA6ZF4Os+X+C2L6bBYz62HCohq64jtY4PjI/WLRO4sn8CHJ/GMgYwZ7nmLM859fGV\nuV/xetEe/EAL4BddE5kwrGfz3ttGTrmt2LSX14p24xD4s/+6rklMGN5wzCBN3Z2U54NNe/nb7t04\nDvgMXN89iZ8Pb8yyBiFP5/7n/hpnEeqtUFo2uH3KO+I4zhwCa+pkZmY2+b+9tJR4Wsb4qKn1Exvj\nIy0lvlkhKyoqSExMpKamhk2bNgEwdepUli9fzmuvvUZycjIHDx4kMzOzviy7dOlCx44d2bBhA6+8\n8gplZWVMmjSJ4cOHM3nyZB5//HHuv/9+du7cSVxcHO3btwcgKSmJzp07M3HiRJ544gk+/vhjJk+e\nTEFBAdXV1fTu3ZuioiI2btzI4MGDiY+P56qrruL1118H4LnnnuO8884jPT2d++67j61btwIwfvx4\n5syZQ25urncL/PAB+CYPircFLmVFJy4xraBjIrTvBu26QLvO0L4rdBsCbTvz7ucVzN37DeX+Nhwx\nbci48iLu+tHAkMbt0LWUxXlZgd/fFj5uHJUGYfgLMqlFKTmbsur/3fwmNfTjJvlK2dhgzPvDMCZA\nsikld+OJ9/iBEeEZN1xCNYUSB3QCyoHPgURgfrDHSe3ViflT0k782dvMH4zP52PNmjX07NmT8vLy\n+q/HxsZ+60O94/ePfy05OZl27drh9/sxxlBXV0d1dTWO45Cenk5sbCyLFi3itttuO2sGv99PZWUl\nO3bsoE2bNjiOc9oPFE8e63SPe0JlOezNht3rA5f9nwWmMboOgq4DoctA6HcFxCVDxyRo3eE7X65P\nq1J25WRR4wQK5tK+3UK+CMH6/XXDuF5a1nAK1RRKGTD02N2XQzHGcam9Op3zD2XDhg1s27aNQYMG\n1ZfmK6+8Qnl5Oddeey3PP/88u3fv5umnn+a+++7D5/NRXFzM7t27AbjhhhtYs2YNixYt4g9/+APF\nxcXMmjWLjIwMHn30UdLT06moqKCoqIiKigrKysooKCgAYP369QDs27ePlStX0rVrV3r27ElBQQFj\nx47l6NGj5OTksHfvXg4dOsRTTz3F73//e1avXs2MGTPYsydwyITt27dTUVFRfz/q+Otgbw58+QF8\nsRyK8yFhKCRdDCMnQ4/h0KFHo/+kP5nNgrFRKjbG9dKyhotxgvYBwrnJzMx0Gh7MKi8vj0GDBoU9\nx65du+jTpw+lpaX1O7NkZmaSm5vLokWLwp7HBlvv/XH1Hyb26USq7wv47DX4fBG07w7fuwL6XQlJ\nlwbmokWikDHmMcdxMs/2vIjdE9OWwsJCAAoKChgxYkT91/bv38+RI0do27atzXhRL6ewlLvnvs9E\nZyXdWqzgaKc42oyYBFNWwPl9bMcTiSgq8JOMGTOGk/8qmTdvnqU0HnOwgDbLMnnPt5z360Zyd82v\nSb9oHHf/4Hu2k4lEJBW42Fe2Gz78d8hbQtyAmxm3dxbf1LUnNsbHv/XtbDudSMRSgYs9tVWw9o/w\n6Z9g5O1wbw492p7Pc8NLo3arAZFgUoGLHbs+gcX3Bjb3y1gNnXrVPxTNWw2IBJPrj4Viw5o1azDG\nUFZWdtrHjx+/pLHHP7n99tuZMGFCo8Y+cuQIQ4YMYdasWWd8Tn5+Pn6/n5SUFPbv39+o1w2buhpY\n8Tj8/XZIfxJ++ddvlbeINJ4KvBl69Tpz4bz55pvMnj0bCHz4OXTo0DM+tzGvd7K2bdvSpUuXMz7+\nzDPPsGzZMnw+H0uWLKF79+6Nfu2QK9sNL6bD/i3wq49h4HjbiURczfMFvm7dOqZPn864ceNYunQp\nu3fvxhjDjBkz6Nu3L0VFRac857ja2lqSkpKYPn06H374IS+//DLPPPMMmzdvZuPGjYwePZolS5aQ\nl5fHgw8+yJgxY9i5cyePP/448+bNY9SoURw6dOiUTEVFRdx5551MmDCB2tpa7r33XubMmcNdd931\nrS1knnzySYYNG8b8+fOJi4ujrKyMP//5z6xdu5Y1a9YwZMgQcnNzefXVV8nMzOTmm2+mtLSUJ554\ngqSkJKZNm8akSZPC8j5TtA7mXgFDJsANrwV2YxeRc+KOOfDMjkF6nf855UsNj4OyfPlyxo8PrBU2\nPLZIu3btvvWcwYMHAxATE8M999zDihUruPTSS7n55ptZtWoVcXFx9Qe+gkDRXn755aSnp9O2bVuM\nMSQnJ5OVlVW/R2ZDZWVlLFy4kJkzZ5Kfn89LL73E4cOHSU5OZvLkyfXPS0xMBKBnz8DBeeLi4khI\nSGD06NGMGTOm/nkPPPAAb7/9Ns8++ywvvvgiSUlJxMfHM378eCZOnBic9/a7bPk7LHsIJrwA/dND\nP56IR7ikwE8t3mBpeByU0+2V6vf7v/M5kydPJjMzk5tuugmfz3fa45dUVlaSnZ3Nr371K/x+PyUl\nJfXHXDndmD169GD27NlMnTqVTz75pP41Tz7+ic936h9QTT1+SsiPnZLzMqz+v3Dr4sCBo0QkaDw/\nhXL8OCixsbEUFhbWH9+k4bFFTn7O8b01i4qK6NKlCxMmTKif605LSyMvL4/CwsL6459MmzaN999/\nn4cffpiSkhJWrFjB2rVr6d+/Pzt27KCwsJCvv/6aqqoqALZu3crmzZu5/vrr6dWrFzfeeCMvvfQS\nV1xxBUOHDqW4uJiioiJSU1M5cuQIy5cvp6qqigMHDnDZZZexfv36b2WcOXNm/dTPHXfcwZ49eygr\nK2Pnzp313xcSWS/Ah0/Bbe+ovEVCQMdCOUf79u1jy5YtpKdHz9RAUN77ja/Amn+H25dCXFJwgol4\nRGOPheL5NfBzUVFRwTXXXEPHjkGao48WeW/Dyifh5rdU3iIh5I458AjVvn17srOzbceILLvXw9v/\nCje9AZ372U4jEtUieg08UqZ3vOSc3vPyr+C1W2DCn6DHsOCFEpHTitgCj42NpbKy0nYMz6msrCQ2\n9iwnfT2dmkpYeCNccif0/0nwg4nIKSJ2CqVz587s2rXLdgxPSkhIaPo3vfvbwOnLvn9/8AOJyGlF\nbIHHxcXVnxFHIlzeEtixEn71UbNPaSYiTRexBS4uUf4VLPlXmPRXaK2tcUTCKWLnwMUFHAcW3QWX\nZEDSJbbTiHiOClyaL/evcPSg5r1FLNEUijRPRTF88Cjc+HdooV8jERu0Bi7N897DMHSStvcWsUir\nTtJ0BWtg9zq4K8t2EhFP0xq4NI2/Dt57BK58Alq2s51GxNNU4NI0m16FVh1g8M9sJxHxPE2hSONV\nHYJVM+CGhdphRyQCaA1cGu+T/4S+Y6HHcNtJRAStgUtjHS6BDXNh6oe2k4jIMVoDl8ZZ+58wZGLg\ngFUiEhG0Bi5nV/FN4OTE/7LWdhIRaUBr4HJ2H8+Ci/43dOxpO4mINKACl++0Of9LKrP/wmd9JtuO\nIiInUYHLGeUUlvLR/BksqhrJ9X/dSU5hqe1IItJAyAvcGPOAMeZWY8w/jNHGw26S/cUeJpnl/Fft\neGpq/WQVlNiOJCINhGMNPAYYCxQ4J50x1xiTYYzJNsZk5+TkhCGKNMVPa1eRS3920ZPYGB9pKfG2\nI4lIA+Eo8FhgEfALY0ybhg84jjPHcZyRjuOMTE1NDUMUaTR/Hcn5L5J41W+5P30A86ekkdqrk+1U\nItJAODYj7AC0Aj4AasIwngRD3mJo15WBl1zJQNtZROS0Ql7gjuM8eOzm30I9lgTRujkw6m7bKUTk\nO2grFDnV1/+E0p0w8CrbSUTkO6jA5VQbXoTU26BFrO0kIvIdVODybZXlsPUNGHGr7SQichYqcPm2\nz/4GKWOgQ4LtJCJyFipwOcFxAtMnF0+xnUREGkEFLicUZYG/FnpfbjuJiDSCClxOyH0Vht+k06WJ\nuIQKXAKqD0Pe2zB0ku0kItJIKnAJ+HwxJKXBed1tJxGRRlKBS0DufBh+o+0UItIEKnCBgzvhm8+h\n/09tJxGRJlCBC2xeABf8AmJa2k4iIk2gAvc6xwkU+LAbbCcRkSZSgXvdng0Q0xoShtpOIiJNpAL3\nui2vw4XXadtvERdSgXtZXS38cxFccK3tJCLSDCpwL9v1IXTsCfF9bScRkWZQgXvZljcCW5+IiCup\nwL2qphK2LYELJtpOIiLNpAL3qi+XQ/cLoUMP20lEpJlU4F615e/68FLE5VTgXlR9BHashEH/y3YS\nETkHKnAv2rECegyDdvG2k4jIOVCBe9Hni7X2LRIFVOBeU1sFX7wHg66xnUREzpEK3GsK1kCXQTpx\ng0gUUIF7Td5iGKzpE5FooAL3krpayF8KA6+2nUREgkAF7iWFn0DHJOjUy3YSEQkCFbiXaPpEJKqo\nwL3C74e8JTDoZ7aTiEiQqMC9Ys8GaNMJOveznUREgkQF7hX578AgfXgpEk1U4F6R/y4M+KntFCIS\nRCEvcGNMjDHmXmPMJcaYlqEeT06jZAdUlkHCcNtJRCSIYsIwxu+Azscu2WEYT062/V3oPw58+oNL\nJJqE4190KrAYGA8Ma/iAMSbDGJNtjMnOyckJQxSPyl+m6RORKBSOAn8BGEVg7fvzhg84jjPHcZyR\njuOMTE1NDUMUDzpaCl/lQp8xtpOISJCFfArFcZxlwLJQjyNn8MUH0Pv70LKt7SQiEmSaFI12+Us1\nfSISpVTg0ay2OnD2nf4/sZ1EREJABR7NitZCfD8d+1skSqnAo5l23hGJairwaOU4gfnv/ipwkWil\nAo9WxdsCJd5tiO0kIhIiKvBolb8UBowDY2wnEZEQUYFHK81/i0Q9FXg0qvgGivOh1/dtJxGREFKB\nR6Pt70HfsRCjgz+KRDMVeDTarukTES9QgUebmkrY+SH0u9J2EhEJMRV4tCn8GLoOhnbxtpOISIip\nwKPN9vd07BMRj1CBRxPHOXH2HRGJemctcGPMyJPu3xm6OHJOireBA3QdZDuJiIRBY9bA7zLGjDXG\n/NQY8ynw51CHkmba/m5g+kR7X4p4QmMKfCfwLPAOUAf8OqSJpPm2v6fpExEPaUyBPwbsAC45dnvY\ndz9drDhyEL7+Z+D0aSLiCY05J+YNjuMsPH7HGFMTwjzSXF8shz4/gNjWtpOISJicdQ28YXkfu786\nZGmk+Y7Pf4uIZ2gzwmhQVxM49+X30m0nEZEwUoFHg6IsOD9F574U8RgVeDTQzjsinqQCjwbafV7E\nk1TgbleyA6oOQfehtpOISJipwN1u+3vQPx18+lGKeI3+1bud5r9FPEsF7maV5bB3I6T80HYSEbFA\nBe5mO1ZAchq0bGc7iYhYoAJ3M219IuJpKnC38tfBF++rwEU8TAXuVns2QPvuEJdsO4mIWKICd6tt\n78DAq2ynEBGLVOBu5DjHCny87SQiYlFYCtwYM8gYsy4cY3nCge1QWwkJOreGiJeFvMCNMd8D+gLd\nQj2WZ2xbAgPG69yXIh4XjjXwe4DbgM7GmAsbPmCMyTDGZBtjsnNycsIQJUpsW6rpExEJfYE7jvNr\n4AHggOM4W056bI7jOCMdxxmZmpoa6ijRoXwflHwBvXTuSxGva8w5Mc+Z4zi7gN7hGCvqbV8G/a6E\nmJa2k4iIZdoKxW22LdXmgyICqMDdpepQ4PRp/a6wnUREIoAK3E2+/ACSL4XWHWwnEZEIEJY5cDl3\nOYWlnLdmIa36/pBetsOISETQGrgL5BSWctvcj+n69Ufc/PH55BSW2o4kIhFABe4CWQUlXOzfzJdO\nT/bWxpFVUGI7kohEAE2huEBaSjwJLdazzH8psTE+0lLibUcSkQigAneB1J7tGNpmM+XDH2L+oMGk\n9upkO5KIRAAVuBvsXENMl/7cNu4y20lEJIJoDtwN/rkIhkywnUJEIowKPNLV1UD+OzD4Z7aTiEiE\nUYFHuoI1EN8POibaTiIiEUYFHuk+fwsGa/pERE6lAo9kdTWBg1dp+kRETkMFHskKVkN8X4hLsp1E\nRCKQCjySffYaXHi97RQiEqFU4JGqqgK2vwcXTLSdREQilAo8Um17B5LToF1n20lEJEKpwCPVZ3+D\nizR9IiJnpgKPRBXfwN5sGKAzz4vImanAI9HWNwLl3bKt7SQiEsFU4JFI0yci0ggq8EhTnA/l+6DP\nGNtJRCTCqcAjzca/wLAbwNfCdhIRiXAq8EhSWw2bF8Lwm2wnEREXUIFHkvyl0HVQYPd5EZGzUIFH\nko1/gRG32k4hIi6hAo8UpYXw1SYYdI3tJCLiEirwSJE7Hy68DmJb204iIi6hAo8EdTWw8RUYcYvt\nJCLiIirwSJD3NpzfB7pfYDuJiLiICjwSrJsNl061nUJEXEYFbttXufA/e2DAVbaTiIjLqMBtWz8H\nLr4DWsTYTiIiLqMCt6miGLYt0bbfItIsIS9wY8xYY8x6Y8xboR7LddbPhiEToV287SQi4kLhWAP/\nBBgH9Dn5AWNMhjEm2xiTnZOTE4YoEaTqEGS/BJfdZzuJiLhUyAvccZxq4MfARGNMq5Mem+M4zkjH\ncUampqaGOkpkyZ4HKT+E81NsJxERlwr5J2fGmJuAx4EqYNSxa2+rrYKsP8ENr9lOIiIuFvICdxzn\nVeDVUI/jKpsXQLcLIOEi20lExMW0FUq41VbBh0/BmGm2k4iIy6nAwy3nv6HrYEi6xHYSEXE57T0S\nTtWH4aP/Bze+bjuJiEQBrYGH0/o5kJymuW8RCQqtgYfLkYOw9jm4fantJCISJbQGHi6rpsMFE6HL\nANtJRCRKaA08HL7+J3z+D7h7ve0kIhJFtAYeao4Dyx6CMQ9B2/NtpxGRKKICD7WtbwTmv1Nvt51E\nRKKMplBC6XAJvPsw/HKBjvctIkGnNfBQeve3gTPNJ460nUREopBWC0Nl+3uwZz38y1rbSUQkSqnA\nQ+HQ17D4Prh2LrRsZzuNiEQpTaEEm98Pb2XAiFugz+W204hIFFOBB9snswJHHBzzkO0kIhLlNIUS\nTDtWQtYLkLFKW52ISMipZYLlwBfwZgZc9zJ0TLSdRkQ8QFMowXDkICyYBD/+N+h9me00IuIRKvBz\nVXUI5l8HA8YHPrgUEQkTFfi5qDkKC34J3YbAlY/bTiMiHqMCb67qw7DwRmjfDa5+BoyxnUhEPEYF\n3gy523ey77lxHDBx8PPZ4GthO5KIeJAKvIm2bP2MtvOvYcnBJC7Pv5acPYdsRxIRj1KBN0XBavr+\nYwIL637I9NobqK6FrIIS26lExKO0HXhj1NXCJ8/Aujns+dGz/HWpjxbGT2yMj7SUeNvpRMSjVOBn\nU7ID3poKsW0gYxX9OyYyv0cpWQUlpKXEk9qrk+2EIuJRKvAzqamET5+DT5+HMdPgkqngC8w4pfbq\npOIWEetU4CdzHMh7G5b/HrpdAHeuhPP72E4lInIKFfhxfj9sexvW/Edgm+6rnoZ+V9hOJSJyRirw\no2WweQFkvwSxbeFHv4P+47RjjohEPG8WeG01FKyCrW/C9mWBNe2rZ0Gv0SpuEXEN7xR4RTEUrIYd\nK2D7u9B5AAz5OaQ/Ae272k4nItJkri7wnMIzbM7n98PBHbA3J3Ap/BTKiqD396HvWPjR76FjT3vB\nRUSCwLUFnlNYyk1z19K5tphNMft4dFQsSXV7AidW+HoLtO4IPVMDlwt+EbjWWXJEJIq4ttGyCkro\nW7eL/2r5FAVOAqV7LiBp2MWBaZGug6F9F9sRRURCKuQFbowZC/QHkhzH+T/Bet20lHiebZHC96uf\nIzbGx/yfpIF2rhERDwnHwawygHzgl8aYtg0fMMZkGGOyjTHZOTk5TXrR1F6dmD8ljfvTBzB/Spr2\njBQRzwlHgbdscPtb2+g5jjPHcZyRjuOMTE1NbfILp/bqxN1j+6m8RcSTwlHgfwaGA286jnM4DOOJ\niHhCyOfAHcdZDiwP9TgiIl6jEzqIiLiUClxExKVU4CIiLqUCFxFxKeM4ju0MABhj5gJ7mvGtqUDT\nNiJ3Ly1rdNKyRqdzWdZEx3GmnO1JEVPgzWWMyXYcZ6TtHOGgZY1OWtboFI5l1RSKiIhLqcBFRFwq\nGgp8ju0AYaRljU5a1ugU8mV1/Ry4iIhXRcMauIiIJ6nARURcSgUuIuJSKnAREZdydYEbY8YaY6Ya\nY560nSXUji3remPMW7azhIMxZpAxZp3tHKFmjIkxxtxrjLnEGNPy7N/hXsaYB4wxtxpj/mGMMWf/\nDvcxxnQ3xiw2xowxxvyHMeZBY8ygUI3n6gLnO07XFoU+AcYBfWwHCTVjzPeAvkA321nC4HcEzhl7\nFVBrOUuoxQBjgQInSjd/cxxnP9ABGAhcDBQDt4RqPLcX+BlP1xZtHMepBn4MTDTGtLKdJ8TuAW4D\nOhtjLrScJdRSgcXAeGCY5SyhFgssAn5hjGljO0yIhaWbQn5GnhDzzOnajDE3AY8DVcCoY9dRyXGc\nXxtjegMjHcfZYjlOqL1A4OeZDXxuOUuodQBaAR8ANZazhIQxJg7oBJQT+HkmAvNDNl6U/iUjIhL1\n3D6FIiLiWSpwERGXUoGLiLiUClxExKVU4CLHGGO6GWO8sO25RAkVuAhgjDkfeBNIsJ1FpLFU4CIB\ndwGjgbuNMX1thxFpDBW4SMDHx66fdxxnh9UkIo2kAhcRcSkVuEiA//gNY0wLm0FEGksFLhLwGfAF\ncCsQ7QdakiihY6GIiLiU1sBFRFxKBS4i4lIqcBERl1KBi4i4lApcRMSlVOAiIi6lAhcRcan/D/1Y\n6pSLt2bIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f41cd96a5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# everything after a '#' is a comment\n",
    "\n",
    "## we begin importing libraries we are going to use\n",
    "# import all (*) functions from numpy library, eg array, arange etc.\n",
    "from numpy import *\n",
    "# import all (*) interactive plotting functions, eg plot, xlabel etc.\n",
    "from matplotlib.pyplot import *\n",
    "# import the numerical integrator we will use, odeint()\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# time steps: an array of values starting from 0 going up to (but\n",
    "# excluding) 10, in steps of 0.01\n",
    "t = arange(0, 10., 1.)\n",
    "# parameters\n",
    "r = 2.\n",
    "K = 10.\n",
    "# initial condition\n",
    "x0 = 0.1\n",
    "\n",
    "# let's define the right-hand side of the differential equation\n",
    "# It must be a function of the dependent variable (x) and of the \n",
    "# time (t), even if time does not appear explicitly\n",
    "# this is how you define a function:\n",
    "def f(x, t, r, K):\n",
    "    # in python, there are no curling braces '{}' to start or \n",
    "    # end a function, nor any special keyword: the block is defined\n",
    "    # by leading spaces (usually 4)\n",
    "    # arithmetic is done the same as in other languages: + - * /\n",
    "    return r*x*(1-x/K)\n",
    "\n",
    "# call the function that performs the integration\n",
    "# the order of the arguments is as below: the derivative function,\n",
    "# the initial condition, the points where we want the solution, and\n",
    "# a list of parameters\n",
    "x = odeint(f, x0, t, (r, K))\n",
    "\n",
    "# plot the solution\n",
    "plot(t, x, '.')\n",
    "xlabel('t') # define label of x-axis\n",
    "ylabel('x') # and of y-axis\n",
    "tt = arange(0, 10, 0.01)\n",
    "# plot analytical solution\n",
    "# notice that `t` is an array: when you do any arithmetical operation\n",
    "# with an array, it is the same as doing it for each element\n",
    "plot(tt, K * x0 * exp(r*tt)/(K+x0*(exp(r*tt)-1.)))\n",
    "legend(['approximation', 'analytical solution'], loc='best') # draw legend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a much better approximation now, the two curves superimpose each other!\n",
    "\n",
    "Now, what if we wanted to integrate a system of differential equations? Let's take the Lotka-Volterra equations:\n",
    "\n",
    "$$ \\begin{aligned}\n",
    "\\frac{dV}{dt} &= r V - c V P\\\\\n",
    "\\frac{dP}{dt} &= ec V P - dP\n",
    "\\end{aligned}$$\n",
    "\n",
    "In this case, the variable is no longer a number, but an array `[V, P]`. We do the same as before, but now `x` is going to be an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of x: (50, 2)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f41cd652ba8>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEFCAYAAAD69rxNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXecZGWV97+3Qld1dZ7pCUyeASaQ\nhxmUoEgSAwYMKOiKKK6LadVd2SAbmN3XsO6q77uoKLqrYgJRRBFUkAwDyszADAyTc+4cKqfn/eOp\n21PdXeHe6rrPvT39fD+f+cxMdXXV033rnuc85/zOOYYQAo1Go9FMPnxuL0Cj0Wg0taENuEaj0UxS\ntAHXaDSaSYo24BqNRjNJ0QZco9FoJinagGs0Gs0kRRtwjUajmaRoA67RaDSTFG3ANRqNZpIScHsB\nJh/5yEfEvHnz3F6GRqPRuM6aNWv+RwjxkWrP84wBnzdvHrfeeqvby9BoNBrXWbNmzUErz9MhFI1G\no5mkaAOu0Wg0kxTPhFA0Go1GFQMDAxw5csTtZQAQDoeZN28ewWDQ9vdqA67RaKYcPT09LFq0iMbG\nRlfXIYSgt7eXgwcPsnjxYtvfr0MoGo1mypHJZAiHw24vA8MwmD59Oslksqbv1wZco9FMSQzDcHsJ\nwMTWoQ24RqPROMSaNWswDIPnnnsOgJ/+9KesWbOmbq+vDbhH+ef7XiaX1+PuNJrJzM0330x7ezsb\nNmwAIBQKccstt9Tt9XUS04Pk8oIfPbePv3/TcppD+hJpNE6y6B8eqMvr7P3yVeMei0QifOADH+A7\n3/kOH/vYxzAMg0Cgfve0tg4eJJnJAZDK5LQB9xjP7+3jzLlthIN+t5eiqROlDG89uemmm7jtttv4\nyle+wo033ljX19YhFA+SMA14Nu/ySjRjWXP/Zl4+NOj2MjSTiNNOO42LL76Y7du309nZWdfX1u6d\nB0mktQH3KvF0Tl8XjW0+9rGPsXTp0rq/rjbgHuS4B55zeSWasSTTOX1dNLa59tprHXldHULxICMe\neEZ7el4jkcnp6+JBHt16jFcOD7m9DOVoA+5B4jqE4lkSGR1C8SIPbDrKn/f0ur0M5WgD7kGSOoTi\nSfJ5QTKTH7k+Gu+QzORITsGNVRtwD2LGwJP6qO4pkll9MvIqUzW0pQ24BzkeQtGenpdI6OviWeLp\n7MgGO5XQBtyDjKhQpqBH4WXiOrnsWRKZvCevy/r16zEMg89//vNcfPHFPPjgg3V9fW3APUhSJzE9\niRn7noqentdJpnOevC6rVq0C4Morr+T888+vax8U0Abck+gQijfRJyPvEs9kPX9dent7mTNnTl1f\nUxfyeBBdSu9NdIWsd0mk87V74Le21WcRt5ZvsfDoo48ya9YsvvCFL9TnvQpoA+5BkpkcLeGA5z2K\nqUZcyzs9S3IiKpQKhrdeXHbZZVxyySV1f10dQvEg8XSWjkiDNhQew8xNaHmntxBCEE9nPXm/vPDC\nCwBs2rQJIerf31974B4kkcnTHgnqo7rHSGRyNAb9njQUU5l0Lk9eeDM3sXLlSkcMt4n2wD1IIp2j\nrTGoDYXHiKdzdOiN1XMk0/J6eFGF4jTagHuQRKYQQvGgRzGVSWZytOnr4jnimSzgTQ/cabQB9yCJ\ndI72SHBK9nbwMol0jvbG4JT09LxMIp3DMOx74E6GNuwwkXVoA+5BEpk87Y1BUrppkqeIZ3J0NAWn\npKfnZRKZQsjRxnUJBoMkk0kHV2UNIQS9vb2Ew+Gavl8nMT1IIp2lPdKgY60eQ56MGkhlh91eiqaI\nRDpHR6SB/nja8vd0dnayd+9e5xZlg3A4zLx582r6Xm3APUgikyuoULQH7iWSmRzTmvTG6jXM++Xo\noHWPur29nfb2dgdXpQYdQvEgZgxcGwpvkchIT0/rwL1FcW7CK3FtVWgD7kGkR6HVDl4jns7Rpk9G\nniORydEcDuI3DDI5bcA1LpIpFCW0hALaUHiMZMED1ycjb5FI54gE/YQCvil3zzgSAzcM43NAN/A+\nYBMQA24XQhxz4v1OJJKFar9w0K8NhccwC6wyuTxCCAzDcHtJGgoVsg3ynklm8rTUJuiYlDjlgQeA\nS4GZwIzCY1ePfZJhGB81DGOdYRjr1q9f79BSJheJtPwwSm9CG3AvEU/niDT4Cfr1tfES8VH3zNTy\nwJ0y4EHgPmAR0Fh4bJy7IoS4QwixWgix2mx8PtUx+22EAn6tA/cYyUyRodD5Cc9QfGqdaglmpwx4\nKxBCGvFuIAz82qH3OqEYMeBB7eV5jcSo8JbeXL1CPC2vS8MU9MAdiYELIW4u/PNuJ17/RMY8DjYU\njuk61uodEhkZQtHhLW8xNgY+ldAqFI+RLHgTPp8xYsQ13iCezhEuqB2SOrzlGcx7RsfANa5jehOA\n9vQ8RC4vyOTyhAI+mZ/Q18UzmKfWcNA/5XIT2oB7DPPDCBTi4FPLo/AqZqLMMAzC+rp4ikRGq1A0\nHsFMlAEFJcrU8ii8iikhBH1dvEZCq1A0XiE5yoDrEIpXSGZk/BvQCiGPkdAxcI1XGOXpabmaZ0iM\n2Vh1EtM7mOog7YFrXCeRLvL0tAfuGRJFG6tuc+AtEkXqIK84PJsODrDlyJDj76MNuMdIjlWhTDGP\nwqvEx22s3jAUGm/qwO974TDP7Oxx/H20AfcYOoTiTUZvrNoD9xLm6chLG+tAIk1bY9Dx99EG3GMk\nMqM9Pa94FFMdHQP3JkIIktkc4YC3PPCBeIb2SIPj76MNuMcYayi84lFMdYr1+VOxYMSrJDN5Gvw+\nfD7DU3UTA/E07RHtgU85EmP1xvqo7gnGb6z6ungBU4ECEA54yANPZOjQBnzqYWpaQeuNvURy3HXx\nhqc31Rm1sXrofhmMZ2hr1CGUKUcikyM8SoWiDYUXGFuJ6RVPb6qTSGdHJZe9kJsQQjCYyOgk5lRE\nh1C8SfHGqnuheIdEOj+md5D798twKkso4KMh4Lx51QbcYxQfCcMe+UBqxrY40BurVxh1v3jEAx9U\npEABbcA9x7hmVtrT8wSjchO6wMozxNNZGhvkXBqveOBSQuh8+AS0AfccibSuxPQi8czoNr9JvbF6\nAnkykmYs7JE5sgMJNRJC0AYcgHxecM+6A24vAyHE6EIej3gUmrEeuNaBe4W4B1VbA/EM7QoUKKAN\nOABDyQx//8tNCCFcXUcmJ/AZEPTLy6JDKN6huJReJzG9g+yDIkMoshLT/esyEE/Tpj1wdURTWfIC\n16VhxV4e6BCKl4ins55MYv7dLzby4oEBt5fhGuNyEx64LtID1wZcGbGU3LWHUxlX11E8DxN0Kb2X\nSGTyo3ITXvD0AF4+NMSRgYTby3ANmTMqxMC94oEnMnRoFYo6oqkscNyQu0WxAgXMboTuexQa78oI\nu6MpEh4wWm4hS+kLKhQPeeA6hKKQ2IgBz7q6jmJJFOgQipdIpMfGwN2/Lvm8oC+WJp6e2gbcTPoH\nfAZCCLI5d6/NQDytQygqMQ131GUDXiyJAh1C8RLxdJZI0PT0vCFX64+nyeUFialswIti4IZhyDCK\ny5vrQEIX8iglVrgB3PbAi8uCQY/uAuiPpd1eAiAT3OFCrFXqwN2/Lj1R+buZ6h54ZGzeyOXNdSCe\nVtKJELQBB7zjgUulw5gQigcMhZu8+9tr2dsTc3UN2VyebF72nQZo8PtIZ/Ouy067h1MAUzsGXjTq\nDvCEBz6Y0DFwpXgqidkwJok5hW9OkAmhvri7XriZXDYMAwCfz6DBA5trT7RgwNPuOh5uUlK55eI9\nI4SQSUwdA1eHV5KYpWPgU9sDH05liSZdDm0VFYuYeOHa9ERTNIcCUzuEkh4dQnF7rFo0laUh4CMU\n8Fd/ch3QBhxpuFvCAQ+EUI5LosAbRsJN0tk86WyeYZcNeDKdH9Eam3ghkdkdTTF/WoT4FD6ljZPe\nupz4V1nEA9qAAxBN5ZjdGnbdAy+WRIEupTevx3DS3QKreCY7ykiANzbXnuE0C6Y1kpziHnh4TO2E\nmx74oEIFCmgDDsjk4azWMDGXY4nJMaX0Qb9BNi/I5d1NlrlF1CPJZakBHx1C8UI/lO5oigXTIlM7\nhFJKheK2B64ogQnagAPSQMxsDRF1OYkZHxPPMwzD9Q+km5iGe8gLMfDg+BCK271zeoYLBlyHUEb+\n73YMXGUrWdAGHJBH9VleCaE0jD2qT93WpV4JoSTHGAnwxmDjnkIM3A0Vyhcf3MLanT3K33cs8bS3\n+gf1KxpmbKINOFI+OKsl5P5RvZSh8ECsFeAbj+5QfmMMmyEUlz3wsUYC3G9zYJbRz+uIuKID33Z0\nmD297urzc3lBJpcnVDR7Mhx01+EZjGsPXDkyhOIBD3xMCAW84ellc3m+/scddA2llL5vLJXF7zNc\nV6HIcu2xMXB3q2T742maQgHaGoOulNIPJzMMJrxxMjL1+eC+Bz4QzyirwgRtwAGIpbPMag25b8BL\neOBhD3S+OzKYJJcXypNl0WSWWS0h19v8ymEOY2Pg7hqKnmiaGS0hIg1+V5KY0VTWdQM+NmcEXoiB\nq5vGA9qAA9LTm9kSdj2JOVYSBQUP3OUY+IH+OIBylU40lWV2W9gbIZRxoS13DUVPNEVncwPhoJ9E\nJqe8rD+azDLkAQ983P3iAQ9cVRk9OGTADcMIGIbxKcMwXm0YxtcNw7jZMIwVTrzXRJE9LWBaU4Mn\nPPBxIRQPaMEP9smBAXHFG1w0leWktkb3QyhlKzHd9MBTdDaH8PsMGvw+5ZvJsEc88HEnVrc9cIWt\nZME5D/wWYCnwGWAV0A1cP/ZJhmF81DCMdYZhrFu/fr1DS6lMLJWlKRQg0iANpZua60S5ZJnLIRTT\nA48r9sBjqSwntYU9IiMcbyjcvC7dw9KAA0Qa/EoTmUIIT4RQSjs8LnvgiQwdTZM/hLIK+A3wRuCc\nwmPG2CcJIe4QQqwWQqxetWqVQ0upTDSVpTkUwDAMIg3ultOXV6G464Ef6IvjM9S3LR0JobgdA0+X\n0oG7G9rqjqaY0WIa8IDSzTWWziEE7hvwEiFH9z1wtaX0gepPqYnbgQuAu4AQMA/4iUPvNSFi6SxN\nIfkhaAr5iaWyyjqJjaW0B+6+DvxAf4JFnU0uxMBzzGgJkc1JuVjQ707KpqSMMOjuXMye4TQndzYD\nsipUpRLFzEkMJdw+GWXLnFjduS5CCAYTaVonuwEXQvwO+J0Tr11vzBAKQFMo4GocvKQHHvSRdDsG\n3h9n1cIO9THwZIbmUIDmcIDhZJZpCo+mxZSOgbsbQumJpuhskb8P6YErNOCpDNOaGjzggec9pUKJ\np3MEfL5xpwInmfIqlGgqR1Ph5mwOuRdCEUKMa2YF7h/Vk5kc/bEMSzqblXvgsVSOplBAXhcX4+Cl\nKjHd7oViJjEBGhXHwIeTMjcxnMyQdzNn5DEVSr/iIh7QBrzggRdCKA0B14Y6pLJy4ovfNzpV4Lan\nd2ggwUntYVrC6vtODxfyEy3hIEMultOXkxF6Komp1APP0h4JEmkIjFTLukEiPb5LpJseuGxkpfaU\nqA34mBCKWx54qfg3uJ/EPNAXZ15HIxEXwkuxEQMecFVKOHbqC8jr4lYM3Cyjn95shlDUFvNEk1la\nQkHaGoOuasG9pkIZTKhNYIIFA24Yxu2GYZxjGMZnCpK/m1QsTBWmkQBoCbsXA4+XOKaD+70dDvQn\nmN8RocmFir9oKktzOECLixsrjJ58biJbHLjk6SUyNIUCI1NfwkG/UhXKcFJel9bGoKtx8EQ6X+K6\nuO2Be8yAAz3AqcBXgceQipIThlg6V+SB+13rCV7ZA3fPgB/sizN/WoRIg/rNLTrKA3fPUCRLeOBh\nF9VBZhWmSaTBr/Q0YIa22hoDrhrweCZbonunmxurN2PgZwB3InXd30Rqu08YokUeuJshlFKJMnC/\nmdXB/gTzOhppCqn1wNPZPPm8IBTw0RIOuuqBl4yBu3hdiuPf4IIKJSlHELodQkmmc0RKFVi5FNoa\nUNxKFqzJCK8DThNCbDAMY2bh/ycMsZRsmATQ7IKXaVLKSID7ybID/dIDF0JtLxQzN2EYxoiM0C1K\nx8DdO6r3FBXxADQGFcfAUxk6m0O0uR1CKZObcM0Dj6dHbawqsGLAXwu8yzCMILKacglwiZOLUkk0\nlSVS5IH3xtKurKPUhxHclxGaScy+WFqpDjw6JjfhpqEoOdDBxWTZWA+8scFPn8LPbTSVZXFnM61h\ndw14vOSoO7XhpGIG4hlOmdms9D2tGPAfAUcBU/A5w7nlqKc4iemmDrxUogzcPapHU1kSmRwzmkOk\nMnmlHvhoAx7kUH9C2XsXk8nlyQs5n7QYN3uhmK1kTSINfg71q/uMDBWSmG574OUnJbmXXPZiCOVr\nwP8AGaQH7k7TEoeIj0piuhdCKVUWDO6GUKT3HcEwDBpVS9UKChSAlpB7IZRERsZZi4cGgLtH9Z5o\nisWdkZH/Kw+hJLO0hAK0RYIcPZpU9r5jKd8/3yUZoUdVKF8GuoB+oA94yNEVKUZ6esd7objngY8v\nCwZ39cYH+xPM72gEZJGTSqlatEif76YKJZkeP6cUzBi4eyGU0R54gERG8enIAx54qR41Qb9BNi9c\n6SraH0/TobiQx4oH/nVgU9H/T3doLa5QXMjT7KoHPr4sGNz3wOdPk55eOOgjnc2Ty4tx1aJOYHp5\ngKsqlLLJZReP6sVl9ACNDT7lHnhzSOrA3Wz1WyrsaBjGSH4i0uBUr77SDCQ86IELIf4W2IHsKrhT\nCPF3jq9KIbGiXigyhOKOV1WqLBjcNRQH+uPM75AG3Gy3q8oLL25x0OxyCKVkgVXAPbnaOAMeVN3M\nytSBeyAGXuJ05EY5vRCCwXhGeSdTK5WYXwSeBr4NPGkYxpccX5VCol5JYpYoCwZ3Y3oH+qQG3ERl\nyba8LvJmcLOUvqw6yKWNdWwZPbhQyJPM0Bp2v5S+1ExMcEchlMjk8PlQ2okQrMXAM0AY8AMtQNTR\nFSlECEEslSUyEgMPuFiJmS8dQnFxJubB/uMhFFCb5C3OTbgdAy/lgTf4pQFXPYtyIJGhMegfKaMH\ntRurOY2nKeR3XUZYLuzohgcuBzmob3dsJUh0OnAbEEca8pmOrkghqWwew2DkZmgK+YkmswghxqkO\nnKa8CsUdT08IIWPgHccNuEpDEUtlRxJ1zaFAYQqM+utSKlEG4CvMokxlS2+8TjG2iAekwVLVjTCe\nzhEK+An4fSMeuBvXBSpIb13wwN1oJQvWPPCbgXbgSqCz8P8TgmIJIUhDbhi4YjATZY+D7oRQBuIZ\nfD5j1ITtJoWVqsUqlIDfRyjgI6a4mRaUD6GAO2GUnjFFPGBurApPRgV5Z0PAR9DvznXJ5PIjaxiL\nGx64GxJCqOCBG4bRDgwiwyZfKPrSamC3w+tSQiyVHUlgmphhAtWxrHLJMrc88OIEpklEYT+UaCo3\nkpsARoY6FD+mgnLXBYo3V3U3bnc0RWfLWAMeUDbQYbhIHQSMJDK9dV18yhPMAwl3QiiVPPBNwD8B\nLwIvFP35mYJ1KaE4gWni1lCHeIkBreBeDPxAX4L50xpHPdbUoC5HYI5TM3ErDl7umA7utDnoiaaZ\n0Tw2hCI3eRXTcaIp2cjKxK1EZqKMPh8KHrhip8eNVrJQOQb+HqSnHQfWFR4zkB74CUGxVM3ELSVK\nMlNatyrjrOrjvwf6ZRVmMZEGv7J+KLExHricyqP+ulQPoajd7McW8YCUeIYDcqxak8OecDR5PIQC\nuCYlLBdyBLc88PSocKMqyl5tIcRzID8cQognCv82gCZFa3Oc4jiriVs9wcsdCQN+Hz7DIJMTNAQU\nGvC+OKeOacyjUqUzPObatITd2VgreeBhFzoSji2jNzETzE4b8OExJ6NWl3qCVwyhuOSBuzF0u1IM\n/CzgdmCxYRjvL/rSTE6QoQ7FRTwmzS5V/Um1Q+mIlplVL5WwcYqD/QkuWz5acKRahVJ8VHcrhJLM\n5Ogoc2O6ksSMjk9iguxIqEILPlykzwdcm8pTLuQILnng8TRLOtX7tmUtghDCjIEfBjYW/rwAfErN\n0pwnlh7vgTeH/K6U0ydLtMY0CbnQ+e7AGA04SAPuhgoFoCUUdKWYp1wpPbhjKMoZcFWbqznMwcSt\nGHiyTOEb6Bj4CEKIxwzDuEwIMWQ+ZhjGHOeXpYZYUbGIiUqpXDGVjoRhxUqUfF5wqH90FSZItUNP\nVE3f6bEJ5uawVKGopmIM3IU+NaVi4GB2JHT+9zP2urgZA/fSxupGK1mwVshzjWEYH+W4VqoVOMW5\nJakjVjIGHiDqkgqlUkxP5QeyO5qiORQYl1SVY9WcNxLF49RMXFOhVNpYFScxS5XRmzQ2qCnmiaay\nTC8KKbU1BtnTE3P8fccSz1RWoajeWAfjGTqavFnI8wlgK7KE/iXgPkdXpJBoanzSx42OhPm8IJ3L\nEw5WioGr+0Ae6Iszb1qpRFlASdFGrFAsUqy6cU2FUtHTU2soBkuU0Zuo0oIPJzPjVCiuhFBKzMM0\nccMD74+nPacDN/kf4G7gCeAp4GxHV6SQWCkduAsGPJnNEQ6MHxpgotqAF/cBL6Yp5Ceu4HcTLVFg\n1eKSvDNRppQe1Pdq7ylRxGOiauDGcNIjIZQKoS3VMXAhhCutZMFaCCWO7Al+KfA44N6I9DoTS2XH\nJUKaQ34O9qs1FOX6bZiEFLcuLe4DXowqD3xssQi4HELxSCl9d5kEJkBEUT+UsdfGTRVKpY1V5f2S\nzOQxUN+JEKz1A/++EGKZEOKwEGIpcLGCdSmhZCWmCx54pWM6qDcUpcroQd1UnlL6fLcm05eau2gi\ne4IrNOBlEphgeuAKrk1ytIzQVQ+8bG5CrT5/IOFOIyuorAPfUOZLncACZ5ajllKFD24kMcs1pjdR\nHwNP8Nazx4uNIiE1lZilDLhbU3kqJ5fVJjFLldGbNDb4iSvwOkuV0g8m3NlYp5fT5yvuRtgfyygf\npWZSKYQSBXaVePyEUKBAGU/PBQ+8kpEA9R0JS5XRg7peKNExDZPAvaEO1WSEKj09qQEvbSgiwQBJ\nl2LgbrSUjaezJfM04I4HrnoSj0klA/4GIURi7IOGYbQ4uB6llE1iKi6lr2QkwEyWqflAZnN5jg0l\nmdMeHvc1VR54qR41cjK9O2qHSrFWlWvqGU6xcGFHya9FGvwcG3J+LcPJzCgPPBz0Q6EFs8oYcLkB\nKKDeA3erlSxU7oWSADAM49ExXzoJWOHkolRRupmV+sn0leJ5oPaofmQwSWdzqLRULSj7xDjtbUXH\nlGuDDKGo9sCFEMSrxFpVJzHLxcDDCkIox6fxjDYbZhxcpQGvFHZU74G700oWrMkIpyH7gg8CMWCc\nVz5ZmTRJTIXJsnIJTJCNtYJ+5+Px0RIVsuGgj1xekFZoMDM5gQEE/RX0+cpDKO6pUBKZ3MgQh2Lc\nSGTG0+MVZCaqPXC3yujBmozwo0Cy8G8D+EvnlqMOIQSx9PgWrk0h9SXblVpjgloVysG+BPOmlY4t\ngpqBF9Fklpmt41umNhekhNPLGLF6UzW0FfSRVJnEHE6X1YFHFFRijlWgmLSG1XckLDcPE1zwwONp\n2j2YxDR5DijuFP+CQ2tRSiqbx28Y4zr8NTUEiGdy5PMCn09NUiZRoSwY1CYxK3ngYPbcyDHdwTXI\nJmPjO7uZLWWVGXAPnYyEEPTGUmWVFypUKMOpLK3h8SajrTHIYFyxAa/WC0WxB77IhU6EYM2AfwJ4\npfDvNLIr4aSnVPwbwO+TzfHjmZyyMVGJCmXBID+QquLyB/riXHRKZ9mvq+iXPlbpYKK6I2E1D1xl\nLxSzjL6c19kY9JNQcV3KGXAXPPBSA1DAJR24SyoUKzHwHyDL598FnAeocX8cJlaiD4qJ6ji4FRWK\nKk/vYH+iZBWmSUTByLlS6iBQX8xjyQNXFNrqHi5fRg/yujhdSl9uJqlbBryyB64yhJJxZRoPWPPA\nnwTmA93A1cD1VBmrZhjGCqThfwTYC2wXQjw+gXXWnVIJTBNTiTJL0VoSVbx92Q9cjad3sEQb2WKa\nQgpirWWuTavicnpr8k4116VSGT0UuhE6vJZoqvTw4rbGIEOKJZ5yJmaZ5LLC+wXk6citQh4rHvgL\nQojZQogzhRALgIcBDMMoKUg1DONU4GRgEXAtsB34qzLP/ahhGOsMw1i3fv36WtZfM6WGOZg0hxV7\n4NWSmIo8inxe0BNNMbNlvAbcJKKgmKdUl0hQLyX0kgdeqQoT1CQxy4VQ3OiHkighQDBRrQ7qj3uw\nlL6IVYZh/AEpIfQhR6wtQ45Ve1WJ538SmAvMKHr9ktlAIcQdwB0At956q/MjtYsopWc1aWpQ2/mu\nuqFQY8AHExkiDf6Ko9uaFPTciKYy45pZgfqB09b0+YoM+HD5Kkw4nlx2kmgqS2u4hAqlMcgrR4ZK\nfIczCCFk4r/M5zQU8JHO5ZUJEQbi7unArRjwZmAIMFfYB3QA7aWeLIT4tGEYi5Bhll8BZwHfm+hC\n600slaWpjNcry+nVHcGqH9XVdCPsjVU+pgNEFPxuyuUnVHckrJrEVKgOqlTEA2oGOlSKgavsCZ7K\n5gn4fATK6PMNwxhxeipdv3qQzOQQULaXv9NYMeAXIFvKLgD2CyGSAIZhTCv3DUKIvcgQimcpNY3H\nRHUSs2ovlKBPSX/jnmjpaS/FKPHAyxiK5nBAqVwtkc5Wvy6KjuqVyuhBep3ZfJ5sLl/WsE2U4VS2\n5OR11UnMas3f4LiU0GkDLr3voNI+MMVYudKXAseALcBRwzDeCSCE6HNyYU4TS5VPHDYpPqpb6kao\nwgOPppneVMUDd1iFksrmyIvR49RMWsJBhlWHtjxwXaByFSZIr9PpqTzl5J2qDXi10BaokxIOJNKu\nJTDBmgd+JXANkAEagbcA9zq5KBWU04GD+sn01T6Qqnpu9MZS1T3wkJ/u4ZRja4ilcuPGqZm0qpYR\nZvLeSmJWCKGA/Jwk0jlaSsSp60GpQRug3oDHqyT9QV3/oP6YexJCsOaBh4DpyGHGLUB5icIkIlpB\nheJKCMUDKhQZQrHggTsYa42VGKdm0qy4I6EVGWE6l0cI5/Pv1XTgIJUoTiYyo8kKMkKFPcET6fJl\n9CZhRa1+B10s4gFrHvhdwPcoLZF2AAAgAElEQVSBmUAX8GFHV6SIWCrL3PbSeufmUIDDA8mSX3OC\nSlNfQF0pfW80xbLZlbsFOz0XczhZ2suDwlAHpTLCbMWwhc9nEPT5HG+lWq2M3iTisBZceuDjjVWk\nwU8mlyedzVdUMNULSzFwRR64m42swJoHPgvphRtABFjo6IoUEUvlynp6yisxLcVaFYRQLMTAG4MO\ne+AVTkaqhzpU88BBzemoWhm9idODjcvFwA3DUKoFtxJCUeWBy2HG3o6B/y1wC9L7bgM+BHzbyUWp\nIFohBt4UChBVONQhns4SCVaqxFQTQrEaA3dShVJOgQIuhFDSlWPgUOzpOeeFWQmfgNkPxVkDXu50\nZMbBq8Xp60GlToQmymLgLhbxgDUD/hDwIscLeVYYhnEWsFII8UMnF+ck8UqVmIqTmMlMvmxZMKgM\noaQrFouA8yqUSi0OWlWrUDJZCx648x0Jq5XRm0QclnhWvDYKPfBqIUdQ54F3D6c4eUaz4+9TDisG\n/G+Azxb930B65QKYtAa8XLk2FGY/KjIU2VyebD5PQwXtrrokZqpqCMVxD7xSj5pCiwNVFXbVKmRB\njadnRYEC0OigjLDcNB4TlcU8XlKhdA+nmKng1FEOKwb8r4GXSjx+Vp3XopRyHe9A7WR6sy1mpUKA\nUMBHpqB2cKpgIJ3NE0/nqg5nbXLYA69UYOX3GYQLY92cksoVY0VvrGKwcfdwqmIfFBMnp/IkM3kC\nvvH9801UNrTykgqla6hy7yCnqWrAhRDfKPOlJ+q8FqVUMhSy54aiD6OFeJ5hGCOjzJxSO/TF0nQ0\nNVT1bJ1WOpRrmGRiDnVQY8Crl2KrOB31VCmjN3EyiTmcylT8nbc1qpvKYym5rMgD7xpOjpsepRJ3\nCvg9QDSVpbmiCkWNB55M52msEP82cVqJIsMn1bPpTit0YiXmYRajsiNhIl09Bq5iqEN3lUZWJk62\nlK2UwAS1U3mqDUABNSejdDZPNJVlmosqlClpwIUQxFJZImUrMdWV0sczlRUoJmGHexz3xtKWEmVm\nOCebc+bmKDWRvhiVShSrIRSnk5jVyuhNIkHn8hOV1EEgE8xTzQM3c0aqRi+WYkoa8GQmPzJhvRTm\nBPSMQ0aqGNmYvnpYxOmjem+0uoQQZDjHnBvqBJXknaBWC55IWw2hOG8orIZQEmknN9YqHrjCJKYX\nYuBdw9aui5NMSQMeS1f+MEojpUZKKL08iyEUJz1wC0U8JpGQn7hDIaZy/TZM1Brwyt0IwZz+4nwS\n05qMMEAi48zvplpuQqUBT2a8oULpGkq6qkCBqWrAq3h5oC6MUmmySDFOx/R6LBTxmDQ5OJWnUi8U\nUDfYeGRoQFVPz9ncRD4v6ItVb/ML0NjgcyyJWW1jVWnArcg7wwpCW13DKVcTmDBFDXi0ipEAdYlM\nK3FWcL4a00oRj4mTHrg1FYrzhiKdy+P3GWXDbCayV7tznxM5JSlAKFD9M9IYdG6w8XAyQ0ulGHhj\nkCFFJ6N4pnrYUYkHPpxihosSQpiiBrxSL3ATVT3BrWhawfkQSl/MRgjFSQ+8SnhL1WR6K14eOJ/E\nlFWYFjfWBr9jQ5ajFkIoqgp5khZUKCpi4N3DOoTiCpU04CbNihpaJSzE88D53tNWk5jgbMl2NbWD\nKhmhFaUDOL+x9thIlDnZTraaOkhlKb1XVCiyiEcbcOVUy6iDLBlXYsCrdCI0cV4Hbk1GCM5WY5ab\nh2miKolp1QN3etiG1T4o5lqcK+SpEtoKyTJ+p+SlxcQtJJeVeODRFDNbdQhFObFUtqrXqyyEYiFR\nBqbawbk+F1Y6EZo45YGnsjkEpcepmbQo0oFLL89KctnnWNgCrCtQoFAl6+DJqFIM3OczaA4FlMTB\nk1YqZBV54FpG6AKxdGUvDxSGUCw05gFndeDxdA4Dw5IaBpxL8Jred6V+L8pCKGk78k5nT0bWQyjO\nNbOqpkIBdUoUL8zEzOel02OlR42TTE0DbimE4uzgAhOrKpRw0LkBur0WptEX45QHXi3+Dcd7oTiN\n5Rh40OEkpsVGVuBwL5Qy49SKUWXA4+lsVWfD8aR/PE1zKKBkAlElpqwBt+KBq1KhWFY7OOTpSQ24\ndU/Cqc3NSm5CqlAUeHmWY+AOJzGjKTpbrG2uTg50qCbvBDVKlHxekMrmK4bZwHkP3O0uhCZT0oBH\nqzRMApRVYsZtqR2c+UD2RtN0WmhkZSJjre4YcGVJTMsxcIcLrKIpZjRbMxSmN5h24HMSTWVpqaBC\nATUeuGm8q/UfcdoDd7sLocmUNOBWPHBVScykLb2xUyEU6wlMcG7ghZXromoqj70YuMNJTIseOJj9\nUBzaXKt44K0KWspaCZ9AQR3kpAfugT4oMEUNeDRVvXxdpQ7cWqzVQQ88lrYVQomEnIm1VpOqgTSY\n8hjtbH7CcoWsg6GtfF7QH7deYAXO9GsXQljKT6jQglu/Ls6rg3QIxSUsJzEVlNLH7ejAnYqBW+wF\nbuJUL5RYhR7tJoZhKAmjWA2hhB3cWAcSGZpsJsoaHWgpm8zItgLV1qEiBp7M5AhbOhk5rM/XHrh7\nyIHG3tCBWxnQCs4ONu61UcQDBRWKA5tbtXJtk+ZwgKjTBtxGaMspT8+OBtyksaH+/VDkNJ7q10VF\nDDxuufmb3FiFEI6so8sDZfQwRQ24JbWDwkIeq0dCp2J6dop4wFShOCAjtBADBzUdCWWFrAVPz0EP\nXCYw7RkJJ0IoVsInoGYuptWN1eczaPA7d228UEYPU9SAVyvXBnWl9FYmbIPDMXAbvcDBuZ4bUulg\nwYCHAww73JHQ1sbq0MlISghrMOB1vjZWEpigxgNPWOhEaBIKOuf0yFayOgbuCl7SgSctT+RxLoTS\nY6OVLBS6EbqkQgE1UkLpgburdrA6C3Pseuq9udrxwB034BY6EZo4dc8IIXQIxS2EEMTSWZos9EKJ\np3OOxdAAcnlhqxLTCb1xPi8YiMuJ9FZxU4UCasrpvaB26LY4Sq0YGUKp7+9mOJWtOJHeRMVcTKuq\nLXDunhlOZfEbhiVnw2mmnAFPZHIE/T4CVRr1B/0+Aj7D0SKNwwMJZraEqg4NAOe8iYFEhuZwwNIa\nTMzhufXe3KpNpDdpCQeIOh1rzViMgTuodqglielEeGu4SiMrExWT6a02fwPnwlvdHgmfwBQ04FYS\nmCZOh1H29MRYPKPJ0nOdkhH22pQQAiMDoeu9ucmjenVPT06md9YDj6dzNAatqR3SOWfUDnYaWZk0\nBgP1j4EnM5ZORq2NQaKpLPm8c6dWq83fwLly+q4h95tYmUw5Ax63kMA0aXK4mGdPT4zFnRYNuEMJ\nmZ6ovSIeEyeUKNUm0pu0KKjGtNJyGKTaIeiTRrze9NhoZGXS2OBzJolp4Z7x+wyaGpzdXK2qUMA5\nD7xrOMkMD5TRwxQ04FalauC8Flwa8GZLz3UqhNIbs58oA2e04Fb6bYDzSUwhBPt74yyYFrH0fBkH\nr78BtzPMwSTSECBe55i81dwEmLMxnQujWO0dBM554LIKUxtwV7AaZwVodlhKuLsnxhKrHrhjIRR7\nEkITJ6oxY5Y9cGc7EnYPp2gI+Cwndp0YHpDPC/otTqMvxomOhNGktSQmOK9E8YYH7o0yepiKBjxt\nzwN3angvwJ6eKItcN+D2inhMnFCiWNUbO+2B7+yKcvIMaycjcGawcX88TYvN5DI4owO3msQE5w34\nsEVJIzgZA/eGhBAcMuCGYVxqGMafDcP4lWEY/2kYxs2GYaxw4r3sErUZA4861A8llc1xbCjFvI5G\nS8+XgwPqv5Yem42sTJoaAnXtuWF6SqGAtRi4k6GtXd1RTp5pw4A7UGRVS/gECkMd6l2JaSPx39oY\nYMBBJcqu7ihLbCX+HfLAT/AY+DPAGwt/zgO6gevHPskwjI8ahrHOMIx169evd2gpo4mlqmvATZod\nKlgB2N8bZ157o2UPy0kP3E4vcJNIg7+uzb6sFouAqUJxzkhID9yakQBn+qH0DNtXoIAZQqlzctli\njxqAOe2NHBqI1/X9TYQQ7Dg2zKmzWiw939kY+AkcQhFCpIHLgfMB81M4rgO7EOIOIcRqIcTqVatW\nObGUcVit9gNnVSi7bShQAAI+g7wQdZ/63TsBFUo9PXAr7Q1MHA+hdEc5xY4H7sDm2lOjBx5xpJmV\n9c11yYxm9vTE6vr+JocGEjSHA7Q1WovHO+mBe6ETITgXQvkL4D+Au4AXgHnAT5x4L7vEUjkbnp7f\nsaO6HQkhyDaqoYC/7nK13hoSZVDwwOtoKIZT1WcumrSEgo52I9zVFbMZA6+/oai1XakTczGjFrsR\nAizpbGJ3tzMGfMexKEstet/gjAeezORIpHN0RKxtIk7jSC2oEOLHwI+deO2JEktnmWYxZNAUCtAT\nTTmyjj3dMc6a32bre8zS4Ih9e1uWnmiKzhpUKFJGWF8P3PLGGpbJ5XxeVB2tZZfhZIbBRIa57dZy\nE1Doh+IZD7z+4Rw7KpTFnU2OeeDbjw3bMuBObaydzQ0YRn0/d7Uy5VQo9nXgziQx7XrgUH8teCqb\nI5nJ0dpofx+PNNR3sHHMhtbY7zNoDPqJOqAQ2t0tr4udjUG2+nXGUNilsc7NrIQQDCetyTsBZreG\nGUpmHDm5bjs2zNJZdpLL9ffAu4ZTzPBIGT1MQQNuTwfubAx8icUiHpN6V2P2xdJMa6rNm2gK1dcD\nH7axsYJzQx12dtmLf0NBIeSACqWWEEq9e6Gksnl8PsOSOghkZeqi6U3sdcAL33EsajmBCU554N6R\nEMIUNeBNFtqEgnNJzKFkhng6yyybUqR6J8tqLeIBZzxwq1pjgGlNIUfCW7u67WnAwZlhGz02pySZ\nyKHGddxYbWjATZbMaGJ3nQ14Pi/Y1R3lVBubqxMx8C4PVWHCFDTg9kIoziQx9/bEWDS9ybbnW+8Q\nSk+NRTxQ8MDraCiiSXse+NJZzWw7Oly39zepxQMPO1CJWXMSMygn8tSruZbV4qpiFnc2safOicyD\n/QnaG4OWY/HgXAzcKxJCmIIGPJ62LldzqhuhnS6ExTjhgdfi5YE51KGOOnCbIZRls1scMeCyiMd+\nbqKenl6u0KPdarK9mIDfR8BXv8+JTGDaNeDN7O6J1uX9TbYfG2bpbOvhE3Bm2EbXkHckhDAFDXjU\nRgzcqRDK7m7rPVCKqXcMvDdmv5WsSb0rMa2OUzNZPruFrXU24JlcngP9CRZNt2vA6+vp9cfTtDYG\nbZfRmzTWsZx+OGld3mmyZEb9lSjbbCpQwBkP3CuTeEymnAG3U8jT4pAKpRYFCtQ/hFJrEQ/IXij1\njoHb8cCXz26tuwHf1xvnpLaw5YEBJvUe6tBdQxvZYiJ1LKeXRTz2NM9LCiGUevZI33Fs2Fb8GxyM\ngXukjB6mpAG32wslU/dm/bUb8PqGUGQv8Al44HVWodiJtZ7UFiaVzdFbx0Tmzq4op9hMYIIZA6/n\ndUnR2VK72L+eHngtIZT2SAMBv0FPNF2XNQBst1nEA0554DoG7hpCCOJpeyqUtsYgB/sTdV3DxAx4\nHT3wGnuBQ/3lanbknSArU5fXOQ5ut4mVSb3nYtYySq2YeraUtdPIqph6FvTk8oLdPbUkl+ufm+iP\n2RsA7jRTyoAPJbOEg378Noo0zp7XzsaDA3VbQ3dU9ppur6Gcst5JmYnICOvdC8XqOLVils9uZUs9\nDXiNHnioztelJ1qHEEqdrk0tKhSQicw9dUpk7u+LM6MlZHuIcL37tPdGU7RHglXn6arEOytRwIb9\n/Zw51175+qtm+9h4oH4GfE93bd431N/Tq7UXONS/F4rVcWoj7HuWFTPDbDs6VLc11KJAgfqfjHqi\naTonkChrrONUnmG7IZR9z8LQ4bpqwbcdHWbpTHvhE4BwnXMTsomVd8InMMUM+NqdPVx0Sqf1bxg6\nwoefewMH9u6s2xpqDZ9AfSv+hBCyF3iNHngo4COXF2Tq0Fwrnxd0Daesy+ZyWfjJu7kg9mjdQihC\nCHZ122tiZVLvXigTDaFEgn6SdVSh2Crk+cPn4fEvsbiOTa3stJAtJhSsf2jLSxJCmGIG/JmdvVx0\nynTr37DvGXz5NGcc+03d2rhOyIDXMYkZTWUJ+gzL8wXHYhgGkTr13dhydIiWcICT2iw2kDr2EuTS\nzN9zN9uPRcnVYQr6saEU4WBtoa16n4x6aiyjN6lnR0JbIZRUFLq3wiu/5pTWbN1i4Nu7orZ6oJjU\n3wP3loQQppAB74ulOdAX56x57da/ad8zcNa1XON/jJ3HBuuyDjtzMMdSz6P6RCSEJpE6VWM+ub2H\n1y2dYf0b9q2Fs68jED3CqxoPsb9v4gME7I5RK8YJGeFEEmX1nMqz9cgwC63q4g8+D7PPhFOuYPGh\n+9nfF6/L5rqjBg041N8D7xryVhk9TCED/uyuXlYv6rBXHLFvLZx/E6lQJ10vPFCXddRahQn1nb3Y\nG6s9/m3SVKdqzCe3d3PxqTYN+OKL4dzruSH0WF3i4LUqUEAainSmfgndiXrgkTpN5TnQF6cnmuJs\nq07P/udgwQVw3kcIrv9fZjQ1cGiCCq5MLs+enphtBQrI++XK/NOIvU9PaA0mXuuDAlPIgK/dZTP+\nHeuBocMw60wOLnkP07fdNeE15PKC/X1x25V+JqGgj5MG1kMddOk9E1CgmNTDA4+lsmw6OMAFJ1sM\nbQkB+5+VhmLlBzg//gQ7Dx6b0Bqgdg04QNjv44u9n4EdD094HbKMPsO0CTR9l4ONJ77RP7ati0uW\nzbSu2tq/FhZeKK+NL8BVLTsmXFK/rzdWU3EVgN+Avwn8AvHHf5vQGky6hpPM9FArWZhSBryXC0+2\nYcD3rYX5rwZ/gOZV72XB8AYYOjKhNRzqTzCjOVTThxFgZmIXN+78JOx+bELrALMPig0jkR1fMFOP\nfijP7urlzHlt1iViPduhoRna5kLbXAZmrKZt128mtAaQHngtXh5A6+BmTs7tgrW3TXgdfbE0bY0T\nk6qFG/zE63AieGRLF5evmGntybkMHNog7xnDgPNu5G2Z3084Dr7dZgvZURzdRIAcDB2EIxsntA7Q\nHrhrHB5IMJjIsNxOM5x9z8CiiwBYtvAkHsydT2b9jya0jt090ZoTmABLD/yCow0L4dlvTWgdYFNC\nmE3DV5fDbz8L6eM3ZFMd9MZP7ujmdUstGgmQG+vCC0f+mz/3g5zX++sJrQEKMfAaDXjHjnv5WeDt\n0L0Njr0yoXXUQ+kQqUMhTyyVZd3ePl57qkWn58hG6FgEjYVwy1nv5dTYenqO7J3QOrbbHOIwis33\n8YjvQuJnXQ9/vmNC6wDvdSKEKWLAn9nZwwVLptsbv7XvGVgoDXgo4OdPHW8lv/6HkK/9aDoRBQrp\nGPMOPcC3Zv0bHHlRGosJ0GtHQrjvGWibB5kEfPu1cHAdAJHQxHuCP7m9m4uX2jwZLbhg5L8zV15F\nS26AxP71Na/BnCBzUi3H41yWph338Ssug9Ufhj99u+Z1QA2j1A6uk55vUVitHoONn97ZwzkL2q23\nb93/LCw4//j/w630LLqKJft/OaF12J2DOYIQsPlXPBF8LYPL3wdb7od4X83rEEJoGaFbrN3Vy4V2\n5IOJfujbAyedM/JQ06LVDBvNEwpfTMiAv3wvQ53ncsA3B1bfCM9NzAu31Qt82+/g9KvhHd+Gy/8F\nfnYtPPZFWoJiQv1Q9vXGiKZyrJjdav2b9j87srECBINB/tj4RuLPfK/mdezqirJkhr0xase/+VFE\n20J25mbD6g/BK/dNyFDYUqAIAfd8CO65Ab5+Ojx4M+x+nEggP2EP/NEtXVy2fJb1b9j37KiNFcA4\n70YuHnpAhldqpJYuhAAc3QQiz96GU0g0dMCyN8MLtZ+ghxJZGvy+mmW3TnHCG3AhhExg2ol/738O\n5q6CwPEb6ewFHTzS+CbY8MOa1zIRBQrrv0/30vdJudp5N8LmX0Gst+a1WO4FLoQ04MveLP9/+tVw\n09NwaAOf2vNxAv21FzmZ3rdlwzmwH7JJmH7yqId3z3sHLbt+C6nainp2dcdqTmCy6S7yZ75XytWa\nZ8Kyqyb0GbGlQDm0AQIh+PRG+MCvoHkW/HENV/7uYt5/+As152zyecGj27q4fLnF0FY+X9hYLxz1\n8KxTV7NPzCD9Sm0KrnQ2z4G+OEtquWc23wenv4NQMCD7obzqL+H570G+to2tazjJDDtdCF++F7b/\noab3ssMJb8B3dccI+HwsnB6x/k37noFFrxn10Dnz2/hhdDXsehyiXTWtpdY+4BzZBMPHiC+89Lih\nWPFWWPe/Na0DbMgIj20Gnw9mLD/+WMtseP89bJ71dq56/gbo2lrTGp6wrf8ueHljJhnNmb+YXc0r\n4aV7alpHzRrw5BDs+CP+s95FKpuXXStf/Vfw5+/JatEasBVCeeVXckM1DJixDC7+HHz0MTZd9Vvi\nIggP3VLTGl4+PEhLOMAiq5/V3h0QaoHWOaMe9vsMHoq8hcxz361pHXt6YsztaLQ8j3OEQviE068+\nPi1p7ipomlmzUbWdwFx7G/ic99ZPeAO+dlcPF5483d74sr3PjPMmlnQ2cyAWJLX0KnjxJ7bXkczk\n6I6mmNtusdqwmPXfh3OvJ9TQcLxg5PxPSI+ihDqkGplcnqODSWuGwvS+x/7+DIMdC97D2rkfkuXT\nNqWN6WyeP+3u5TV2pJ37144Kn5gsm93Cr/1Xwrrv1ySxrFmBsuU3sOg1+Js7CfoN0rk8zDkH2ufD\n1vvtvx42kphCwOZfw2lXj/uSv2MedzR+ROYLCvkKOzyypYsrVtgJn6wdFz4xOTD7CgI9W6Bnh+11\nbD9WWw8Ujm4CBJx0zuh5pa/6aM3JTFmFaTFH0r1NSpAXX1LTe9nhhDfgz+zssRf/Tg3LcuC5q0c9\n7PMZnDG3jc2zr4YNd9pOZu7rjTO/o9G+PCw1DC//Es79wOhS+lmnwcwV8ms2+dULhzhjbptFA/4g\nLH1jyS81Nfh5pv1qGNgHOx6ytYYN+/tZPKPJXjXovrWwcLyhWD67lXv6T0UkB+HwBlvrABkDr0mB\nsvEuOPu9wJhqzFffBH/6jv3XA/b0xpllJZl6eIMM8c06fdyXIg0B+jNBuPQW+MMttje1R7d2cZnV\n8AkUwielDfiCGR28PPNtNZ0Wd9QwRg2Q3vdp8mQSDvpJmtXLp18Nx16G7u22X9JWFebGu+Csa8Bv\nv4ujXU5oA57LC57b3WdP/33gTzJ5GRx/E509v52n4osh0Ah7n7K1lj09URZ31mAkXvoFLHottM4Z\nP5Hngk9KSaGNGzSby/PNx3by15efWv3Jw0ehb/e404hJJBRgKOuDN3xJeuFZ6w38n7BbfRnthuFj\nMOuMcV+a1Roih0HsjPdLL9wG6WyegwMJeyE2gIED0hgUNrdwcdn28rfA4EE4/KKtl3xmZw+D8TSv\nWjyt+pM33zdipMYSaZCDjTnnfZCOSgWGRY4NJdnfF2fVwg7rC9/3LCwo/RlZ3NnE78NvhI0/GyVB\ntcK2WiSEQozEv4HRHnggBOd+EJ63H9LZfHiIeR0WTs/5PGy6G86+zvZ71MIJbcBfOTzEjJaQNY/G\nZN/aEf33WM6Z38bGQ4Ow6gZY/wNba9ndE6stGbP++7DqQ0CJmZinXA65NOx50vLL3b/pMLNaw5y/\nxMKpZPvv4ZQrwF9aSjbSd3rplVIDbOPGkAlMGwZ8/7Mw/1Ul44qGYbBsVgsvz3yrDGvYUIHs640x\nt72GOOtLP4fT3i6NAmPaHPgDcN5HbEkKhRB85Q/b+Ozrl1Zv9yCEVLucPj58ArI7YiKdk7+rK/8d\n/vivljfXR7d2cfHSGdZbTgwehEwMOks7BIs7m9gw1CpDLDYlljVJCI9sRIZPzgYY7YGDlHpu+rnM\nX1hk8+FB1u7q5V2r5lV/8t6nIDKt5MnICU5oA/7Mrh4uslqibbL3mZJxVpAe+MYDA4izroE9T8CB\n5y2/bE19wA9tkJLGky8DSnQjNAy44OOWJYW5vOC2R3fyaSveNxTi328q++VRvVDe8EV46quyBUEV\nuodT7O+Ls3KBjcZiFY7pIIccvzwYhnP+Au77uOUQ167uGhKYQsDGu+Gsa0ceGndtzr1ehp8sJrwf\nfuUYqUyOt541p/qTD78AvmDJ0wiMmZZ08mUwbQms+x9L63hkiw31CRzvf1Imx7R4RhO7u6Pwpv+A\n526HA3+29LLJTI5DA/YHTMuN7R0j6xnlgYOs4F3yOuklW0AIwRce2MJnrjjVmiZ+413KvG840Q34\nzh4utJMkyyTg6EvS0yvB7NYwfp/BwUQI3vYNuOeD8lhvgZo04Ov+Vx75fPIylRxqfNZ7ZaLKQpLo\ngZeO0BFp4EIrm1o6JjezU64o+5RRk19mLIMzr4HHvlD1pZ/e2c2FJ0+331iszMYKsPykVrYcGYYr\nboV4Lzzzfy29rKzAtHldjrwo5YxFhSsNYztFRqZJQ2IhpJPLC7760HZufsMya5JK0/suYzQbC15n\n3uwE+Pp/hyf/SzoDFUhmcjy3u9d+Z8gyCUyAGc0hsjlBf3A2vPX/wS8+bOmEtLs7xoJpERoCNj4j\npvqkKLE7vyPCs7vHyG3NZKaF0OMjW7roHk5x7Xnzq79/OgZbH4Az3m19zRPkhDXgqWyODfv6OX+x\nDQ/84PMyOdhQ+oY2DEN64QcHYPmbpZf18+urHk+FEPbbyCYHZThg5QdGHjK9vFFDloON8lj43O0V\nXy6fF9z2yA7++vJTrSlydj8Oc1ceL40uQVNoTC+U1/29jLcefbniSz+xzWb4JDkkN6g5K8s+Zdns\nFrYdG5KJvWt+II/ru5+o+tI1NbHaeLfcOIt+j02hwPjOe6++SW7CVT4fv9l4iOZwwFri0IzxllCf\nmPh8xviE9/KrpBGvwLO7e1lxUgsdVgdrgPTAK5yMDMNg8Ywm9vTG5BpWvBV+/YmqxnNHVw0JzCMb\nAWMkfAJw42sXs35fPxNNniUAABFdSURBVI9vKzoJLbwIfIGqifdMLs8Xf7eFz1+1wpr4YMtvYcGr\nocWGgmeCnLAG/IX9A5w8s5m2iI05i2P6bJTinEIYBYCL/056Wr//h4rf89WHtrNoesReGe6mn8OS\nS0Z9GHw+g6DPR398TGXbeR+RapSD5cvJf7/5KJFQgIut9rbY9uDx4p0yjJu9GJkmjfgf/rHsDZrP\nC57a0WMvgXnwz9J4B8r//pbOamFnV1QO3mibC++8A+79SynnKsNgPMPze/tZZsdQ5LLw8i/g7GtH\nPfzxS07m1t9spi9WZKxnroB5q+E3nyyrC09n83z94R3c/IZl1jbWIy/K2PbsMys+rTE45tpceouU\nv/btKfs9tqsvE/2yuGr22RWftriziT3mdJ4r1sjkeJWw39ZaxqiVOJlEGgJ88R1ncsuvXiZmVg0b\nBlz5f2SorUJI56d/2s/c9kYusepsbPzZuM+F05ywBtx290GAvU/DwtdUfMpZ89rYeLAw3MHng3d8\nRyYRN9xZ8vl3/Xk/9286zB3Xr7auRc9lpedWSF4Wc/0FC/nYj9ePPq63zJJl7j99D+x6dNz35POC\n/35kB5++/BRra8jnYftDZeWDJk2leqGs+pBUjGwtXX33ypEh2iJB5k+zU1hVWj5YTHMowMyWMHt7\nC8Mdllwij8r33FCylDuezvLhHz7PG8+YbW9O6q5HZcJ2TDXo5Stm8daz5/A3P3/xeOgC4J3flSGd\nez5YUrN/97oDLOpsspZUhorqk2LG9UNpmQXnfxweWVPy+UIIHt1qo/sgwP4/wdxzq8rlRk2oDzTA\nNd+Hp75W1uF4dlcvdz9/wN5aRop33jHuSxcvncGrl0zjvx4q6h90yuVw9e3ws+vkKWIMg4kMtz26\ng1uuWmHtnhk6LHMTVZyeejP5Dfih0h8COf/SRvgkm5YXYMGrKz7trLntbD40eHzEWrgVrv0p/HHN\nuKKJJ7Z3818Pbef7N5xnvboun5dHzNa5sPh14778+TevYFpTAzffs2m0oVj2Jnjvj+Dej8oy3iIe\n3nKMgN/g0mUWb4hD66GpE6Ytrvi0SIN/fC8UfwDe+CVZBVii1N+2fBBK9tkoxbLZLaNnZL7mb6Cx\nAx7+l1HPS2fz3PTjDSzubOKWN1u8QUEaifXfl+GTEnzuDcsYTmb59pO7jj/YEJGfD8OQxiJ9fHpQ\nIp3jG4/u4OYrl1l//wrqk2IaTSlhMRd8QhrdjeMTeNuODWMYcKodPfz+6idWGGPAQW6Ab/k6/OIG\nSIweGP7c7l4+8dMNfON9KznDzsZqhk9mn1Xyy/981Wncv/EIL+wvygMsvRLe+R24630y31PENx/b\nyetPm8Vyq316Nv0cTnubDGkqZHIb8FivTIr84kbp9RV4bGsXW48Os3qhBT2tyeEN0qsKV/7QtEWC\nzGoNs7O7qFH9jKXwtttkPLyQ1Hzl8BCfvftFbv+Lc1liNcYqBPz+72VhzHvuHEleFuPzGXz9vedw\naCDBV/4wpiPhwgvhA/dJTfafv1t4Sel9//VlFmPfUAiflFefmEQK08/zY8dmnXwpnPEuuP3CURrk\nZCbHQ5uP2kuSZZLy5iyTWC5m+eyW0dN5fD55Mtn6wMimlssLPvvzFwkHfHz5nWda78OSy8KvPwnR\nY2UNeNDv47brVvK/T+/hz3uKEnWBELz7B9A0A378rhEJ253P7mXVwg7OnGfRUFUxUsXIoQ5jDHhD\nE1z3U3jyK/KeKUpqmuoTWxXLpgKlCks6m8dPqD/tbfKEVxQP/9PuXj7+kw1847qV9k7PuawsXa+Q\n2O1oauCf37KCf/jlS6SL1UKnXAHv+h/4+Qdgj6zt2N8b5551B/js65dae38hCuETdeoTk8ltwJum\nw8eehdaT4PYLGPrTj/jUTzfwr7/ZzHc+sMpe57B95eWDYzm7OA5uYiY1734/XXs2ceMPn2fN207n\nvEU2NpHHvihvivfdLT23MoSDfr57/Woe2nyUHz23b/QXZ58BH/odPPtNePzLPLrlGHkBrz/NRmxz\n2+9gaXUD7i8ky5Kl5nRe/i/wnh/Cw/+K+MWHeXTDK1zxtSeY1xGxVxl7eIPcIEPV46HLZrewdeyU\n+sYOuRk++DnEuh/wr/eupz+W5r+vW2m9KjYdh7vfD9Gj8MH75amrDHPaG/nPd5/Np+96gd5oUcjE\nH5BH9pnL4c63M9TfxR1P7uZvrBoJkN73aW+vGj4B+Rkp2VJ2zkr4q4JW+faLGN78EP92/yt896nd\nvONcCzpnk0xCJqvnra761EWdEfb2xMZv9Ff+H6kjv+eDbHrhT3z8Jxu47bqV9pRjw0fhzrdBcgAu\n+nTFp77t7DnM7WjkO0/sGv2Fky+Vie97Pgi7H+c/fr+VG1+z2Hrp/JGN8vcx//zqz60zk9uAAzRE\nEK//dx5e+Q2O/O6/+Oyxf+ShGxbZG58GFfXfYzl7XhsvHigx5PjivyO15EqCd76Vn7b8N2/tOGD9\n/dfeJmN4f3Fv1VMAwLSmBr7/ofO47ZEd/PGVMVLGaYtJXv8gwy/+isR9n+FzF02z7ln17ZYx27mr\nLD29ORTg7ucPcHigxOzDBeez591/4Pf74Kz7r+K75x3hm+8/117RzL61Zav8xrJ8dut4Aw6yP8l1\nd7PnqZ/xmVfey/dXrCcsLPaQiffBj94hr8l1d5VVKBVz6fKZvP2cuXz25xtHGy2fD676GtGTLiD+\nnTdw4/zDnGJVmTRSYVg9fAIwqzXM1x7exs/XHWBwbNK7IULy9V/mgSW3ELvnJl6/76s8/MnzOGe+\nDV3+ofVyM7Lw+2gJB2kJBzg6lBz9hUAIbvgtB8LLmPfrd/P7eT/gotbqdQQj7H0a7rhEhhrfd4/c\nrCtgGAb/fvUZ/O8ze9jZNeZzsvhi8tfcSfruGzh1953cuNJGAnXjXTJ5WeLE7DT+W2+9VfmbluLx\nxx+/9ZJLLrH9fXt7YnzyZxt44oify973tyyLxAjc/wkQeekpRaaDUeYXKwT07oRXfgMv/Bje8IWK\nnu/ItyGPvysXdLCzK8rmw4Ns2DfAs3v6+K+t09m9+FquPq0V46F/kjddUydMO7m857T+h9KA3/Bb\neZqwSHukgfMWT+NTP3uB85dMRyD4zcbD/PcjO/iX3+3jxbYreGfkRVZt/iLGrkelTrV1LoQqhHRe\n/Jk0ViuusrSGpbNaeGxrN1/63Rbu3XCQA/0JfIZBazjA//3jDv7lt9tY/pqrOf/iNzDrqX+WBjnY\nJG+2SvHCfF5KB9f+t9SXz6juqbaGA/znQ9s4d0EHe3pibDs6zOYjQ2w8OMA9O/J8q28V77vmWpq2\n/VLG6HNpWTEXKONpDR6CO98uO1Ne9TVbvS1evWQad/15Pz2xNCtOauGp7T3c+ew+vvDgFr65bx5z\nO5r4QPou/E9/VXqh4VbZza/cZ+ToJhkGumKNJQ/8kmUzaA0H+cPmo6y5/xXW7esnlxfM7Wjk9y8f\n5aYfr6cnOJfz3/kpzhp8nKa1/ykNkL+h8j0D8tpsuBNaTpLJQAv8ccsxktkc/bE0O7qi7Dg2zPZj\nUdYdjPHptWFefc3nWNY4CA98VuaiZiyTIadSCCE/F3+4Ba7+ljz5WnRQWsNBGoN+vvXYLt52zhw2\nHhzk1y8e5vbHd/EvTwyxqWElnzxpCx2P/1OhXYaAjoXlFVC5jAwDvfHL8lRTJ9asWfPErbfe+ni1\n5xmihu5tTnDrrbcKu5tJXyzNlV9/gr+6+GQ+dNGi40fi3l3w5H9KiVD0mJRczTlXZsw7Fstj+b5n\npDHxBaTnvfwqy95NMpPjnd9aiwDaGgO0hoO0NgZpawxyUluYGy4srCWXhS2/hmf+n4x7zj4DWufJ\n6TZtc6Ftvpzx+Mi/wQ0PjFM2WOXhV47x1z97gXDQx+uWzuDS5TN53dIZtJuDcdNx2PWI3Kh2/AFm\nrJAhn8Zp8kYt/vPsN6QUcLm9bHouL9h0cIAntnfz+LZuXjo0yNvPnsM/vGn58UGw6biUj+15UlaZ\nNs+EeefB/POkFG1wv7x5D78oj6WNHXLG4lu+ZimEAvCP925ix7Eo4aCfcNBX+NtPcyjATa87mdlt\nhbV0bYWnvya1wDNPh/YF8k/HQvm3LwC//IiUaF70acsGopgjgwneetvTJNI5Vi7o4KJTOnnNKZ2c\nNqf1+KDg7m3SMG++Vx7DT3u7TB77G+QfX0D+vfUBqSR5vf0BvcPJDI9s6eK3mw7z1I4els9u4R/f\nvGK08mXrg7D1tzKEF+uBeatkSGD+q+RG17VFNnnr3iqbQTW2y9CUhRAKSDXWw68cwzAMfAb4DAOf\nT3rF73/VguNhk1RUtmR49psQKSTSOxYf/7t1jiwWGz4C1/xQdn60ST4veM93nuWlQ4OcPKOZVy2e\nxqsWT2P1oo7jYZN0TIYSX7pH2olTrpD2I9goeyEFwxCMSAdwy/1wo71mbtUwDGONEOLWqs+bzAYc\n5Dis1kolrol+aRAOb5BGo2+PPFIvvEj2PGlfWNPNaQshpEHq3yO9usGDctDq4EF5077re1V1vdXo\nGk4yvSlUfYJ4NiULXHY+LA2qyI/+E2qRKpIJZtMzuXzlSst8Thqvg8/LP0c3yQ1tzjlys52zsq4e\nTVmGj0qjNLD/+J/+fdJAvO7vYeX7J/Tyg4kMoYCv+iBrIaDrFbnJxrqk0cxlCn8K2vIr/12WxU+A\nZCZHg99XOXkb65FN3fY/J69NICw17TOWF/5eZinMNyEySdlnvG+PvG/Mv/v3walXyt9FhbqAaqSy\nOVLZfGXbYRLvk/mHnp2QTch7NpOQ1biZOFz4aTi1fMVyLXjGgBuGcSmwFJgvhPincs+r1YBrNBrN\niYZVA64i6v5RYBtwnWEYowLMhmF81DCMdYZhrFu/vvahtBqNRjMVUWHAixsrjDq3CSHuEEKsFkKs\nXrXKmupBo9FoNBIVBvzbwErgXiGEvY7uGo1GoymL4zN/hBAPAw87/T4ajUYz1Zj8hTwajUYzRdEG\nXKPRaCYp2oBrNBrNJEUbcI1Go5mkeKYS0zCM7wEHa/jWVcBUE5FPtZ9Z/7wnNvrnHc88IcRHqr2Q\nZwx4rRiGsU4IYa0hwwnCVPuZ9c97YqN/3trRIRSNRqOZpGgDrtFoNJOUE8GA3+H2Alxgqv3M+uc9\nsdE/b41M+hi4RqPRTFVOBA9co9FopiTagGs0Gs0kRRtwjUajmaRoA67RaDSTFMfbyTqJ1XFtJwKG\nYcxGZq+/CrwF6AJ+K4TY4urCHKBwXf8DOATs5AT+WQEMwzgf+DywA8gBceB2IcQxVxfmMIZhrAB+\nADwC7AW2CyEed3FJjmEYxqnAT4E8sAnYTh0+05PdAy87ru1EQwhxFGgFlgPnAd3A9a4uyjmeAd5Y\n+HOi/6wAfwLuAy4EZhYeu9q95ThPwaCdDCwCrkUatL9yc00O835gPtCHvIfr8pme7Aa87Li2E5gT\n/mcWQqSBy4HzAXP0+An5swIIqeV9hePGG07gn7fAJ4EbgBmAOeL+hP2ZCwOKlwJJYEXh4Qn/vJM6\nhMIUGtdmGEY70AEMIW/2ecBPXF2UQxiG8RfAvwEp4DFO4J8VwDCMLwHPA19D3txh4NeuLsphhBCf\nNgxjEbAa+BVwFvA9N9fkJIZhfAppsH8AHKFOn2ldyKPRaDSTlMkeQtFoNJopizbgGo1GM0nRBlyj\n0WgmKdqAazQazSRFG3CNpoBhGLMMw5jl9jo0GqtoA67RAIZhTAPuBU5yey0ajVW0AddoJB9HVkJ+\nwjCMk91ejEZjBW3ANRrJ04W/vymE2OXqSjQai2gDrtFoNJMUbcA1Gkne/IdhGH43F6LRWEUbcI1G\nsgnZzvWDQKPLa9FoLKF7oWg0Gs0kRXvgGo1GM0nRBlyj0WgmKdqAazQazSRFG3CNRqOZpGgDrtFo\nNJMUbcA1Go1mkqINuEaj0UxS/j8n16SAF00haAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f41cd6da6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we didn't need to do this again: if the cell above was run already,\n",
    "# the libraries are imported, but we repeat it here for convenience\n",
    "from numpy import *\n",
    "from matplotlib.pyplot import *\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "t = arange(0, 50., 1)\n",
    "\n",
    "# parameters\n",
    "r = 2.\n",
    "c = 0.5\n",
    "e = 0.1\n",
    "d = 1.\n",
    "\n",
    "# initial condition: this is an array now!\n",
    "x0 = array([1., 3.])\n",
    "\n",
    "# the function still receives only `x`, but it will be an array, not a number\n",
    "def LV(x, t, r, c, e, d):\n",
    "    # in python, arrays are numbered from 0, so the first element \n",
    "    # is x[0], the second is x[1]. The square brackets `[ ]` define a\n",
    "    # list, that is converted to an array using the function `array()`.\n",
    "    # Notice that the first entry corresponds to dV/dt and the second to dP/dt\n",
    "    return array([ r*x[0] - c * x[0] * x[1],\n",
    "                   e * c * x[0] * x[1] - d * x[1] ])\n",
    "\n",
    "# call the function that performs the integration\n",
    "# the order of the arguments is as below: the derivative function,\n",
    "# the initial condition, the points where we want the solution, and\n",
    "# a list of parameters\n",
    "x = odeint(LV, x0, t, (r, c, e, d))\n",
    "\n",
    "# Now `x` is a 2-dimension array of size 5000 x 2 (5000 time steps by 2\n",
    "# variables). We can check it like this:\n",
    "print('shape of x:', x.shape)\n",
    "\n",
    "# plot the solution\n",
    "plot(t, x)\n",
    "xlabel('t') # define label of x-axis\n",
    "ylabel('populations') # and of y-axis\n",
    "legend(['V', 'P'], loc='upper right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 3.        ],\n",
       "       [1.0050839 , 2.97163869],\n",
       "       [1.01033622, 2.94355311],\n",
       "       ...,\n",
       "       [2.23180797, 8.7761318 ],\n",
       "       [2.17956623, 8.69839475],\n",
       "       [2.12937238, 8.62112553]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interesting thing to do here is take a look at the *phase space*, that is, plot only the dependent variables, without respect to time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial condition: [1. 3.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f41cd71a7b8>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEFCAYAAAD69rxNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHjFJREFUeJzt3Xl4FEX+BvC3QiCAQMINEuUQFxVC\nwARWQEXEg0VYwAt1lQhqUBTEA3f9gTgkKCqu6AIirFcAwQvlEkUUkRVcOeQYBUFEXbnkDip3Ur8/\najozk/uY7prqfj/PM09310xS3wn6plJT3S2klCAiIvPE6C6AiIjKhwFORGQoBjgRkaEY4EREhmKA\nExEZigFORGQoBjgRkaEY4EREhmKAExEZKlZ3AZY777xTJiYm6i6DiEi7MWPGvCKlvLOk10VNgCcm\nJsLn8+kug4hIuzFjxuwozes4hUJEZCgGOBGRoRjgRESGYoATERmKAU5EZCgGOBGRoRjgRESGipp1\n4J53PBt46uzIf9+z/gykLQBi4yL/vYlIKwa4Dv53gTl3ONPXL18BYxsU/fyI7cAZdZ2phYgiigHu\nhBO/A+OaFP+aa54DOtgU6utnA3PvLvy58S0KtrXtD1w7zZ5aiChiGOB2yc0FMmoX/tzNbwGtejhX\nS7ub1aMwKycCH48Kb9v4lnqEuuszoMmF9tRHROXCAI+09bOAufcUbPdlO19LaXQeqh6hft8HPNsy\nvO3f3cKPez4LdLzL3tqIqFgM8Ej570vAR38Pb3t0BxBXU089FVGjfsFfOCsnAR+PDB4velg9LAM/\nApp2cqY+IgLAAK+43RuAqZeGt0XraLsiOt+nHpZv5gDvDgoev5ZvSmj0QSCmkjO1EXkUA7wifPH5\njl0Y3EVpc516WD78O/DVS8HjjDrB/TuWAGd1dK42Io9ggJfH1zOA+SGj0ccPA0Loqyca/OVp9QCA\n3JzwAH/lyuD+uVcDf3vb2dqIXIoBXlaho+7Ow4CrMvXVEq1iKoX/NfJ6L+Cn/6j97xeH/AwF4Dvs\neHlEbsEAL4vQ8Oaou/RuXxjc3/YJMNOaepHBn+kZ9YER2xwvjchkDPDSyD8l4KW57khreUXw5xf6\nAfAf+4Jh3igJuPsLPfURGYQXsyrJ0YPB8K7bkuEdSY2T1c/Tlw0M/DDYvsevwtwXDywdq68+oijH\nAC/O8SPAM83V/hU+YOhandW4W9POwTC/aVawffn4YJhnl+o+r0SewSmUouScAp46S+3/dRJw4W16\n6/GS864J/qXz3mBg45tqf0Lr4Gv4lxCRPSNwIUQjIcR8IURXIcREIcQOIURXO/qyTWY9tb18FMNb\np2unBkfmoaxR+bx79dRFFAVsGYFLKfcIIWoBqAJgGAAJoMCVnYQQ6QDSAaBXr152lFI+eR+mtQUu\nHaG3FgqyQvzIbuC589T+upnqAZh76QKicrJ7DvwUgLoANgD4MP+TUsppUspUKWVqSkqKzaWUUuhS\nwbv/o68OKlqtxsFRefuQv47GJap/v+l99NVG5CC7plASoEbcZwOYC2AogH/Y0VdELX82uM85VjP0\nmaT+rUYfCrZtXxacYsk5pa00IrvZNYVyGEBy4HC6HX1EnJTA0sBZlQxv88TEBP/dFj4ArHlV7Vuf\nZVw0BOgxTk9tRDbhMkLLmAS1vSFLbx1Ucb0mqDAftS/Y9t8Xg6NyIpdggAPAnDuD+6376quDIiu2\nSnCu/Mz2wXYryA9u11cbUQRwHTgA+N9RW06duFf6MrU9uB34VyDMrW2HO4Fr/qmjKqIK4Qjc+pO6\na/R/xkoRUKdFwXXlq1/m9AoZydsBnpsb3O/2qL46SA8ryGslhrRx9QqZw9tTKNZd44d/o7cO0uvB\nb9V2y0fA7P5q31q9ctdSoEmUnKNAlI93R+BHDwb3E87SVwdFj1Y91Ih85K/Btn9frkbky8frq4uo\nCN4NcOsqg6EngBABQOWqBefJl45VQT6tm766iPLxZoDnnA7ux3jzR0CllD/Id33NDzwpangzvTLr\nqu2jO/XWQeawgrxxu5A2Bjnp5c0At8TV0F0BmWbw5yrIu40KtllBLqW+usiTvBfg1ohp0GK9dZDZ\nuo5QQX7b3GDbmAQGOTnKewFuOfsi3RWQG5zTTQX5fSG322OQk0O8FeDbPtVdAblVvcANrx/cHGxj\nkJPNvBXgM69VW17zhOxS60z139dDW4JtVpCHnvlLFAHeCnAip9RspIL84e+DbRm1GeQUUd4J8HcH\nqe1VT+itg7ylRoOig5xTK1RB3gnwb+aobef79NZB3mQF+UNbg23W1ApROXknwImiQc2GKsiHrQu2\n8YQgKidvBPjHgZMu+k3TWweRxbou+e0fBNsY5FRG3gjwlRPVNrm/3jqI8mt2sQryq0NuuOyLB15o\nV/TXEAV4I8CJol2nISrIW12jjg/9qIJ80Qi9dVFUc3+AH/4lsCO0lkFUKjfPCj9PYdU0FeRbeekH\nKsj9Af58G7X1HdZbB1FZ5L+M7awbVZD/cUBfTRR13B/gRCbLH+TjW3ANOeVhgBOZwJcNPLY/eMw1\n5AS3B/i2T9T2yky9dRBFQqXKKsjv+CTY5osHsnrrq4m0siXAhRCNhBDzhRBdhRDjhRAjhBDn29FX\nsWZep7ZdhjneNZFtzuqggrzN9er4x+UqyH9YqrcucpwtAS6l3AOgFoDzAHQAsA/AgPyvE0KkCyHW\nCCHWrF27Nv/TRFSc618Jnx+f0U8F+ekT+moiR9k9hVIlZL/AOj4p5TQpZaqUMjUlJcXmUohcKv8H\nnWMbcH7cI+yaQkkAUBvAEQCbACQCeMOOvop04je1bXGZo90SaePLBh75MeQ4HpjMO0+5mV1TKIel\nlMlSyiwp5RAp5Vgppd+Ovoo0LlFtB8xztFsirarXUUHefbQ63rdZBfnezcV/HRnJ3atQiLzqkofC\np1VevIjTKi7EACdyM1828HjIWci84qGrMMCJ3E4IFeQ3zQq2+eKB7cu0lUSR4c4AP/CD2vZ4Wm8d\nRNHkvGvCp1Wm9+Fp+YZzZ4BPvFBtL7pbbx1E0Sj/skOelm8sdwY4EZWssDsC/bpJXz1UZgxwIi+z\n7ghkmdKJo3GDMMCJSIX46EMhx/HAksf11UOlwgAnIiUmRgV5x8HqeMXz/JAzyrkvwL95T22HrdNb\nB5Gpej7DDzkN4b4Af3eg2tZpobcOItP5soH0ZSHH8UD2Dl3VUCHcF+BEFDlntg8fjU9ozdF4FGGA\nE1HJfNnAyF9DjuOBn7/UVw8BYIATUWlVrqqCPP5sdfxaD47GNWOAE1HZPOAveIGsn1boq8fD3Bng\nMZV1V0DkbtYFshq2Ucev9+RoXAN3BfiaV9X27z9pLYPIM+5ZwdG4Ru4K8IUPqG1cDb11EHmJNRpv\n0FodczTuGHcFOBHpM2RlwdH4nm/01eMBDHAiihxrNF4rcE/al7pwNG4jBjgRRd6D3wKP7Q8e++KB\nU8f01eNSDHAiskelyuFncT7RCJjSRV89LsQAJyJ7+bKBe1ep/V+/4ZRKBDHAich+9VuFj8Z98cCm\nefrqcQn3BPgfB9Q2/XO9dRBR0XzZQOdhav/tARyNV5B7AnxGH7U9s53eOoioeFdlFlxuyA84y8U9\nAb7Hr7sCIiota7mh5YlGwKz++uoxlO0BLoR4WAiRJoSYJ4QQdvdHRAbxZQP3rFT7Wz/ilEoZOTEC\njwXQDcB2KcNvrieESBdCrBFCrFm7dq0DpRBR1GnYuuAHnJxSKRUnArwygLkArhdCVAt9Qko5TUqZ\nKqVMTUlJcaAUIopa+adUFo3QV4shnAjwWgDiAHwC4JQD/RGRqXzZwF2fqf1V0zilUgLbA1xKOUJK\n+ZaUcqCU8rTd/RGR4ZpcWHBKJTdXXz1RzD2rUIjIXUJDPKM2sO4NfbVEKQY4EUUvXzZw/Wtqf94Q\nTqnk464ATx2kuwIiirQ21xY88YcAuCXAD/ygtldm6q2DiOyR/8QfXzxw8qi+eqKEOwJ8yWi15a3U\niNzNlw0kNFX7TzYG1mbprUczdwT4dwt1V0BEThm+Ebh9kdpfMMzTUyruCHAi8pZmXYDRh4LHHg1x\nBjgRmSkmpuC8uMcwwInIbPlD/PRJfbU4jAFORObzZQMNk9T+2PrADm9cHI8BTkTucM8XwPWvqv2X\nLwc+fkxvPQ5ggBORe7S5Dnh4m9pf+S8gs4HeemzGACcid6lRP3jmZs4JV3+4yQAnIvcp7MxNFyoy\nwIUQLYUQbwkhnhNCVHGyKCKiiHB5iBc3An8ZQHcAgwEMd6YcIqIIc3GIFxfgG6SU9QA0AJBgNQoh\nmtheFRFRJLk0xIsL8A5CiFcATARwqRDiVSHEqwDecaa0Mmpzne4KiCiauTDEiwvwNgAuh7qjfJPA\nthuA1g7UVXYpA3VXQETRzmUhXlyAd5dSNs//AHCZQ7WVzp5v1LbZxXrrICIzuCjEiwxwKeXqItrX\n2VdOOawN3G5JCL11EJE5XBLi5q8DX/Oa7gqIyEQuCHHzA1zm6K6AiExleIibH+BERBVhcIgzwImI\nDA1xBjgREWBkiDPADSGlxOYDmzFlwxQM/XQopm6YiiMnj+gui8hdQkPcgEvRxtrdgRAiFsA9AL4C\nsF5K6Z37HZXDjQtuxOaDm0t83bIdyzBp/aTiX3PjMtStVjdSpRF5gy9bjcBzTgArJwKdh+quqEi2\nBziAkQDqBR5rQp8QQqQDSAeAXr16OVBK9JFSou30tiW+LrNLJvqc0wciZL373qN78cLXL2D+D/ML\n/ZrL3r4s7HjqFVPRuUnnCtVL5AmPHwbGJAAfjwKSbgRqNtRdUaGElNLeDoSYD3U9lScBDJZSfl3Y\n63w+n/T5fGXvwJqrCv3TxwCzNs/CuFXjCrTP7DkTyfWTI9LH3UvuxopdKwp9LqVhCl7v8XpE+iFy\npX1bgMkd1b7D+SKEGCOl9JX0OidG4FMAdIIafW9yoL+oNmfrHPi+9IW1fXHTF4iPi/yHJi9d+VLY\n8aDFg7B6jzrBdu2va5GUpW4CO77rePRo1iPi/RMZrX4roO1NwMY31UAxCgeJto/AS8vtI/DCpkr8\naX5N1QBTN0wtdA5dZ01EUSl0RYpDOVPaEThXoThg6NKhYeHtT/NrD8rByYMLrSMpKwlJWUmIll/s\nRNqFhvaSx/XVUQgGuM2SspKw7JdlAICMzhnag7swVpBPuGxCXlvb6W2RlJWEXJmrsTKiKGGF+Irn\ngdzouXwHA9xG1hwzoEKy37n9NFZTsiuaXgF/mh9Lb1ia15Y8PTnsfRB51oDAaq+MOnrrCMEAt0n+\n8DZJ/er14U/z46PrPsprs6ZWiDyrRdfgfpScqckAt4HJ4R2qSY0m8Kf5kdklM68tKSsJmw54fjER\neVXofPihn/XVEcAAjzC3hHeovi37hr2X/gv7czRO3jVyj9q+UPIJeHZjgEfQVe9elbfvlvAO5U/z\nY/1t6/OOk7KSsPfoXo0VEWlQuVpwX/NUCgM8Qk7nnsbuP3YDcGd4WyrFVAp7f93f6c7ROHlP6FSK\nxlUpDPAIaT+jPQB45vR0f5ofS65fknfMECfPGRq4KojGVSkM8Ai4YcENefspDVM0VuKsRmc0ChuN\nM8TJU+qeE9z/4TMtJTDAI+C7g98BcPfUSXHyhzhP/iHPsKZSZvTV0j0DvIKsUeefav9JcyV6+dP8\niBXq2mjJ05Nx7PQxzRUROew/zznepTsCvHYz3RVgzl/n6C5Bu3UD1uVd1bDjGx2RE0WnHBPZxhqF\nfzrG8a4Z4BVgjb5bJrTU0n80Gt91PK44+woAQLsZ7TRXQ+Sw7z5wtDuXBHhzrd2/3+d9rf1Hmwnd\nJqB2XG0A/GCTPMIahb95i6PduiPA67RwvMvjp4873qdJlt+0PG+fIU6e4uClmBng5dThjQ4AgA0D\nNjjetylCV6f4Vvr0FULkhEd3qu2YBMe6dEmA65tCiRHu+BHaxQrxOd/zQ15yubgajndpdvrknFbb\nhLP11kHFuqPNHQA4lUIecvyII92YHeB/7FPbuJqOdtv97e4AgAV9Fzjar6mGpwzP2//t5G8aKyGy\nmfVh5lNnOdKd2QGevUNLt3uPqSvwNYtvpqV/E311y1cAgM6zO2uuhMg9zA7w33/VXQGVUvXK1XWX\nQOQ6Zgf4yd91V0BlMLPnTACcCyeXq3+e2jpwJrLZAc5TtY2SXD9ZdwlE9hvyX7V91v7rI5kd4Bqv\nelenavTcmZqIoogQant0v+1dxdreg4tkL1iA3f98Fm/uOY3KjWOQXX0B4nv31l0WEXmU2SNw6zed\nA7IXLMDux0ZD7tmLGAA5u/dg92Ojkb2ASwnLYnl/dYr9Z//TcwF8IjdxJMCFEOcLIb5yoi+77J3w\nPOTx8OufyOPHsXfC85oqMlPtquoiVy/7X9ZcCZH5bA9wIcS5AM4B0LCQ59KFEGuEEGvWrl1rdykV\ncnr37jK1U/FiYzh7R1RRTozA7wNwO4B6Qoiw9WNSymlSylQpZWpKSjnuJengh5ixjRuXqZ0KN2er\nuibKpO6TNFdCZD7bA1xKeT+AhwHsl1Iae9PIBg8Mh6haNaxNVK2KBg8ML+IrqDC+L30AgJpVnL38\nAZEbOfJ3rJTyJwDNIv6NHbwSoLXaZO+E53Fy1y4cqAUkP5bBVShEpI3ZE5EOz6PG9+6N+N69884k\n9DO8iSi/LR+qbZr9K9TMXkYYG6e7AioD6xffkuuXaK6EyEazb1Lb5pfa3pXZAV6jwMIWilKnc0/n\n7Tc6o5HGSojcw+wArxe41kCus6fU16pSCwCwcudKR/s1WfsZ7QEAC/st1FwJkY0cvj6T2QFeva7a\nOnxZ2RU3rwAADP5ksKP9mir06oNNazXVWAmRzTIC10gatc+R7swOcOtU+n2b9dZBRUqdmZq3H3qT\nYyJXi63iSDdmB7hl73eOd3nbBbcB4LWti5OUlYQTOScAMLzJA3zxatvzWce6dEmAb3K8y0c6POJ4\nnyYJ/cW2ccBGjZUQOUDK4H7Huxzr1h0Bvm+Llm57t1DrwDkKDxf68/Cn+SEcvGokkRZjEtQ2fZmj\n3Rod4HPX7QQAHPnFjy5PLc07dsqTlzyZt//O1ncc7TsajV89vkB4E7ne1K7B/TPbO9q1sQE+d91O\nPPqeCoha4hh2Hj6GR9/zOx7i1vRAxpcZefO9XpMrc5GUlYTpm6YDAHq16MXwJm/I3gnsXq/2fdmO\nd29sgI9fvAXHToWvuTx2KgfjFzs7nSKEwEMpDwFQKy5k6FyYByRlJSF5evBel/40P8ZdMk5jRUQO\nmnCB2j7g/OdwgMEBvuvwsTK12+n2NrejeXxzAEDb6W2x/5j998LTSUqJpKyksOmSBX0XcNRN3mKt\nOkm9A4hvoqUEYwP8zIRqZWq32/y+89E1Uc2FdXu7G3q/774LXS3fsRxJWUloO71tXtvUK6bCn+ZH\ns/hm+gojcpoV3gDQ6zltZRh7NcIRV7fKmwO3VKtcCSOubqWpInWTgjV71mDg4oH46chPSMpKwoYB\nGxDj4GVvI+1UzilcOPPCAu2f9/8cdarW0VARkUZSBlecAFrmvUMZG+B92wf+ZJmnNk0SqmHE1a2C\n7ZqkNkqFP82fN71gzQ+bNL3wy2+/oOd7PQt9zqT3QRRRR3YBz52v9kUl4PGDeuuBwQEOBEI8EOAr\n/nG53mLy8af5sf3wdvSZ1wdAcG30i91fxCWJl+gsrQApZdi0SH7rblvHe1iSt71/N7Bhttr/8z3A\nX57SW08A/6+0UYuEFvCn+bHt0Db0m98PADDk0yF5z1+aeCkmd5/seF0jvxiJ+T/ML/Y1HGkTBYTO\ndw9bB9Rpoa+WfBjgDmhZu2VeIIau3LA+FAzVtl5bzOg5o8Lz5l//+jXSPkor9evX37YelWIqVahP\nIlf5ZRXwypXBY83z3YVhgDssdGS7es9qDFo8KOz5jfs3hq2rtsPH132MxjUa29oHkdFCR90JTYHh\n0Xk9Hwa4Rh0adSgwVSGlxOiVozF329xyf99qsdWw6NpFqFetXkVLJPKW9wYDG98MHo8+BMRE7yoy\nBniUEUIgs0smMrtk6i6FyDt2rQemhVzTpP1tQJ9J+uopJQY4EXnXsUPA083C26Jwrrso5gd4g9bA\n3m91V0FEJgld0215/HDwLl+GMD/Ab31X/UPs3wbUa6m7GiKKZt9/ArxxXXjbYweASmZGoZlVh6p1\nptpOSjHqTx8ictDTzdR0SSgDR9z5mR/gRESFOfEbMC6xYLuLBnq2B7gQohuApwHslFL2s7s/IvK4\n0DXcliszgC73O1+LzZwYga8A0APA0vxPCCHSAaQDQK9evRwohYhcKf/6bYsLpkmKY3uASylPCiG6\nA7hWCBEnpTwR8tw0ANMAwOfzeetWNkRUMXPuBPyF3Iv2rqVAkxTn69HAiSmUWwFkADgBoFNgaw8p\nXf3blsjzCpseAYCrnwQ63etsLVHAiRH4TAAzbe3k8cPqIutjElz1AQWR5+3bCkzuUPhz170CJF3v\nbD1Rxh2rUDjqJnKH3Bwgo5g7PT20BajZyLl6opw7ApyIzFRSYDdJBe761Ll6DOO+AD91DKis58bG\nRFSC7J3AhAuKf02UXwEwmrgnwK158CcacR6cKFq80A449GPxr3H5Uj87uSfA+R8AkV7v3A58+37x\nr7lwAPDXiY6U4wXuCXAAaNoF+HkF8NwFwIObdFdD5E6njqm/dEtj1F4gNs7eejzMXQE+cJFaJ3pk\np+5KiMwnpbon5I7VpXv9ozuBuBr21kRh3BXgAND0YuDnL1SQcy6cqGRHDwLPNC/968/vDfS399QO\nKh33BfjAD4Jna335ItBpiN56iKLByaPA+JbAqT9K/zX1zwfuWckVIVHMfQEOBFekLH4UaNEVaNha\nd0VE9vv2ffVBYln1mQy0vzXi5ZD93BngQgD3rQEmpQJTOgPX/BPocKfuqojKLzcXWPx/wFdTyvf1\nV2YCXYZFtibSzp0BDgD1zgWGrQf+1Q744CH14Jw4RaNjh4HXelbs3q41GwP3rgKq1opcXRT13Bvg\nAFCnuQpta07cFw80TgYGL9dbF3lDYfdfLK/L/g/o+gjPd6Aw7g5wiy8beP9uYMNsYPeGYKDzlF0q\ni8P/UzcO+N/KyH7fQYuBsy+K7PckT/BGgANAv5fUY/YtwJYPVFtG7eDzaQuB5pfoqY2clXMKWP0y\n8NE/7OujYzpw1ViexEK28k6AW26epbY5p4HMusH2rEJu6Xb/BqB2M0fKolKSEti/FVg6Ftg837l+\nk24A/vIMUL2YK+cROcx7AW6pFBv+oebn44HPxoa/5oXk4r9Hx8HqTiCVvPtjLNbJP4Cti4FN84BN\nc3VXozTtom5wm5iquxKiCmPyWLqOUI9Q3y0C3ry56K9ZNVU9nBBXK/CoAUAAJ44Ax7OBk7870380\nqHsucPFwIOlGILaK7mqItGOAF+e8niUvPSwp5CPlxBH1iHaV4tTP7fzewLlXAXE1dVdE5FoM8Ioq\nTcgTEdmAa+iIiAzFACciMhQDnIjIUAxwIiJDMcCJiAzFACciMhQDnIjIULavAxdCdAPwJwBnSSlH\n2d0fEZFXCCmlvR0IMRvAVACvAEiSUh4NeS4dQHrg8BiAT0v5bVMArI1knYbh++f75/t3t0QpZYm3\nEXMiwOcAmAgV4G2llGW4q2qR33ONlNKzVyPi++f75/v37vsP5cQc+EsA2gN4LxLhTUREiu1z4FLK\nJQCW2N0PEZHXmLoKZZruAjTj+/c2vn8C4MAcOBER2cPUETgRkecxwImIDMUAJyIyFAOciMhQxt1S\nzaun5gfe99MAdgLYBmAvgIVSys1aC3OQEOJ8AK9DnbH7E4CtUsplGktyjBAiFsA9AFYBuAnALnjo\n318I8TCAfQBuAbARwB8Apkgpf9VamGYmjsDTAWwBcLMQorruYhy0AkCPwKMD1H/MA7RW5CAhxLkA\nzgHQDCrAtgIYrLMmh42EGrgMhzqV3FP//lCDzW4AGgCoH2jrq6+c6GBigFcJ2RfaqnCYlPIkgO4A\nLgIQF2j2zPsHcB+A26H+540PtHnp/acAmA/1C7xdoM1L778ygLlQv8CrBdq89P4LZdwUCjx6ar4Q\n4lYAGQBOAPgMQCKAN7QW5SAp5f1CiGYAUgG8D6AtgJd11uSwKQA6AXgT6he4p/79AdSCet9zoaZP\nqgKYp7WiKMATeYiIDGXiFAoREYEBTkRkLAY4EZGhGOBERIZigJOnCSEeEUJIIUSWEKKmECJGCPGQ\nEOJ1IURcyd+BSB+uQiFPE0LEAPgewH4p5Z8DbX8DsEBKeURrcUQl4AicPE1KmQu1xrqjECJFCCEA\nVGF4kwk4AifPE0LUgbrGzKzAY5uU8me9VRGVjCNw8jwp5UGoMxxvAnABw5tMwQAnUiYDqA7ga92F\nEJUWp1CIAoQQk6WU9+qug6i0GOBERIbiFAoRkaEY4EREhmKAExEZigFORGQoBjgRkaEY4EREhmKA\nExEZ6v8BIJoZhYT1JdgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f41cd837048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# `x[0,0]` is the first value (1st line, 1st column), `x[0,1]` is the value of \n",
    "# the 1st line, 2nd column, which corresponds to the value of P at the initial\n",
    "# time. We plot just this point first to know where we started:\n",
    "plot(x[0,0], x[0,1], 'o')\n",
    "print('Initial condition:', x[0])\n",
    "\n",
    "# `x[0]` or (equivalently) x[0,:] is the first line, and `x[:,0]` is the first\n",
    "# column. Notice the colon `:` stands for all the values of that axis. We are\n",
    "# going to plot the second column (P) against the first (V):\n",
    "plot(x[:,0], x[:,1])\n",
    "xlabel('V')\n",
    "ylabel('P')\n",
    "\n",
    "# Let's calculate and plot another solution with a different initial condition\n",
    "x2 = odeint(LV, [10., 4.], t, (r, c, e, d))\n",
    "plot(x2[:,0], x2[:,1])\n",
    "plot(x2[0,0], x2[0,1], 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Congratulations**: you are now ready to integrate any system of differential equations! (We hope generalizing the above to more than 2 equations won't be very challenging)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For more info:\n",
    "\n",
    "* [Python tutorial](http://docs.python.org/3/tutorial/index.html) (chapters 3 to 5 are specially useful).\n",
    "* [An introduction to Numpy](http://nbviewer.ipython.org/github/iguananaut/notebooks/blob/master/numpy.ipynb)\n",
    "* [Another one](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-2-Numpy.ipynb), covering a little bit more ground.\n",
    "* [The matplotlib gallery](http://matplotlib.org/gallery.html): all kinds of plots, with sample code to use."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}