rabbits-and-foxes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rabbits and foxes\n",
    "\n",
    "There are initially 400 rabbits and 200 foxes on a farm (but it could be two cell types in a 96 well plate or something, if you prefer bio-engineering analogies). Plot the concentration of foxes and rabbits as a function of time for a period of up to 600 days. The predator-prey relationships are given by the following set of coupled ordinary differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dR}{dt} &= k_1 R - k_2 R F \\tag{1}\\\\\n",
    "\\frac{dF}{dt} &= k_3 R F - k_4 F \\tag{2}\\\\\n",
    "\\end{align}\n",
    "\n",
    "* Constant for growth of rabbits $k_1 = 0.015$ day<sup>-1</sup>\n",
    "* Constant for death of rabbits being eaten by foxes $k_2 = 0.00004$ day<sup>-1</sup> foxes<sup>-1</sup>\n",
    "* Constant for growth of foxes after eating rabbits $k_3 = 0.0004$ day<sup>-1</sup> rabbits<sup>-1</sup>\n",
    "* Constant for death of foxes $k_4 = 0.04$ day<sup>-1</sup>\n",
    "\n",
    "Also plot the number of foxes versus the number of rabbits.\n",
    "\n",
    "Then try also with \n",
    "* $k_3 = 0.00004$ day<sup>-1</sup> rabbits<sup>-1</sup>\n",
    "* $t_{final} = 800$ days\n",
    "\n",
    "*This problem is based on one from Chapter 1 of H. Scott Fogler's textbook \"Essentials of Chemical Reaction Engineering\".*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Solving ODEs\n",
    "\n",
    "*Much of the following content reused under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth. Partly based on David Ketcheson's pendulum lesson, also under CC-BY. https://github.com/numerical-mooc/numerical-mooc*\n",
    "\n",
    "Let's step back for a moment. Suppose we have a first-order ODE $u'=f(u)$. You know that if we were to integrate this, there would be an arbitrary constant of integration. To find its value, we do need to know one point on the curve $(t, u)$. When the derivative in the ODE is with respect to time, we call that point the _initial value_ and write something like this:\n",
    "\n",
    "$$u(t=0)=u_0$$\n",
    "\n",
    "In the case of a second-order ODE, we can write it as a system of first-order ODEs, and we would need an initial value for each equation: two conditions are needed to determine our constants of integration. **Check that you understand how to do this!** The same applies for higher-order ODEs: if it is of order $n$, we can write it as $n$ first-order equations, and we need $n$ known values. If we have that data, we call the problem an _initial value problem_.\n",
    "\n",
    "Remember the definition of a derivative? The derivative represents the slope of the tangent at a point of the curve $u=u(t)$, and the definition of the derivative $u'$ for a function is:\n",
    "\n",
    "$$u'(t) = \\lim_{\\Delta t\\rightarrow 0} \\frac{u(t+\\Delta t)-u(t)}{\\Delta t}$$\n",
    "\n",
    "If the step $\\Delta t$ is already very small, we can _approximate_ the derivative by dropping the limit. We can write:\n",
    "\n",
    "$$\\begin{equation}\n",
    "u(t+\\Delta t) \\approx u(t) + u'(t) \\Delta t\n",
    "\\end{equation}$$\n",
    "\n",
    "With this equation, and because we know $u'(t)=f(u)$, if we have an initial value, we can step by $\\Delta t$ and find the value of $u(t+\\Delta t)$, then we can take this value, and find $u(t+2\\Delta t)$, and so on: we say that we _step in time_, numerically finding the solution $u(t)$ for a range of values: $t_1, t_2, t_3 \\cdots$, each separated by $\\Delta t$. The numerical solution of the ODE is simply the table of values $t_i, u_i$ that results from this process.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler's method\n",
    "*Also known as \"Simple Euler\" or sometimes \"Simple Error\".*\n",
    "\n",
    "The approximate solution at time $t_n$ is $u_n$, and the numerical solution of the differential equation consists of computing a sequence of approximate solutions by the following formula, based on this equation:\n",
    "\n",
    "$$u_{n+1} = u_n + \\Delta t \\,f(u_n).$$\n",
    "\n",
    "This formula is called **Euler's method**.\n",
    "\n",
    "For the equations of the rabbits and foxes, Euler's method gives the following algorithm that we need to implement in code:\n",
    "\n",
    "\\begin{align}\n",
    "R_{n+1} & = R_n + \\Delta t \\left(k_1 R_n - k_2 R_n F_n \\right) \\\\\n",
    "F_{n+1} & = F_n + \\Delta t \\left( k_3 R_n F_n - k_4 F_n \\right).\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is one way to do it. Not necessarily the best, but it'll work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "step_size = 1\n",
    "k1 = 0.015\n",
    "k2 = 0.00004\n",
    "k3 = 0.0004\n",
    "k4 = 0.04\n",
    "end_time = 600.\n",
    "step_size = 1\n",
    "times = np.arange(0, end_time, step_size)\n",
    "foxes = np.zeros_like(times)\n",
    "rabbits = np.zeros_like(foxes)\n",
    "rabbits[0] = 400.\n",
    "foxes[0] = 200.\n",
    "for n in range(len(times)-1):\n",
    "    delta_t = times[n+1] - times[n]\n",
    "    rabbits[n+1] = rabbits[n] + delta_t * (k1 * rabbits[n] \n",
    "                                            - k2 * rabbits[n] * foxes[n])\n",
    "    foxes[n+1] = foxes[n] + delta_t * (k3 * rabbits[n] * foxes[n]\n",
    "                                      - k4 * foxes[n])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9+P/Xe7KTfQNDAgQUQUA2A2JR3BWXVmt7rV5v\n1dbq9dZW7XKr9nvvte3v67eb1lZb9Wpri13cW7UWtaKo1QoaFNmRAAGCgYQl+z7z+f3xOZMMIcsk\nmZkzk3k/H488zpnPnHPmcyCT9/nsYoxBKaVUfPK4nQGllFLu0SCglFJxTIOAUkrFMQ0CSikVxzQI\nKKVUHNMgoJRScUyDgFJKxTENAkopFcc0CCilVBxLHOwAEUkF3gJSnOOfMcbcKSLfA64Hap1Dv2uM\nWe6ccwdwHeAFbjbGvOKknwT8DkgDlgO3mEGGLBcUFJjS0tIh35hSSsWzNWvWHDDGFA523KBBAGgH\nzjLGNIlIEvC2iLzkvHevMebuwINFZAZwBTATGA+sEJHjjTFe4EFs4FiNDQJLgZcYQGlpKeXl5UFk\nUymllJ+I7ArmuEGrg4zV5LxMcn4Genq/BHjCGNNujNkJVAALRaQIyDLGrHKe/h8DLg0mk0oppcIj\nqDYBEUkQkbVADfCqMWa189bXRWSdiDwqIrlOWjGwJ+D0Kiet2Nnvnd7X590gIuUiUl5bW9vXIUop\npUIgqCBgjPEaY+YCJdin+lnYqp0pwFygGrgnVJkyxjxsjCkzxpQVFg5apaWUUmqYgmkT6GaMqROR\nlcDSwLYAEXkEeNF5uReYEHBaiZO219nvna6UUoPq7OykqqqKtrY2t7MSVVJTUykpKSEpKWlY5wfT\nO6gQ6HQCQBpwLvBjESkyxlQ7h30W2ODsvwD8SUR+hm0Yngq8Z4zxikiDiCzCNgxfDdw/rFwrpeJO\nVVUVmZmZlJaWIiJuZycqGGM4ePAgVVVVTJ48eVjXCKYkUAQsE5EEbPXRU8aYF0Xk9yIyF9tIXAn8\nu5OpjSLyFLAJ6AJucnoGAXyVni6iLzFIzyCllPJra2vTANCLiJCfn89I2k4HDQLGmHXAvD7SvzjA\nOXcBd/WRXg7MGmIelVIKQANAH0b6b6IjhpVS4VdVDh+/4nYuVB80CCilwqu+Cn59Nvzpcti2wu3c\nRERpaSkHDhw4Kv3aa6/lmWeeOSq9vLycm2++GYA33niDf/7zn2HPo9+QegcppdSQbX6xZ7/8UZh6\njnt5CSFjDMYYPJ6RP0uXlZVRVlYG2CCQkZHBpz71qRFfNxhaElBKhdf21yF/Kiy4Hna8Ad5Ot3M0\nbJWVlUybNo2rr76aWbNmcd1111FWVsbMmTO58847jzj2Jz/5CSeeeCILFy6koqKiO33FihWUlZVx\n/PHH8+KLNkC+8cYbXHzxxVRWVvLQQw9x7733MnfuXP7xj3/w9NNPM2vWLObMmcOSJUtCfk9aElBK\nhdf+DVB6Kkz6FLz/CFR/BCVlI7rk9/+6kU2fNIQog9aM8Vnc+emZgx63bds2li1bxqJFizh06BB5\neXl4vV7OPvts1q1bx+zZswHIzs5m/fr1PPbYY9x6663df/ArKyt577332L59O2eeeeYRAaK0tJQb\nb7yRjIwMvv3tbwNw4okn8sorr1BcXExdXV1I7xm0JKCUCqfWOmjYC2NPgAkLbdonH7qbpxGaNGkS\nixYtAuCpp55i/vz5zJs3j40bN7Jp06bu46688sru7bvvvtudfvnll+PxeJg6dSpTpkxhy5YtA37e\n4sWLufbaa3nkkUfwer0DHjscWhJQSoXPgY/ttvAEyCqGlCyoHfiPXjCCeWIPl/T0dAB27tzJ3Xff\nzfvvv09ubi7XXnvtEaOZA7tu9rff1+veHnroIVavXs3f/vY3TjrpJNasWUN+fn4obgXQkoBSKpzq\ndttt7iQQgcLpUDPyIBANGhoaSE9PJzs7m/379/PSS0eOfX3yySe7t6ecckp3+tNPP43P52P79u3s\n2LGDadOmHXFeZmYmjY2N3a+3b9/OySefzA9+8AMKCwvZs2cPoaQlAaVU+PiDQLYzndjY6bDlb+7l\nJ4TmzJnDvHnzmD59OhMmTGDx4sVHvH/48GFmz55NSkoKjz/+eHf6xIkTWbhwIQ0NDTz00EOkpqYe\ncd6nP/1pPv/5z/P8889z//33c++997Jt2zaMMZx99tnMmTMnpPchgyzs5bqysjKji8ooFaNe/AZs\n/AvcVmlfv/sAvHIHfLsCMoY2Q/DmzZs54YQTQp/HUaCvfxsRWWOMGbQFXquDlFLhU1/VUwoAWxIA\nqN3sTn7UUTQIKKXCp6kGMsb1vC5w6r8PbHMnP+ooGgSUUuHTchDSC3peZxZBQjLUBbX8rYoADQJK\nqfAwBpprjwwCHo+tHjpc6Vq21JE0CCilwqOjGbraYEzBkem5pXBYSwLRQoOAUio8WpxZNNN7B4FJ\nWh0URTQIKKXCo9kfBHp1Bc2ZBK2Hoa0+8nkaofvuu48TTjiBq666yu2shIwOFlNKhYc/CBxVHTTJ\nbg/vgqLZkc3TCD3wwAOsWLGCkpISt7MSMloSUEqFR7/VQaV2G2NVQjfeeCM7duzgggsu4J577uHS\nSy9l9uzZLFq0iHXr1gFwyy238IMf/ACAV155hSVLluDz+aitreVzn/scCxYsYMGCBbzzzjsAvPnm\nm8ydO5e5c+cyb968I6aLiBQtCSilwqPZWfy8dxDInmi39VXDv/ZLt8O+9cM/vy/HnAgX/Kjftx96\n6CFefvllVq5cyfe//33mzZvHc889x+uvv87VV1/N2rVr+eEPf8iCBQs47bTTuPnmm1m+fDkej4db\nbrmFb3zjG5x66qns3r2b888/n82bN3P33Xfzq1/9isWLF9PU1HTUFBKRMGgQEJFU4C0gxTn+GWPM\nnSKSBzwJlAKVwOXGmMPOOXcA1wFe4GZjzCtO+knA74A0YDlwi4n2eSuUUsPTfAAS0yA5/cj0MXmQ\nmDqyIOCyt99+m2effRaAs846i4MHD9LQ0EBWVhaPPPIIS5Ys4d577+XYY48F7EIygdNMNzQ00NTU\nxOLFi/nmN7/JVVddxWWXXeZKNVMwJYF24CxjTJOIJAFvi8hLwGXAa8aYH4nI7cDtwG0iMgO4ApgJ\njAdWiMjxxhgv8CBwPbAaGwSWAi8d/ZFKqZjXfODoRmGws4lmjYeGT4Z/7QGe2N22fv168vPz+eST\nnvvz+XysWrXqqCf922+/nYsuuojly5ezePFiXnnlFaZPnx7R/A7aJmCsJudlkvNjgEuAZU76MuBS\nZ/8S4AljTLsxZidQASwUkSIgyxizynn6fyzgHKXUaNNyANL7mfc+q3hkQcBlp512Gn/84x8BuzRk\nQUEBWVlZ7Nq1i3vuuYcPP/yQl156idWrVwNw3nnncf/993efv3btWsBOE33iiSdy2223sWDBgkEX\nmAmHoBqGRSRBRNYCNcCrxpjVwDhjTLVzyD7AP0FIMRA44XWVk1bs7PdO7+vzbhCRchEpr62tDfpm\nlFJRpLUO0nL7fi9rvF1xLEZ973vfY82aNcyePZvbb7+dZcuWYYzhuuuu4+6772b8+PH85je/4Stf\n+QptbW3cd999lJeXM3v2bGbMmMFDDz0EwM9//nNmzZrF7NmzSUpK4oILLoj4vQTVMOxU5cwVkRzg\nLyIyq9f7RkRCVrdvjHkYeBjsVNKhuq5SKoLaGyC7nzrurGJorAafFzwJkc3XCFRWVnbvP/fcc0e9\nv2LFiu79k046ifXrbeN1ampq9yIzgQJLB24ZUhdRY0wdsBJbl7/fqeLB2dY4h+0FAuaOpcRJ2+vs\n905XSo1GbfWQmt33e1njwdfV04NIuWbQICAihU4JABFJA84FtgAvANc4h10DPO/svwBcISIpIjIZ\nmAq851QdNYjIIrGLal4dcI5SarRpa4DUrL7f85cQ6vU50G3BlASKgJUisg54H9sm8CLwI+BcEdkG\nnOO8xhizEXgK2AS8DNzkVCcBfBX4NbaxeDvR3DNo/yb42Ux49ivg87mdG6ViS1cHdLUOXBKAIbcL\naI/yo43032TQNgFjzDpgXh/pB4Gz+znnLuCuPtLLgVlHnxGFVtwJDVWw/mmYdiHMusztHCkVO9ob\n7DalvyDg9AkZQhBITU3l4MGD5OfnYysTlDGGgwcPjmiQmY4Y7kvLIahYAYtvhQ1/hrV/0iCg1FD4\nJ4frryQwJh8SUoYUBEpKSqiqqkJ7DB4pNTV1RIPMNAj0ZfcqMD6Yep5tvFr9v9DZCklpbudMqdjQ\nHQT6aRMYxoCxpKQkJk+eHILMqUA6gVxfqt4DTxIUz4fSU8HXCZ986HaulIodg5UEwDYOa8Ow6zQI\n9GX/Riicbp/8SxbatD2r3c2TUrGku02gn5IAjHzqCBUSGgT6UrsFCqfZ/fR8yJ8Ke95zN09KxZJg\nSgJZ46HxE+195zINAr11NEPdblsS8Jtwsi0JaPc0pYLT5pQE+msTANtDyNcFzTX9H6PCToNAb4cr\n7TZ/Sk9a0RxoOQiN+1zJklIxp60eEEjO7P+YYXQTVaGnQaC3ut12m1Pakzb2BLut3Rzx7CgVk9ob\nbHuAZ4A/Mf4BY9o47CoNAr0ddpa8y5nYk+YPAjUaBJQKykDzBvn5p47QxmFXaRDorW43JI05ckm8\n9AK7OEbNpv7PU0r1GGjeIL8x+ZCQbBuHlWs0CPTWsNcWU3sPSy+cDjWRX/BBqZjkrw4aiAhkFmlJ\nwGUaBHprqoGMY45OHzvDdh3VHkJKDa69AVIGaBT2i/EVxkYDDQK9Ne2HjLFHpxdMhY4m7SGkVDA6\nmo9eYL4vMb7C2GigQaC3pv2Q2UdJIM+Zs+TQjsjmR6lY1N4EKRmDH5c1HhqqtYTtIg0Cgdqb7NN+\nXyWBvGPt9tD2yOZJqVjU0TTwGAG/rPHgbbcz9ypXaBAI5B+52FebQPYE8CRqSUCpwfh8tjoo2JIA\naJWQizQIBGrcb7d9lQQSEiG3FA5qSUCpAXW2ACbINgH/qGFtHHaLBoFATf4gMK7v9/OmwKGdkcuP\nUrGoo8luk7UkEAs0CATyB4G+GobBtgsc2qGNWEoNpKPZboPpIpo+FsQDjdXhzZPqlwaBQE37bb1/\nWl7f7+dNgc7mnmChlDpae6PdBlMdlJBo2+C0Osg1gwYBEZkgIitFZJOIbBSRW5z074nIXhFZ6/xc\nGHDOHSJSISJbReT8gPSTRGS98959Em2rRTftt9ND9DfplX9mUW0XUKp/Q6kOAh0r4LJgSgJdwLeM\nMTOARcBNIjLDee9eY8xc52c5gPPeFcBMYCnwgIgkOMc/CFwPTHV+lobuVkKg5RCMKej//TwnCGgP\nIaX6110dNJQgoCUBtwwaBIwx1caYD5z9RmAzUDzAKZcATxhj2o0xO4EKYKGIFAFZxphVxhgDPAZc\nOuI7CKWWQzAmt//3syc63US1JKBUv7qrg4YSBLRNwC1DahMQkVJgHuBfcPfrIrJORB4VEf9fz2Jg\nT8BpVU5asbPfOz16tB7qvz0AbP1l9oSehWeUUkcbTnVQR2PPamQqooIOAiKSATwL3GqMacBW7UwB\n5gLVwD2hypSI3CAi5SJSXltbG6rLDq7lEIwZIAiAnT5Cg4BS/RtydZCOFXBTUEFARJKwAeCPxpg/\nAxhj9htjvMYYH/AIsNA5fC8wIeD0Eidtr7PfO/0oxpiHjTFlxpiywsLCodzP8BkDrYcHLgmAHTCm\nQUCp/rU7JYGkIHoHgY4VcFkwvYME+A2w2Rjzs4D0ooDDPgtscPZfAK4QkRQRmYxtAH7PGFMNNIjI\nIueaVwPPh+g+Rq6tHox38JJAbqkNFq11EcmWUjGnoxES02z1aTD8QUDHCrgimP+lxcAXgfUistZJ\n+y5wpYjMBQxQCfw7gDFmo4g8BWzC9iy6yRjjdc77KvA7IA14yfmJDi0H7TaYkgBA3S5IywlrlpSK\nScHOG+SX6TxPanWQKwYNAsaYt4G++vMvH+Ccu4C7+kgvB2YNJYMR03rYboMpCYCtEiqaE84cKRWb\n2puCGyjml5hiu2ZrdZArdMSwn38q22BLAtouoFTfgp1GOpCOFXCNBgG/VicIDFYSSM2GtFydSE6p\n/nQEuaBMoKxiHSvgEg0Cft0lgQEGi/nlajdRpfo11OoggKwirQ5yiQYBv9ZDdjbD1CAae7WbqFL9\n62gKfqCYX9Z4+x3sbA1PnlS/NAj4tRyyAaC/yeMC5ZZC/R7wdoU9W0rFnKH2DgIdMOYiDQJ+rYeD\nqwoCGwR8XVp8Vaov7cNsGAYNAi7QIODXVh98v3/tIaRU34yxg8WG2iaQqQPG3KJBwK+tLrj2ANAg\noFR/utrA+IZRHeQfMKal60jTIODXVm+7fwYju8ROKX1Yu4kqdYT2Ic4g6peSCSnZWh3kAg0CfkMJ\nAp4EyJmoJQGleusY4loCgXTAmCs0CICtxxxKEADtJqpUX/wlgaFWB4EzVkCDQKRpEADbN9nbMbQJ\n4TQIKHU0/1oCWhKIGRoEwJYCYOglAZ1SWqkjDXVVsUBZxdC0H7ydoc2TGpAGARh+EAA7pbRSyuoY\nSXXQeMBA476QZkkNTIMADDMITLZbrRJSqkd376AhjhMAHSvgEg0CEBAEghwxDJA7yW51NlGleoy0\nTQB0rECEaRAAO1AMhlYSSM22aw9oSUCpHiPtIgpQr0EgkjQIwPCqg0B7CCnVW0czJCRDYvLQz03L\ntYvT11eFPl+qXxoEYHglAdAgoFRvHc3Daw8AEIGcCXaGXhUxGgTAlgSSxgz96UWnlFbqSO3DWEsg\nUPYEqNsduvyoQQ0aBERkgoisFJFNIrJRRG5x0vNE5FUR2eZscwPOuUNEKkRkq4icH5B+koisd967\nT0T6WsA+8lrrhl4KAMibrFNKKxVoOAvKBNKSQMQFUxLoAr5ljJkBLAJuEpEZwO3Aa8aYqcBrzmuc\n964AZgJLgQdEJMG51oPA9cBU52dpCO9l+IY6ZYSfziaq1JE6hrG0ZKDsCXYQZntj6PKkBjRoEDDG\nVBtjPnD2G4HNQDFwCbDMOWwZcKmzfwnwhDGm3RizE6gAFopIEZBljFlljDHAYwHnuGvEQUC7iSoF\nDG9VsUA5E+22TksDkTKkNgERKQXmAauBccYY/6iOfcA4Z78YCPwfrHLSip393unua6sPfi2BQFnF\nzpTSlSHPklIxaaRtAv4goFVCERN0EBCRDOBZ4FZjTEPge86TvQlVpkTkBhEpF5Hy2traUF22f23D\nbBPQKaWVOtJI2wSyJ9itNg5HTFBBQESSsAHgj8aYPzvJ+50qHpxtjZO+F5gQcHqJk7bX2e+dfhRj\nzMPGmDJjTFlhYWGw9zJ8w60OAu0mqlSgkbYJZIyz4wy0JBAxwfQOEuA3wGZjzM8C3noBuMbZvwZ4\nPiD9ChFJEZHJ2Abg95yqowYRWeRc8+qAc9wznLUEAmkQUKrHSMYJAHg8duU+LQlETGIQxywGvgis\nF5G1Ttp3gR8BT4nIdcAu4HIAY8xGEXkK2ITtWXSTMcbrnPdV4HdAGvCS8+Oujia7Juqwg8Dkniml\nh7IegVKjjbfLrjGckjmy62RP0IbhCBo0CBhj3gb6689/dj/n3AXc1Ud6OTBrKBkMO/96AMP9Ax44\npbQGARXPOkYwg2ignAmw7dWR50cFRUcM++cNSska3vn+IKCziap4N5IFZQJlT7SLy3S2jTxPalAa\nBPxBYNglAWdKaW0XUPGuexrpEJQEQEfiR4gGgeHOIOqnU0orZXWvKjbCNoHuAWO6al8kaBAYaRAA\n7SGkFIxsVbFA3WMFtHE4EjQItDvj3oYzYtgvb7IGAaVCVR2UVQySoN1EI0SDwEgbhkGnlFYKAhqG\nR1gdlJBoxwronFwRoUGgrd72ZkgIZshEP3JLdUpppULVRRRs6Vp73EWEBoHhzhsUSGcTVaqnTWAk\ns4j65U7W71OEaBAYyZQRfrqugFI9bQJJISoJ+Efiq7DSIBCKIJBVDJ4kDQIqvnU02QDgCcGfldzJ\ndqulgbDTIBCKIKBTSis18hlEA+U5QUDbBcJOg0AoggDoWAGl2kMYBLSdLWI0CGgQUCo0Rrq0ZKCU\nTEgv1JJABMR3EBjpWgKBcku1IUvFt5GuKtZbrg7CjIT4DgIjXUsgUP6xdnto+8ivpVQsCnUQ0LEC\nERHfQSAU8wb5FRxvtwe2jfxaSsWika4q1lvuZDsAs6s9dNdUR9EgAKGrDvIkwoGPR34tpWJRe1Po\n2gTA6SFk4LDOJhpOGgQgNEEgIQnypmgQUPGrozn0bQKgPYTCTIMAhCYIgK0S0uogFY+MgY7G0FYH\n6ViBiNAgACEMAlPh4HadTVTFn64228kilCWB9EJ7Pe1sEVYaBGBkawkEKjgefJ3arU3Fn/YQrS8c\nSMQ+WGnpOqwGDQIi8qiI1IjIhoC074nIXhFZ6/xcGPDeHSJSISJbReT8gPSTRGS98959IiKhv50h\nCsVaAoG6ewhpu4CKMx0hnEE0kFaxhl0wJYHfAUv7SL/XGDPX+VkOICIzgCuAmc45D4hIgnP8g8D1\nwFTnp69rRlZbvZ3waiRrCQTKP85uNQioeBPKtQQCFUyFhqqekoYKuUGDgDHmLeBQkNe7BHjCGNNu\njNkJVAALRaQIyDLGrDLGGOAx4NLhZjpkQrGWQKC0HMgYp08uKv50Ly0ZhpIAwMGK0F5XdRtJm8DX\nRWSdU12U66QVA4GrQ1c5acXOfu/0PonIDSJSLiLltbW1I8jiIEI1ZUSgguO1JKDiTzjaBEAHYUbA\ncIPAg8AUYC5QDdwTshwBxpiHjTFlxpiywsLCUF76SGEJAlNtEDAmtNdVKpp1NNptqKuD8qaAePTB\nKoyGFQSMMfuNMV5jjA94BFjovLUXmBBwaImTttfZ753urnCVBNrqoDmMJRilok27EwRSQ9TJwi8x\nxY7G1yAQNsMKAk4dv99nAX/PoReAK0QkRUQmYxuA3zPGVAMNIrLI6RV0NfD8CPIdGuEIAoXT7LZm\nc2ivq1Q0a2uw21D1tAukPYTCatBuMSLyOHAGUCAiVcCdwBkiMhcwQCXw7wDGmI0i8hSwCegCbjLG\neJ1LfRXb0ygNeMn5cVc4gsDYmXZbswmmnB7aaysVrfwlgZTM0F+7YCpsXwk+r13FT4XUoEHAGHNl\nH8m/GeD4u4C7+kgvB2YNKXfhFMq1BAJljIUx+bB/Y2ivq1Q0a2+wjcLh+CNdcDx426Fud89UEipk\n4nfEcCjXEggkAmNn2JKAUvGivSE8pQDQHkJhFr9BINTzBgUaNxNqtoDPF/prKxWN2hvD0x4AOhI/\nzDQIhLo3A9iSQGcz1FWG/tpKRaO2MJYExuTBmAKo1c4W4RC/QaD1sN2m5YX+2uOcxmFtF1Dxor0x\nPA9UfuNmwH6tYg2H+A0CLc5MGGm5Ax83HIXT7VZ/aVW8CGebAMC4Wbbbtc87+LFqSOI3CPhLAmPC\nUBJIybCrItVoSUDFiXC2CYANAl2tusBMGMRxEPCXBMIQBMBWCWlJQMWLtoYwBwF/FeuGgY9TQxa/\nQaDlECSkQFJaeK4/doZdEamzNTzXVypa+Ly2I0Q42wQKp9s5hLSdLeTiNwi0HrZVQeFa22bcDDsO\nQaePUKNdu3/KiDC2CSSlQv5UDQJhEN9BIFxVQQBFc+y2+qPwfYZS0SCcU0YEGjdTq4PCIH6DQMuh\n8PQM8sudbAeiffJh+D5DqWgQzsnjAo2bCXW7ej5PhUT8BoHWQzAmjEFABIrmQvXa8H2GUtEgYiUB\nZ+oxrWINqTgOAofDWxIAGD/P9hDqag/v5yjlJn+bQDimYAmkPYTCIj6DgDFOdVAY2wQAxs8FX6c2\nZqnRLVLVQdklNtDsWx/ez4kz8RkEOprsH+dwDBQLNH6e3WqVkBrN2ursNi0nvJ8jYjtcaDtbSMVn\nEOieNyjM1UE5kyA1R39p1ejm/z6lhjkIAIyfb0vWWsUaMvEZBFrCPFrYT8SWBj7RkoAaxVrrICkd\nEpPD/1nF820pfp+2C4RKfAYB/5QR4a4OAtsuULMZOtvC/1lKuaGtLvxVQX7j59vtJx9E5vPiQJwG\ngQhVB4EtCWjjsBrNWusiUxUEtnE4fSzsXROZz4sD8RkEIlUdBFCywG73rA7/ZynlhtbDkSsJiNgq\nob1aEgiVQYOAiDwqIjUisiEgLU9EXhWRbc42N+C9O0SkQkS2isj5Aeknich65737RMI1aU8QIlkS\nyBoP2RNhz6rwf5ZSbmiri8x3yW/8fLvUpH+QmhqRYEoCvwOW9kq7HXjNGDMVeM15jYjMAK4AZjrn\nPCAiCc45DwLXA1Odn97XjJzWw5CcEZmGLICJJ8Pu1XZ8glKjTSSrg8CWBDDa4SJEBg0Cxpi3gEO9\nki8Bljn7y4BLA9KfMMa0G2N2AhXAQhEpArKMMauMMQZ4LOCcyIvEQLFAE06Gpn123hOlRptINgyD\nNg6H2HDbBMYZY6qd/X3AOGe/GNgTcFyVk1bs7PdO75OI3CAi5SJSXltbO8wsDiDc8wb1NnGR3e7W\ndgE1ynS1Q2dLZEsC6fl2DE7V+5H7zFFsxA3DzpN9SOs5jDEPG2PKjDFlhYWFoby01VwLYwpCf93+\njJ1hh9Rru4AabVojNFq4t0mfgl3vahVrCAw3COx3qnhwtjVO+l5gQsBxJU7aXme/d7o7mg9AxtjI\nfZ4nAUrKtCSgRp/uKSMiWLIGGwRaDsCBbZH93FFouEHgBeAaZ/8a4PmA9CtEJEVEJmMbgN9zqo4a\nRGSR0yvo6oBzIssYaKqB9DCUMAYyYRHUbOp5clJqNPD/PkeyOghg4qfsdtc7kf3cUSiYLqKPA+8C\n00SkSkSuA34EnCsi24BznNcYYzYCTwGbgJeBm4wxXudSXwV+jW0s3g68FOJ7CU57A3jbIx8EJi4C\nDOx5L7Kfq1Q4RWryuN7yj7WDxnb9M7KfOwolDnaAMebKft46u5/j7wLu6iO9HJg1pNyFQ/MBu41k\ndRDAhIV2Yfudb8Lx50X2s5UKl0iOuQkkYquEdr8b2c8dheJvxHCT03wR6ZJAUpotDex4I7Kfq1Q4\nNTu998Z9XSKSAAAbVElEQVTkR/6zJy2G+j1Qtzvynz2KxF8QaHaCQKRLAgBTzrCrIvkDkVKxrvkA\neJLCv6pYXyadYrdaJTQicRgEnCeXSJcEwAYBgB1vRv6zlQqHlgOQXmCrZyJt7AwbfCrfjvxnjyLx\nFwSaagGJ7DgBv6I5theFVgmp0aL5oDvfJbBdr0tPs98nHS8wbPEXBJpr7DoCCYO2iYeeJwEmL9Ff\nWjV6tBywI3jdctzZtl3gwMfu5SHGxWEQqHWnKshvyhnQUAUHK9zLg1Kh4vb36Vink2LFCvfyEOPi\nLwg07nenUdjv2LPsdtvf3cuDUqHiZnUQQO4kKDheg8AIxGEQqIbM8e59ft5kKDwBtrozVk6pkOlq\nh45Gd6uDAI47x/YQ6mx1Nx8xKr6CgM9ng0CWi0EAYPqF9pe2pfcM3UrFEP/ASzdLAmDbBbraoFKn\nkBiO+AoCzbXg63I/CEy7CIxXq4RUbGtxgkC6y0Fg0mJITIWKV93NR4yKryDQ+IndZha5m4/x8yDj\nGNi63N18KDUS/pKAmw3DYEfjT14CW5Zrr7thiK8g0OAEgSyXg4DHA9OWQsVrtl5VqVgULdVBADMu\ngfrdUK1LTg5VnAaBfhc1i5xpF0FHkw4cU7Gru2R9jLv5AJh2IUgCbHJnhvpYFl9BoLHa/qK4XXwF\nmHK6HfK+/hm3c6LU8DRU2xXzUjLczokdADp5iQ0CWiU0JPEVBBo+sU8tngS3cwKJKTDzs7DlRWhv\ncjs3Sg1d4yfut68FmnEJHNoB+ze6nZOYEmdBYG90/dKeeLldpFsbiFUsaqh2v30t0PSLQTxaJTRE\n8RUE6nbbEYbRYuIpkFUC655yOydKDZ3bAy97yyiE0lNh/dNaJTQE8RMEvF1QtwdyoigIeDxw4udh\n++vO7KZKxQifFxr3RVdJAGDuVXB4p64xMATxEwQa9toBWrmlbufkSLO/YPO17km3c6JU8Jpr7e9t\nNFWvApzwGUjOhA//4HZOYkb8BIHDlXYbbUFg3AyYcDKUP2qntVAqFnR3t46i6iCA5DEw6zLY9By0\nN7qdm5gwoiAgIpUisl5E1opIuZOWJyKvisg2Z5sbcPwdIlIhIltF5PyRZn5IojUIACy4Hg5th51v\nuJ0TpYLTWG230VYSAJj3RdvhYuNf3M5JTAhFSeBMY8xcY0yZ8/p24DVjzFTgNec1IjIDuAKYCSwF\nHhCRyPXVrNsFnsToGCjW24zP2FGX7//G7ZwoFZz6KruNxu9TSRkUTIPy32oDcRDCUR10CbDM2V8G\nXBqQ/oQxpt0YsxOoABaG4fP7drgSskvcWVFsMIkpMP9q21XU/+VSKpodroSkdHfX5uiPCCy8Hj75\nAPa853Zuot5Ig4ABVojIGhG5wUkbZ4xxyorsA8Y5+8XAnoBzq5y0yDi0IzqrgvzKvgQIrHrQ7Zwo\nNbhDO+33yY0F5oMx91/tet7v3u92TqLeSIPAqcaYucAFwE0isiTwTWOMwQaKIRGRG0SkXETKa2tD\n0HXS54Paj6Fw+sivFS45E2130fJH7WpNSkWzwzvtAknRKjkdyr4Mm1+0D4CqXyMKAsaYvc62BvgL\ntnpnv4gUATjbGufwvcCEgNNLnLS+rvuwMabMGFNWWBiCeX4aqqCzGQqnjfxa4XTqN+3qSKu1NKCi\nz4Gmdv7yYRXff349nQd2sK45l5rGNrez1b+FN9h2wFUPuZ2TqDbsICAi6SKS6d8HzgM2AC8A1ziH\nXQP4x3C/AFwhIikiMhmYCkSmwq52q90WnhCRjxu2sdPhhE/D6v+F1jq3c6MUAFv3NXLz4x+y6P+9\nxjee/IiVa9aTZDp4akcip/14JU+V7xn8Im7IKoLZl8MHy+wUF6pPIykJjAPeFpGPsH/M/2aMeRn4\nEXCuiGwDznFeY4zZCDwFbAJeBm4yxnhHkvmg1W6x22gvCQAs+Ta0N8C7v3I7JyrOHWru4L+eW88F\nv3iLlVtquPqUUv76tVN5/UsTAbjpsnMpK83lO8+s48n3d7uc234s+U+7muBbP3U7J1Fr2F1ljDE7\ngDl9pB8Ezu7nnLuAu4b7mUNRXd/KWx/X8n7lYS7e+RZzPHnc91o1n53nYXZJTiSyMDxFc+zsov+8\nH066FrKjsAueGtU6unw89m4lv3htGy0dXr64aBK3nnM8uenJ9oAPKwEoKp3OsrmT+dLv3ue/ntvA\nicU5zBif5Vq++5Q3GeZfY0sDn/p6dLdjuGRUjhg2xvD5B9/ltmfX88bWGoo7KtmTMIE/rd7NZ375\nDrc9s472rsgUQoblnO/ZIfmv/39u50TFEWMMr23ez9Kfv8X//dtm5k3M5eVbTuP7l8zqCQAANZsh\nMQ1yS0lM8HD/lfPISk3iv55bj88Xhf3yl/ynbRt444du5yQqjcogICL88LITeeXWJbx/+xKmmkpm\nLzidNf99LjeefixPlu/hy797n46uKJ2mIbcUFv0HfPQ47P3A7dyoOPDx/kaufvQ9rltWDgK/vXYB\ny760gKnjMo8+eP8G237lrMuRMyaZ2y+Yzge763hpw74I5zwIWUX2+7TuSdi9yu3cRJ1RGQQAlhxf\nyLRjMpH9G8HbAcXzyUhJ5PYLpnP3v8zhnYqD/M/zG9zOZv9O+xakj4W/3gLeTrdzo0apA03tTr3/\nP/hoTx3/c/EMXrl1CWdOH4v0NwZg/0YYN/OIpMvmlzClIJ1frqzAROMo3dO+badtf/Gb+n3qZdQG\ngW5719ht8UndSZ8/qYT/OONYnnh/D3/fGIVPLmCXnrzobti3Dt75hdu5UaNMW6eXX62s4IyfvsHj\n7+3hqpMn8uZ/nsmXT51MUsIAfxaaauwMomOPDAIJHuE/zjiWzdUNvF1xIMy5H4aUDLjgR1CzEVZr\nl9FAoz8IfPKhnZcne8IRyd8453hOKMriv5/fQGtHlLYPzLjE/rz5Y6jZ4nZu1Cjg9RmeXVPFmXe/\nwU9f2copx+bz928s4Qe96/3741+6sVdJAOAzc8eTl57MH1dFaU+h6RfD1PNh5f+DA9vczk3UGP1B\nYM9qWwroVbRNTvTw/c/MZH9DO4++s9OlzAXhwrshJROe+RJ0tLidGxWjvD7Dcx/u5dx73+RbT39E\nYWYKT9ywiEeuLuPYwiEsFF/9kd2Om3XUWymJCfzLSSW8unk/+xuicBCZCHz6F5CYCs9eB10dbuco\nKozuIFBfBQcrYPKSPt9eODmPc2eM48E3tnOgqT3CmQtSxli47GHbI2P5t93OjYoxXV4ff/mwinPv\nfZNbn1xLcoKHB66az3NfXcyiKflDv+DuVZB/HKT3fe6VCyfi9RmefD+KB5Bd8ksbzF7/gdu5iQqj\nOwjseNNup5zR7yG3LZ1OS0cX//vm9ohkaViOOwdO/w6s/aNON62C0tTexW/f2clZ97zJN578iOQE\nDw9eNZ/lN5/GhScW4fEMY+I3n8+WrCcu6veQ0oJ0Fh+Xz5Pv74nO7qIA0y+CsuvsWJz1z7idG9dF\n4bzKIbTjDUgvhLEz+j3kuLEZXDqvmN+v2sX1S6YwNjM1cvkbitNvs+0by79t53CfttTtHKkoVHW4\nhWX/rOSJ9/bQ2N7F/Ik5fPfC6Zw345jh/eEPdHAbtB6CiZ8a8LArFkzk649/yNsVB1hyfAjm/gqH\npT+yMwk891XbJbukbNBTRqvRWxLwdkLFCphypl3QfQA3nzWVTq/hf9+M4tkGPQnw+d/aEcVPXwu7\nV7udIxUlOr0+Xtm4j+t+9z5LfrKSR9+p5IzpY/nLVz/Fn7+6mKWzhvnk31vl23Y7QEkA4LyZ48gd\nk8QT0TqVBEBiMlz+e8g8Bv70hZ75xeLQ6A0CO960Ty0zLx300NKCdC6dW8wfVu2iJhobtPxSMuBf\nn7b1mn+4rOdLqeJSRU0jP1y+mVN++Dr//vs1rN9bz42nH8s/vnMm9185j3kTcwe/yFB8/LJ9as6b\nMuBhKYkJXDa/hFc37Y/etjaw7Rr/9iyIB5Z9Gg5UuJ0jV4zeILDxL5CSZevTg3Dz2cfR5TM8GM1t\nAwAZhXDtclsl9IfPw8d/dztHKoJ2Hmjml69vY+nP3+Kcn73Fr9/eybyJOfzmmjL+eftZfGfpdMbn\npIX+g9sb7YPVtIuCWkjmyoUT6PQa/vxBlK+UVzAVrvkr+Lyw7GLYt97tHEXc6AwCxsDectsAlJgS\n1CmT8tO5bF4xf1y9Ozq7twXKKoIvLbe/wI9/Af75S11LdZQyxrC5uoFfrazgovv+wZl3v8Hdf/+Y\njJRE7vz0DN694yweubqMs08YR+JAg7xGqmIFeNth+oVBHX7c2EzKJuXyxPt7onMEcaCx020gQODR\nC+y9xhGJ9v+gsrIyU15ePvQTfV47JXNa8EXi3QdbOOueN/i3RZP43meOHgwTddqb4LkbYfNfYdbn\n4aJ7IC2KZ0hVQWlu7+Kf2w/y+pYa3thaQ3W9fSiZOyGHi2cXcdHsIoqyw/C0P5A/fM4OFLt1Q9Dr\ndD9dvof/fGYdj1+/iFOOHUZ31Ehr+AT+eLkdVXzGHXbqFmd+pFgkImuMMYO2eI/e3kGehCEFAICJ\n+WP43PwS/vTebm48/ViOyY7SnkJ+KRnwL4/B2/fAyh/a7nuX/HLALrEq+rR0dPHBrjpW7TjIqh0H\n+aiqjk6vISMlkVOPK+Ab54zljGmFjM1y6ffx0E6oeM32UAsyAAB8es54fvzyFh56c3tsBIGs8fDl\nl+z8Qivvgsp/wKfvG/XTT4/eksAw7TnUwtn3vMl5M8fxy3+dH7HPHbGqcvjz9XY91ZmfhfP+L2SX\nuJ0r1Ysxhk/q21i3p461VXWsqTzc/Uc/wSOcWJzNyVPyOH1qIWWleSQnRkGN7d++BeW/hW9ssH8o\nh+BXKyv46Stb+evXTuXEkuwwZTDEjIEPfw8v32EXpDn9Njjla7ZHUQwJtiSgQaAPv3x9G3f//WMe\n+reTWDrrmIh+9oh0tsI798HbPwMEFlxnF9LIjKF7GEWMMVTXt7F1XyMb9tbzUVUda/fUd/eYSU7w\nMLM4i0VT8jl5ch5lpXlkpERZ4fzgdvjVQrswy8U/G/LpDW2dnPHTN5g6NoMnbljU/8yk0ah+L7x8\nm61uzS21weDEy4dUGnKTBoER6PT6uOSX77C/oY3nv7aYktwxEf38ETu8y06Stf5pu5jG3CvtKmXj\n57mds1HJGMOBpg52Hmhm674Gtuxr5OP9jWzZ10hjW1f3cccWpjNnQg5zJ+QwpySH6UWZpCRGcZ2z\nzwvLPgPVa+HrH0DmuGFd5vfvVvLfz2/k51+Yy6XzYnClvG0r7BQT1R9B3rFw8o0w5wpIjbJV1HrR\nIDBCFTVNfPZX7zA+J40/XX8y+RnB9TKKKod2wNv3wrqnoasVjpkNc6+yvaZyJgx+vurm9RkONLWz\n62ALlQeb2XWwmcoD/v0Wmtp7/thnpiYy/ZhMph2TybRxmUw7JovpRZlkpSa5eAfDsOL7tlR56YMw\n91+HfZkur48vPLyKj/c18sLXT2VyQXoIMxkhxsCWF+EfP4NPPoDkDFvtOusyKF0SlaUDDQIh8E7F\nAb78u/cpzk3jgavmM/2Y6I78/Wqts6WCD5b19IMumgPHL4VJi6FkASTHWGknhNq7vBxq7qCmoZ3q\n+jaq61udbRvVdXZ/f0MbXQFz4SR6hAl5Y5iUP4bS/HRK88cwqSCdaeMyKcpOja1qj958PnjzR3YK\n8/nX2Jk3R3g/ew618Jlfvk16SiJP3LAo9krXgfausXN4bXoBOhphTL79Lk05w05WGSXVr1EbBERk\nKfALIAH4tTHmRwMd72YQAFi94yA3/elDGto6ufH0Y7nu1Mlkp8XYE12gAxX2iWbLi/aX2fjAkwTF\n82110bhZcMwsKDwBkqK8d1QfOr0+Glo7aWjrcrad1LV0cqi5g4NN7RxwtgebOjjY3MGBpvYjqmz8\nkhM9FGWnUpSdyvjsNI7JTqUoJ40JuWlMLkhnfE7awIuvxKraj+Gl78COlbbU+Jn7Q9ZNcl1VHVf9\nejVJCR5+cMlMLjqxKLaDZWebHVOw8c+w/XVoPWzTC6bZ6evHz7UPW+Nm2Z58ERaVQUBEEoCPgXOB\nKuB94EpjzKb+znE7CIBdgu/O5zfyt/XVjElO4JwTxnHm9ELKJuVRkpsWu7/IbfV2DqJd78Dud2Hf\nBuhstu+Jxy7Ek1va85Mz0U5tnTHOTsyXljuiJ0RjDB1eH22dPlo7vLR0dNHS4XV+evZbO7po9qe3\nd9HSabeBf+jrWztpaO2itbP/BYJEIG9MMvkZyeSnp5CfkUxBRgr56cnkZ6RQkJHM+Jw0irJTyUtP\njt3/16FqOWT/iG14Fra+ZKs6zv2enWkzxP8G22ub+NqfPmRzdQNzSrK5fMEEzpw2NjyjnCPJ57Or\nAO54w07nUr3WrsAGgNgR/nmT7ZQbeZMhZ5ItMWSMs9+plD7Wch6haA0CpwDfM8ac77y+A8AY88P+\nzomGIOC36ZMGfr9qFy9tqKauxa5Tmp6cwMT8dCbkppGXnkzOmGRyxiSRk5ZESpKHlMQEkhM8pCR5\nSE7wkJzoIcEjeJwvl0cEj8duBRARPNKzNQYM4DPG7huDwab703zO/2F32gDH+3yGLp/B27319bzu\n8pLWtJvMui1kN2wlvWUPGS1VZLbuZUznoaP+PbySSFNiHi0JmbTKGJo96bRIOs2k0cgYGk0azSaF\nFpNEizeRFl8Szb5EWryJNHkTaPIm0k4S7STRRQJek4AXT8DPka+7SCAtKZExyQmMSUkgOy2JrFTn\nJy2RrNQkm5bW8zrLOSY/I5ncMckkhGIitWjm89k1tX2ddhJFb6cdNNnWAG11NvA319rOA3W77DoV\nh5ypUtLHwknX2IbP9IKwZdHrMzxVvodH397JtpomAMZnp3LcuEwm5qUxPiet+/82Oy2JjNTE7u9O\nUoKHpAQhOcHuezyCCN3fHbsFQbrjV+Br/3Gm1/eEgP0j0g0YDL6A75LPOcj0PscX8D30+fA07yex\n5iOSazeSWLeT5PpKkhsqSWw7ePS/SeIY2lML6EzOoTM5m46kLDqSspnwhXtISBle1Vm0BoHPA0uN\nMV9xXn8RONkY87X+zommIODn9Rm27mtkze7DbK9pYvehFqoOt1DXYqseOrw+t7MYcmNoo1gOUCD1\nHOOpp9DTwFipp1DqyZZmsqSVTFpIp5V000y6aSGRo6tZRk5sjydPgt2Kx6aJ857IMLeDnN/N+b50\nf28Cvj+90474bplehw/1vL7SfPaPvK/L2XbatGAkpkHuJLtATPFJMOlTtm0ogiNk7ZQYjazeeZA1\nuw5TebCZPYdaqW8d3QvBZ9LCeDlAodQzlsMUSj2FUsdYqSOHJrKkmRyayZJmxnx3B6kpw+uUEtMj\nhkXkBuAGgIkTJ7qcm6MleIQZ47OYMf7ohmJjDK2dXupbO2nv9NHe5aO9y0tHl93v6PLZJ3Jz5NN6\n92v8rw0+n/0b5BHnKSbgSSew5OB/wvEf53ESPH0cn+AREhOEBI+HRI/Y191bDwkJ0ne68zroKYmN\nga42O3ahq835aT9y2xmQ7usC47Vbn9f+9Ps6MM2HfSwzQ9wytOMDA0F3FYn0et1X2lDP6/V6wGsJ\nJCTZNp2ERGebZINjQrLdT0iyEymmZvds0wtsdZ7L1V0iPd+jLy3uGZXb0tFFQ2sXDW2dNLR20tje\nRWeXj06vodPro8Prs9suH96Ap29Dz1N54Gv7/pHv2e8E3VV+HgksKdBdUu8umQfs43yXPL2+e/Q6\nzn+dBI+QIPa7k+C89nSncVRagkdoEaHdI0xNCv8AtUgHgb1AYN/EEiftCMaYh4GHwZYEIpO10BAR\nxiQnMiY5KuNr5IhAUpr9UWoI/N+fqJ+2ZZSIdPeG94GpIjJZRJKBK4AXIpwHpZRSjog+rhpjukTk\na8Ar2C6ijxpjNkYyD0oppXpEvM7CGLMcWB7pz1VKKXW0UTjaRSmlVLA0CCilVBzTIKCUUnFMg4BS\nSsUxDQJKKRXHon4qaRGpBXYN8/QC4EAIs+MmvZfoM1ruA/ReotVI7mWSMaZwsIOiPgiMhIiUBzN3\nRizQe4k+o+U+QO8lWkXiXrQ6SCml4pgGAaWUimOjPQg87HYGQkjvJfqMlvsAvZdoFfZ7GdVtAkop\npQY22ksCSimlBjAqg4CILBWRrSJSISK3u52fwYjIoyJSIyIbAtLyRORVEdnmbHMD3rvDubetInK+\nO7num4hMEJGVIrJJRDaKyC1Oeszdj4ikish7IvKRcy/fd9Jj7l7ArvEtIh+KyIvO61i9j0oRWS8i\na0Wk3EmL1XvJEZFnRGSLiGwWkVMifi92xZ3R84Odono7MAVIBj4CZridr0HyvASYD2wISPsJcLuz\nfzvwY2d/hnNPKcBk514T3L6HgHwXAfOd/UzgYyfPMXc/2IWiMpz9JGA1sCgW78XJ3zeBPwEvxvjv\nWCVQ0CstVu9lGfAVZz8ZyIn0vYzGksBCoMIYs8MY0wE8AVzicp4GZIx5C+i9kvsl2F8QnO2lAelP\nGGPajTE7gQrsPUcFY0y1MeYDZ78R2AwUE4P3Y6wm52WS82OIwXsRkRLgIuDXAckxdx8DiLl7EZFs\n7APgbwCMMR3GmDoifC+jMQgUA3sCXlc5abFmnDGm2tnfB4xz9mPm/kSkFJiHfYKOyftxqlDWAjXA\nq8aYWL2XnwPfAQJXoo/F+wAbiFeIyBpnPXKIzXuZDNQCv3Wq6X4tIulE+F5GYxAYdYwtC8ZUNy4R\nyQCeBW41xjQEvhdL92OM8Rpj5mLXw14oIrN6vR/19yIiFwM1xpg1/R0TC/cR4FTn/+QC4CYRWRL4\nZgzdSyK2GvhBY8w8oBlb/dMtEvcyGoNAUIvZx4D9IlIE4GxrnPSovz8RScIGgD8aY/7sJMfs/QA4\nxfSVwFJi714WA58RkUps9ehZIvIHYu8+ADDG7HW2NcBfsFUisXgvVUCVU7oEeAYbFCJ6L6MxCIyW\nxexfAK5x9q8Bng9Iv0JEUkRkMjAVeM+F/PVJRARbx7nZGPOzgLdi7n5EpFBEcpz9NOBcYAsxdi/G\nmDuMMSXGmFLs9+F1Y8y/EWP3ASAi6SKS6d8HzgM2EIP3YozZB+wRkWlO0tnAJiJ9L263jofjB7gQ\n2ytlO/B/3M5PEPl9HKgGOrFPB9cB+cBrwDZgBZAXcPz/ce5tK3CB2/nvdS+nYouv64C1zs+FsXg/\nwGzgQ+deNgD/46TH3L0E5O8MenoHxdx9YHv9feT8bPR/v2PxXpy8zQXKnd+x54DcSN+LjhhWSqk4\nNhqrg5RSSgVJg4BSSsUxDQJKKRXHNAgopVQc0yCglFJxTIOAUkrFMQ0CSikVxzQIKKVUHPv/Aayi\ni8eqxqz3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ac7cba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(times, rabbits, label='rabbits')\n",
    "plt.plot(times, foxes, label='foxes')\n",
    "plt.legend(loc='best') # put the legend somewhere it doesn't overlap things\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's copy and paste the code above, and turn it into a function so we can evaluate it with many different step sizes. In order to compare the different results we need to pick some quantity of interest (QOI). We'll choose the maximum number of foxes, because it's simple to find and presumably interesting, although if you care more about the time at which the maximum occurs, or the number of rabbits, or something else, then use that instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "k1 = 0.015\n",
    "k2 = 0.00004\n",
    "k3 = 0.0004\n",
    "k4 = 0.04\n",
    "\n",
    "def solve_by_euler(step_size=1.):\n",
    "    \"\"\"\n",
    "    Evaluate by simple Euler, with the given step size.\n",
    "    \n",
    "    Returns the peak number of foxes\n",
    "    \"\"\"\n",
    "    end_time = 600.\n",
    "    times = np.arange(0, end_time, step_size)\n",
    "    foxes = np.zeros_like(times)\n",
    "    rabbits = np.zeros_like(foxes)\n",
    "    rabbits[0] = 400.\n",
    "    foxes[0] = 200.\n",
    "    for n in range(len(times)-1):\n",
    "        delta_t = times[n+1] - times[n]\n",
    "        rabbits[n+1] = rabbits[n] + delta_t * (k1 * rabbits[n] \n",
    "                                                - k2 * rabbits[n] * foxes[n])\n",
    "        foxes[n+1] = foxes[n] + delta_t * (k3 * rabbits[n] * foxes[n]\n",
    "                                          - k4 * foxes[n])\n",
    "    return foxes.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "run      step size      max foxes\n",
      " 0           10        9998.8\n",
      " 1         6.31        3757.4\n",
      " 2         3.98        6831.8\n",
      " 3         2.51        4686.9\n",
      " 4         1.58        3801.7\n",
      " 5            1        3379.3\n",
      " 6        0.631        3146.0\n",
      " 7        0.398        3012.6\n",
      " 8        0.251        2932.9\n",
      " 9        0.158        2884.2\n",
      "10          0.1        2854.1\n",
      "11       0.0631        2835.4\n",
      "12       0.0398        2823.7\n",
      "13       0.0251        2816.4\n",
      "14       0.0158        2811.7\n",
      "15         0.01        2808.8\n",
      "16      0.00631        2807.0\n",
      "17      0.00398        2805.8\n",
      "18      0.00251        2805.1\n",
      "19      0.00158        2804.7\n"
     ]
    }
   ],
   "source": [
    "step_sizes = []\n",
    "maximums = []\n",
    "print(\"run      step size      max foxes\")\n",
    "for i in range(20):\n",
    "    step_size = 10**(1-i/5.)\n",
    "    maximum = solve_by_euler(step_size)\n",
    "    step_sizes.append(step_size)\n",
    "    maximums.append(maximum)\n",
    "    print(\"{0:2d}     {1:8.3g}        {2:5.1f}\".format(i, step_size, maximum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot on a semi-logarithmic plot (log scale on the x axis) and see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEOCAYAAABB+oq7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VdW99/HPN2EIYVYkF0EEFVRATUpEtF4rxRGca0Xt\nvQ5QaSvW2lGo99b21dtH2j5Pa7XVlpbUqXXAWkVxSATRWydEUSE4MCqDCRFFIBiG5Pf8cdaBk5Dh\nJNkn5wR+79drv84+6+y118omnF/2WnutJTPDOeeci0JWuivgnHNu3+FBxTnnXGQ8qDjnnIuMBxXn\nnHOR8aDinHMuMh5UnHPORcaDinPOuch4UHHOORcZDyrOOeci40HFOedcZDqkuwJtrU+fPjZo0KB0\nV8M559qV119//WMzO6ip4/a7oDJo0CAWLlyY7mo451y7IumDZI7z5i/nnHOR8aDinHMuMh5UnHPO\nRcaDinPO7QdKS0sZMWIEpaWlKS0nZUFFUpGkDZKWJKQdIKlE0rLw2jvhs2mSlkt6T9KZCekjJS0O\nn90mSSG9s6QHQ/qrkgal6mdxzrn2rLKyknHjxrF06VLGjx9PZWVlyspK5Z3KXcBZddKmAnPNbAgw\nN7xH0jDgUmB4yHOHpOyQ507gGmBI2OLnnAR8amZHAL8Ffpmyn8Q559qxiRMnsmHDBsyM8vJyJk2a\nlLKyUhZUzOwF4JM6yecDd4f9u4ELEtIfMLPtZrYKWA6MktQP6GFmr1hs3eN76uSJn+thYGz8LsY5\n51xMUVERTzzxBFVVVQBUVVXx+OOPU1RUlJLy2rpPJc/MPgr7ZUBe2O8PrEk4bm1I6x/266bXymNm\nu4DPgAPrK1TSZEkLJS2sqKiI4udwzrl2Ydq0aWzbtq1W2rZt25g2bVpKyktbR32487A2KmuGmRWa\nWeFBBzU5INQ55/YZt9xyCzk5ObXScnNzmT59ekrKa+ugUh6atAivG0L6OuCQhOMGhLR1Yb9ueq08\nkjoAPYGNKau5c861QxMnTqSwsHD3+5ycHM4991yuvvrqlJTXZFCR9FVJ3cP+f0l6RNIXWljebODK\nsH8l8FhC+qXhia7BxDrkF4Smss2SRof+kivq5Imf62JgXrj7cc45l+Dyyy8HQBJ5eXnMnDkzZWUl\nc6fy32a2RdLJwGnATGJPZDVK0v3Ay8CRktZKmgRMB06XtCycazqAmZUCDwFLgaeBKWZWHU51LfAX\nYp33K4CnQvpM4EBJy4HvEZ4kc845V9unn34KwNFHH82cOXPo2rVryspKZkLJ+Jf7eGCGmc2R9D9N\nZTKzyxr4aGwDx/8C+EU96QuBEfWkVwFfbaoezjm3vysrK6NXr14pH/gIyd2prJP0J2AC8KSkzknm\nc845lwHKysr4t3/7tzYpK5ngcAnwDHCmmW0CDgB+mNJaOeeci0xZWRl5eXlNHxiBZILKz4H3zGwZ\ngJl9ZGbFqa2Wc865qGTanco7wIwwv9Y3JfVMdaWcc85FJ6OCipn9xcy+SOxx3kHA25L+LmlMqivn\nnHOudbZt28aWLVsyJ6gAhMkdjwrbx8BbwPckPZDCujnnnGul8vJygDYLKk0+Uizpt8A5wDzg/5jZ\ngvDRLyW9l8rKOeeca52ysjIgg4IK8DbwX2ZW3wT8oyKuj3POuQhlXFAxs79K6i1pOJCTkP6CmX2W\n0to555xrlYwLKpK+DnyH2GSObwKjiU2/8uXUVs0551xrlZWVkZWVRVvN0J5MR/13gOOBD8xsDFAA\nbEpprZxzzkWirKyMgw46iOzs7KYPjkAyQaUqzLOFpM5m9i5wZGqr5ZxzLgptOUYFkuuoXyupF/Ao\nUCLpU+CD1FbLOedcFNpyihZIrqP+wrD7U0nPEVsM6+mU1so551wkysrKOOqoo9qsvAaDiqQD6kle\nHF67AZ+kpEbOOeciYWYZ1fz1OrE15AUMBD4N+72AD4HBKa+dc865Ftu0aRM7duxo06DSYEe9mQ02\ns8OAZ4FzzayPmR1IbHS9z1LsnHMZrq2naIHknv4abWZPxt+Y2VPASamrknPOuSi09cBHSC6orJf0\nX5IGhe0mYH1rCpX0HUlLJJVKuiGkHSCpRNKy8No74fhpkpZLek/SmQnpIyUtDp/dJkmtqZdzzu1L\nMjWoXAYcBPwTeCTsN7T+fJMkjQCuITZv2HHAOZKOAKYCc81sCDA3vEfSMOBSYDhwFnBHmDUZ4M5w\nriFhO6ul9XLOuX1NOoJKMo8Uf0JsVH1UjgZeNbNtAJKeBy4CzgdODcfcDcwHbgzpD5jZdmCVpOXA\nKEmrgR5m9ko4zz3ABcBTEdbVOefarbKyMjp16kSvXr3arMyk1lOJ2BLg3yUdKCkXGAccAuSZ2Ufh\nmDIgPlqnP7AmIf/akNY/7NdNd845x57R9G3ZM5DMiPpImdk7kn5J7AmySmKTVFbXOcYkWVRlSpoM\nTAYYOHBgVKd1zrmM1tZjVCA9dyqY2UwzG2lmpxAb//I+UC6pH0B43RAOX0fsTiZuQEhbF/brptdX\n3gwzKzSzwraaqdM559It44KKpDMl3SlpdtjulNTqznBJfcPrQGL9KX8HZgNXhkOuBB4L+7OBSyV1\nljSYWIf8gtBUtlnS6PDU1xUJeZxzbr/X1vN+QePTtNwKDAXuYU/fxQDgeklnm1lrOu//IelAYCcw\nxcw2SZoOPCRpErEJKy8BMLNSSQ8BS4Fd4fh4c9m1wF1AF2Id9N5J75xzQHV1NRUVFW1+p9JYn8o4\nMxtaN1HSg8Saq1ocVMzs3+tJ2wiMbeD4XwC/qCd9ITCipfVwzrl9VUVFBTU1NRnV/FUl6fh60o8H\nqlJUH+eccxFIxxQt0PidylXAnZK6s6f56xDgs/CZc865DJWOgY/QSFAxszeAEyT9G3vGf6wzs7I2\nqZlzzrkWy7igAhCeqjqUPUGlg6RyM4tsDIlzzrnoxYNKJj39dQZwB7CMPeM/BgBHSLrWzHz6e+ec\ny1BlZWV0796drl27tmm5jd2p/A44zcxWJyaGsSJPEpvDyznnXAZKx8BHaPzprw7Unlsrbh3QMTXV\ncc45F4V0BZXG7lSKgNckPcCeCR0HAhOAmamumHPOuZYrKyvjmGOOafNyG1tO+BbgcmLr0p8YNoCv\nhc+cc85lqHRM0QJNPP1lZu8A78TfSzowjHx3zjmXoaqqqti0aVNm9alImi6pT9gfKWkl8IqkDyR9\nqc1q6JxzrlnSNZoeGu+oH29mH4f9/wtMCEv9ng78v5TXzDnnXItkalDpICnePNbFzF4DMLP3gc4p\nr5lzzrkWSddoemg8qNwBPCnpy8DTkn4n6UuSfkZstUbnnHMZKJ1BpbG5v26XtBj4FrF1VToQWyDr\nUeB/2qZ6zjnnmiseVPr27dvmZTf19Nd8YH6b1MQ551wkysrK6NOnDx07tv049bSsUe+ccy510jWa\nHjyoOOfcPseDinPOuchkfFCRdJKkyyVdEd9aU6ik70oqlbRE0v2SciQdIKlE0rLw2jvh+GmSlkt6\nT9KZCekjJS0On90W1n9xzrn9lplldlCRdC+xwY8nE1uf/nigsKUFSuoPXA8UmtkIIBu4FJgKzA0D\nLOeG90gaFj4fDpwF3CEpO5zuTuAaYk+lDQmfO+fcfmvLli18/vnnaZn3C5p4+isoBIZFvNpjB6CL\npJ1ALrAemAacGj6/m9hTZzcC5wMPmNl2YJWk5cAoSauBHmb2CoCke4ALgKcirKdzzrUr6RyjAsk1\nfy0BIqudma0jdufzIfAR8FlYRTLPzD4Kh5UB8TDbnz1T70NsjZf+YVtbT/peJE2WtFDSwoqKiqh+\nFOecyzjpDirJ3Kn0AZZKWgBsjyea2XktKTD0lZwPDAY2AbMk/UfiMWZmkiK7MzKzGcAMgMLCwijv\nuJxzLqOkc94vSC6o/DTiMk8DVplZBYCkR4CTgHJJ/czsI0n9gA3h+HXAIQn5B4S0dWG/brpzzu23\nMv5Oxcyej7jMD4HRknKBz4GxwEKgErgSmB5eHwvHzwb+Luk3wMHEOuQXmFm1pM2SRgOvAlcAt0dc\nV+eca1fKysro0KEDBxxwQFrKbzCoSPqXmZ0saQuQ2GQkYi1UPVpSoJm9Kulh4A1gF7CIWNNUN+Ah\nSZOAD4BLwvGlkh4Clobjp5hZdTjdtcBdQBdiHfTeSe+c26/FV3zMykrPMERF+1BX5issLLSFCxem\nuxrOOZcS48ePp7y8nKi/5yS9bmZNDifxEfXOObcPSefAR/Cg4pxz+xQPKs455yJRU1NDeXl5ZgcV\nSV0lZYX9oZLOk9T2k/Q755xr1MaNG6murs7soAK8AOSEObuKgf8k9sSVc865DBIfo5Kueb8guaAi\nM9sGXATcYWZfJTa5o3POuQyS7oGPkGRQkXQi8DVgTkjLbuR455xzaZDuKVoguaByA7EZhP8ZBiIe\nBjyX2mo555xrrky4U0l2mpbnw7QqmNlKYuuhOOecyyBlZWXk5ubSrVu3tNUhmae/TpS0FHg3vD9O\n0h0pr5lzzrlmiY9RSeciuMk0f90KnAlsBDCzt4BTUlkp55xzzZfugY+Q5OBHM1tTJ6m63gOdc86l\nTXsJKmsknQSYpI6SfgC8k+J6Oeeca6b2ElS+CUwhtlTvOiA/vHfOOZchduzYwcaNG9MeVJJ5+utj\nYmNUdpPUKWU1cs4512wbNsQWy013UEnm6a/5kgYlvD8eeC2FdXLOOddMmTBFCyS3Rv0twNOSbiPW\nBDYOuDqltXLOOdcsmTDwEZJr/npG0jeBEuBjoMDMylJeM+ecc0nLhClaILnmr/8Gbic2NuWnwHxJ\n41taoKQjJb2ZsG2WdIOkAySVSFoWXnsn5Jkmabmk9ySdmZA+UtLi8NltSueIH+ecS6NMaf5K5umv\nA4FRZvaymf2J2EDIG1paoJm9Z2b5ZpYPjAS2Af8EpgJzzWwIMDe8R9Iw4FJiMyOfBdwhKT6h5Z3A\nNcCQsJ3V0no551x7VlZWRu/evencuXNa69FkUDGzG8zsc0ndJHUzsw/M7PSIyh8LrDCzD4DzgbtD\n+t3ABWH/fOABM9tuZquA5cAoSf2AHmb2ipkZcE9CHuec269kwhgVSK75a4SkRUApsFTS65KiWk/l\nUuD+sJ9nZh+F/TIgfg/XH0gc0b82pPUP+3XTnXNuv9NuggowA/iemR1qZgOB7wN/bm3BYazLecCs\nup+FOw9rbRkJZU2WtFDSwoqKiqhO65xzSSstLWXEiBGUlpam5PztKah0NbPd66eY2XygawRlnw28\nYWbl4X15aNIivG4I6euAQxLyDQhp68J+3fS9mNkMMys0s8KDDjoogqo751zyKisrGTduHEuXLmX8\n+PFUVlZGXkZ7CiorJf23pEFh+y9gZQRlX8aepi+A2cCVYf9K4LGE9EsldZY0mFiH/ILQVLZZ0ujw\n1NcVCXmccy5jTJw4kQ0bNmBmlJeXM2nSpEjPv3XrVrZu3dpugspE4CDgEeAfQJ+Q1mKSugKnh3PG\nTQdOl7QMOC28x8xKgYeApcDTwBQzi8+SfC3wF2Kd9yuAp1pTL+eci1pRURFz5syhqqoKgKqqKh5/\n/HGKiooiKyNTxqgAKNZ9Uc8H0r1m9p+SvmNmv2vjeqVMYWGhLVy4MN3VcM7tJ/Ly8nbPy5Wob9++\nu4NBa7344oucfPLJPPPMM5xxxhmRnLMuSa+bWWFTxzV2pzJS0sHAREm9w+DE3Vt0VXXOuX3XLbfc\nQqdOtefgzc3NZfr06ZGVkSkDH6HxoPJHYoMQjwJer7P5n/rOOZeEiRMncsABe/4Oz8nJ4dxzz+Xq\nq6ObQjGTmr8aDCpmdpuZHQ0UmdlhZjY4YTusDevonHPtVlVVFZs2bdr9Pi8vj5kzZ0ZaRllZGVlZ\nWfTp0yfS87ZEMiPqv9UWFXHOuX3RSy+9RFVVFWPHjgXgD3/4A127RjEqY4+ysjL69u1LdnZ20wen\nWFJr1DvnnGuZ4uJiOnbsyM033wywT49RgUaCiqT0zkrmnHP7gOLiYk488URGjRpFhw4dWLRoUeRl\ntIugArwMsUeL26guzjm3T6moqGDRokWcccYZdO7cmeHDh/Pmm29GXk4mBZXGFunqJOly4CRJF9X9\n0MweqSePc8654NlnnwXg9NNjE7vn5+fz9NNPR1qGmWVUUGnsTuWbwL8DvYBz62znpL5qzjnXvpWU\nlNC7d29GjhwJQEFBAeXl5bvHlUTh008/ZefOnRkTVBq8UzGzfwH/krTQzKJ9/s055/ZxZkZxcTFj\nx47d/VRWfn4+AIsWLeLss8+OpJxMWZs+Lpmnv+6VdL2kh8P2bUkdU14z55xrx959913WrVtXa9qU\neFCJsl8l04JKY30qcXcAHcMrwH8SW8b366mqlHPOtXfFxcXAnv4UgJ49e3LYYYdF+gRYewwqx5vZ\ncQnv50l6K1UVcs65fUFxcTFDhgxh0KBBtdLz8/NTcqeSCfN+QXLNX9WSDo+/kXQYUN3I8c45t1/b\nvn078+fPr3fG4IKCApYtW8aWLVsiKau8vJzOnTvTs2fPSM7XWskElR8Cz0maL+l5YB6xJYWdc87V\n4+WXX2bbtm21mr7i4v0qb70VTYNP/HHi2FqF6ddk85eZzZU0BDgyJL1nZttTWy3nnGu/SkpKyM7O\nZsyYMXt9VlBQAMQ6608++eRWl5VJY1QguT4VQhB5O8V1cc65fUJxcTGjR4+mR48ee3128MEH06dP\nn8g668vKyhg8eHAk54qCTyjpnHMR2rhxI6+//nqDKzBKoqCgILLO+ky7U2k0qCjmkKgLldQrjHl5\nV9I7kk4MK0qWSFoWXnsnHD9N0nJJ70k6MyF9pKTF4bPblCmNis65/dbcuXMxs0aX9c3Pz2fJkiXs\n3LmzVWXt2rWLioqK9hNULLaA/ZMpKPd3wNNmdhRwHPAOMBWYa2ZDiK04ORVA0jDgUmA4cBZwh6T4\nogF3AtcAQ8J2Vgrq6pxzSSsuLqZnz54UFja8nHtBQQE7duzgnXfeaVVZFRUVmFn7CSrBG5KOj6pA\nST2BU4CZAGa2w8w2AecDd4fD7gYuCPvnAw+Y2XYzWwUsB0ZJ6gf0MLNXQvC7JyGPc861OTOjpKSE\nsWPH0qFDw13W8c761varZNrAR0guqJwAvCxphaS3Q3NTazrtBwMVwF8lLZL0F0ldgTwz+ygcUwbE\nR/L0B9Yk5F8b0vqH/brpzjmXFu+//z4ffvhhvY8SJxoyZAi5ubmt7lfJxKCSzNNfZzZ9SLPL/ALw\nbTN7VdLvCE1dcWZmkiyqAiVNBiYDDBw4MKrTOudcLSUlJQCN9qcAZGdnc+yxx+6fdypm9gFwCPDl\nsL8tmXyNWAusNbNXw/uHiQWZ8tCkRXjdED5fF8qPGxDS1oX9uun1/QwzzKzQzAoPOuigVlTdOeca\nVlxczOGHH85hhx3W5LHx6Vpirfctk2lTtEASwUHSzcCNwLSQ1BG4r6UFmlkZsEZSfDDlWGApMBu4\nMqRdCTwW9mcDl0rqLGkwsQ75BaGpbLOk0eGprysS8jjnXJvauXMnzz33XJNNX3EFBQV89tlnrF69\nusVllpeX06NHD7p06dLic0QtmeavC4EC4A0AM1svqXsry/028DdJnYCVwNXEAtxDkiYBHwCXhPJK\nJT1ELPDsAqaYWXzusWuBu4AuwFNhc865NvfKK6+wdevWJpu+4hLXVmnp4MVMG6MCyQWVHYl9HKFT\nvVXM7E2gvuftxjZw/C+AX9STvhAY0dr6OOdcaxUXF5OVlVXv1Cz1OeaYY8jKyuLNN9/koov2WrE9\nKZkYVJLpG3lI0p+AXpKuAZ4F/pzaajnnXPtSUlLCCSecQK9evZI6vkuXLhx11FGt6qxvl0HFzP4v\nsc70fwBDgZ+Y2e2prphzzrUXn3zyCa+99lrS/SlxrZ2upV0GlWAx8L/AC2HfOedcMG/ePGpqapLu\nT4nLz89n7dq1fPzxx80u8/PPP+ezzz5rf0FF0teBBcBFwMXAK5ImprpizjnXXpSUlNCjRw9GjRrV\nrHyJ0+A3V3l5OZBZY1Qg+UW6CszsKjO7EhhJ7BFj55zb75kZxcXFjBkzho4dOzYrb+ITYM2ViQMf\nIbmgshFIXPdyS0hzzrn93ooVK1i9enWzm74ADjzwQA455JB9Kqg0+EixpO+F3eXAq5IeA4zYBI++\nYJdzzhF7lBhodid9XEs76zM1qDR2p9I9bCuAR4kFFIiNWl+V4no551y7UFJSwqBBgzjiiCNalD8/\nP5/33nuPbdu2NStfWVkZksi0qacavFMxs5+1ZUWcc6692blzJ/PmzWPChAm0dI3AgoICampqWLx4\nMSeccELS+crLy+nTp0+jU+ynQzJPfxVK+qekN8LU92+3cup755zbJyxYsIDNmze3qD8lrqWd9Zk4\nRgWSm6blb8SeAFsM1KS2Os45136UlJSQlZXFl7/85Raf49BDD6VXr17N7ldpz0Glwsxmp7wmzjnX\nzhQXF1NYWMgBBxzQ4nNIIj8/v0V3KkOGDGlxuamSzCPFN4fVGS+TdFF8S3nNnHMug23atIkFCxa0\nqukrrqCggLfffptdu3YldbyZtes7lauBo4itoxJv/jLgkVRVyjnnMt1zzz1HdXV1ix8lTlRQUEBV\nVRXvv/8+w4YNa/L4zZs3U1VV1W6DyvFmdmTThznn3P6juLiYbt26MXr06FafK7GzPpmgkqljVCC5\n5q+XJDX9Uzrn3H6kpKSEMWPG0KlTp1af66ijjqJz585Jd9a396AyGnhT0nvhceLF/kixc25/tnLl\nSlasWBFJ0xdAx44dGTFiRNKd9ZkcVJJp/jor5bVwzrl2pKSkBCCSTvq4goIC/vnPf2JmTQ6kzOSg\nksydijWwOefcfqm4uJhDDjmEoUOHRnbO/Px8Nm7cyNq1a5s8tqysjI4dO9K7d+/Iyo9KMkFlDvBE\neJ0LrASeak2hklaHZrQ3JS0MaQdIKpG0LLz2Tjh+mqTloQnuzIT0keE8yyXdppbOk+Ccc0natWsX\nc+fO5Ywzzmjx1Cz1ac7aKmVlZeTl5UVaflSSWU74GDM7NrwOAUYBL0dQ9hgzyzezwvB+KjA3lDE3\nvCc8JHApMJxYU9wdkrJDnjuBa4AhYfOmOudcSi1cuJDPPvsssv6UuGOPPRZJSfWrlJeXZ2TTFyS/\nnPBuZvYGkPysZ8k7H7g77N8NXJCQ/oCZbTezVcSm4h8lqR/Qw8xeMTMD7knI45xzKVFcXIwkxo4d\nG+l5u3XrxpAhQ5K+U8nUoNJkR33CuioQC0JfANa3slwDnpVUDfzJzGYAeWb2Ufi8DMgL+/2BVxLy\nrg1pO8N+3XTnnEuZkpISRo4cSZ8+fSI/d35+PgsWLGjyuLKyMkaOHBl5+VFI5k6le8LWmVjfyvmt\nLPdkM8sHzgamSDol8cNw5xHZwwCSJktaKGlhRUVFVKd1zu1nNm/ezMsvvxx501dcQUEBq1ev5tNP\nP23wmOrqajZs2NB+71RSsa6Kma0Lrxsk/ZNYP025pH5m9lFo2toQDl8HHJKQfUBIWxf266bXV94M\nYAZAYWGhP7nmnGuR+fPnU11dHemjxIninfVvvfUWp556ar3HbNy4kerq6owNKsmspzJU0gxJxZLm\nxbeWFiipq6Tu8X3gDGAJMBu4Mhx2JbEVJgnpl0rqLGkwsQ75BaGpbLOk0eGprysS8jjnXOSKi4vJ\nzc3lxBNPTMn5k1lbJZPHqEBygx9nAX8E/gJUR1BmHvDP8ChcB+DvZva0pNeAhyRNAj4ALgEws1JJ\nDwFLgV3AFDOL1+Na4C6gC7HHnFv1qLNzzjWmuLiYU089lc6dO6fk/Hl5efTr16/Rzvp9IajsMrM7\noyrQzFYCx9WTvhGo93EKM/sF8It60hcCI6Kqm3PONaS4uJhly5Zx4YUXprScptZWyfSgkkxH/eOS\nrpXULwxQPEBSy1ekcc65dqayspLLLrsMgPvuu4/KysqUlVVQUMA777xDVVVVvZ/Hg0peXl69n6db\nMkHlSmLLCb8EvB62hamslHPOZZKJEyeyadMmAD755BMmTZqUsrLy8/PZtWsXpaWl9X5eVlZGt27d\n6NatW8rq0BrJjKgfXM92WFtUzjnn0q2oqIgnnniCmprYGoVVVVU8/vjjFBUVpaS8pqZryeSBj9CC\nEfXOObc/mTZtGtu2bauVtm3bNqZNm5aS8g477DC6d+/eYL9KeXl5xjZ9gQcV55xr1NSpU/dKy83N\nZfr06SkpLysri+OOO87vVJxzbl9jZjz//PNkZWXtfow4JyeHc889l6uvvjpl5RYUFPDWW2/tbnJL\ntE8EFUnHSjpP0kXxLdUVc865dHvwwQd57LHH+PnPf757qvm8vDxmzpyZ0nLz8/PZunUry5cvr5W+\nfft2Pvnkk/YdVCQVAUXAV4Bzw3ZOiuvlnHNptWHDBq677jpGjRrFjTfeyJNPPsmwYcOYM2cOXbt2\nTWnZDXXWb9gQm70qk4NKMoMfR5vZsJTXxDnnMsh1113Hli1b+Otf/0p2djbDhw9nyZIlbVL2sGHD\n6NChA4sWLeKSSy7ZnZ7pAx8hueavl8NCWc45t1/4xz/+waxZs7j55psZNqztv/46d+7M8OHD97pT\naQ9BJZk7lXuIBZYyYDsgYrPTH5vSmjnnXBps3LiRa6+9loKCAn74wx+mrR75+fk8/fTTtdL2laAy\nE/hPYDGw96MIzjm3D7nhhhv45JNPKC4upmPHjmmrR0FBAXfffXetp73iQaVv375pq1dTkmn+qjCz\n2Wa2ysw+iG8pr5lzzrWxJ554gvvuu4+bbrqJ447ba97bNlXfNPhlZWUceOCBdOrUKV3ValIyQWWR\npL9LuswfKXbO7as2bdrEN77xDY455hh+/OMfp7s6u4NKYr9Kpo9RgeSav7oQ60tJXOrMgEdSUiPn\nnEuD73//+5SXlzN79uyMuBPo2bMngwcPrnWnUl5e3v6Dipmlbtioc85lgGeeeYaioiKmTp3KyJEj\n012d3QoKCvZq/jrhhBPSWKOmNRlUJP2V2J1JLWY2MSU1cs65NrR582auueYajj76aG6++eZ0V6eW\ngoICHnlEU3VvAAAYFUlEQVTkEbZs2UL37t33meavJxL2c4ALgfWpqY5zzrWtG2+8kbVr1/LSSy+R\nk5OT7urUEu9Xeeutt8jPz6eysjLjg0oy66n8I2H7G7G14wtbW7CkbEmLJD0R3h8gqUTSsvDaO+HY\naZKWS3pP0pkJ6SMlLQ6f3aaw8L1zziVj3rx5/PGPf+S73/0uo0ePTnd19pI4XUt7GKMCLZuleAgQ\nxUPS3wHeSXg/FZhrZkOAueE9YTT/pcBw4CzgDknZIc+dwDWhTkPC584516TKykq+/vWvc8QRR/Dz\nn/883dWp18EHH0yfPn1YtGjRvhNUJG2RtDn+CjwO3NiaQiUNAMYDf0lIPh+4O+zfDVyQkP6AmW03\ns1XAcmCUpH5ADzN7xcyM2Mj/C3DOuST8+Mc/ZtWqVcycOZPc3Nx0V6dekigoKNi37lTMrLuZ9Uh4\nHWpm/2hlubcCP6L2CP08M/so7JcB8aXN+gNrEo5bG9L6h/266c4516h//etf3H777Vx33XWccsop\n6a5Oo/Lz81myZAlr1sS+BjM9qCTTUY+kY4FBicebWYvGqUg6B9hgZq9LOrW+Y8zMJO31xFlLSZoM\nTAYYOHBgVKd1zrVDn3/+ORMnTuTQQw/llltuSXd1mlRQUMCOHTt47rnnyM7O5sADD0x3lRqVzCPF\nRcCxQCl77ixaM/jxi8B5ksYRe5qsh6T7gHJJ/czso9C0tSEcvw44JCH/gJC2LuzXTd+Lmc0AZgAU\nFhZGFqycc+3PT37yE5YtW8azzz5Lt27d0l2dJsWfAJs3bx55eXlkZWX2gr1tvp6KmU0DpgGEO5Uf\nmNl/SPo1cCUwPbw+FrLMBv4u6TfAwcQ65BeYWXXo6xkNvApcAdweVT2dc/uW0tJSzjvvPFatWsXk\nyZMZO3ZsuquUlKFDh5Kbm0tlZSVHHnlkuqvTpExaT2U6cLqkZcBp4T1mVgo8BCwFngammFl1yHMt\nsc7+5cAK4Kk2qKdzrp2prKzk7LPPZuXKlWRlZWXcIMfGZGdnc+yxsZVG3n//fUpLS9Nco8YlE1Ti\n66m8J+ntMC7k7SgKN7P5ZnZO2N9oZmPNbIiZnWZmnyQc9wszO9zMjjSzpxLSF5rZiPDZdeEpMOec\nq2XixIl89FHsOaDs7Gy+973vpblGzTN8+HAAtm7dyvjx46msrExzjRqWTFCJr6dyFnvWpz83lZVy\nzrmoFBUV8eijj7Jr1y4AduzYweOPP05RUVGaa5a8xJmKy8vLmTRpUhpr0zhfT8U5t8+qqKjgW9/6\nFjt27KiVvm3bNqZNm5amWjVPUVFRrSavqqqqjA6Kvp6Kc26fNHfuXI477jiqq6v3mso+NzeX6dOn\np6lmzTNt2jSqqqpqpWVyUEwmqCSup3Iue5rAnHMu4+zcuZOpU6dy+umn07NnT15//XUuuOCC3ZNF\n5uTkcO6553L11e1jVY9bbrmFrl271krL6KBoZvvVNnLkSHPO7ZuWL19uo0aNMsAmT55slZWVZma2\ndetWGzhwoEmyQw891LZu3ZrmmjbPJZdcYjk5OQZYTk6OTZgwoc3rACy0JL5jGxynIulHZvYrSbdT\n/3oq16cu1DnnXPP87W9/41vf+hbZ2dnMmjWLiy++ePdnXbt25cknn2TChAk8+OCDe/3ln+mKiooY\nNmwYa9asIS8vj5kzZ6a7Sg1qbPBjfAbhhW1REeeca4ktW7YwZcoU7r33Xk4++WT+9re/1Tsd0/Dh\nw1myZEkaath67SkoNhhUzOzxsPugmdXqJZLUJ6W1cs65JLz22mtcfvnlrFy5kp/+9KfcdNNNdOiQ\n1JSG7U57CYrJdNQvCFOhACDpK8BLqauSc87trbS0lBEjRlBaWkpNTQ2//vWvOemkk9i+fTvz58/n\n5ptv3mcDSnuSzL/A14AiSfOJzb11IPDlVFbKOecSVVZWMm7cONasWcNZZ53F0KFDmTdvHl/5ylf4\n85//TO/evZs+iWsTTQYVM1ss6RfAvcAW4BQzW9tENueci8zEiRPZsGEDZsbatWv56KOPmDFjBl//\n+tfxVcQzSzJT388EDic2/f1Q4AlJt5vZH1JdOeecKyoq4oknnqg1ALBjx45kZ2d7QMlAyfSpLAbG\nWGyalmeAE4AvpLZazrn9XU1NDXPnzuXaa69l27ZttT6rqqrK2BHl+7tklhO+NQx8ib//zMwydzYz\n51y7tmHDBn71q19x5JFHctppp5GVlbVXB3xGjyjfzzUZVCQNkfSwpKWSVsa3tqicc27/UFNTw7x5\n85gwYQIDBgzgxhtvpF+/ftx777188sknXHTRRe12mpX9TTLNX38F7gR2AWOIra9yXyor5Zzb9yQ+\nEhy3YcMGfv3rX3PkkUcyduxYSkpKmDJlCqWlpbzwwgv8x3/8Bzk5ORQVFdG3b18kZfyI8v1dUhNK\nmtlcQBab9v6nwPjUVss5ty+JPxK8dOlSxo8fz5w5c3bflfzoRz/afVeyfv16fvvb3zJsWO3FZuMj\nyocNG8acOXMyekT5/i6ZcSrbJWUByyRdB6wDuqW2Ws65fcnEiRMpLy/HzPjwww8555xz6N27N1Om\nTOGaa67ZK4jUp72MKN/fJXOn8h0gF7geGElsFcgrU1kp51zmqa/5qjFr167lwQcf5PTTT+fhhx9m\n+/btQGxm9E6dOnHLLbfUe1fi2jclPNjVNgVKOcALQGdid0oPm9nNkg4AHgQGAauBS8zs05BnGjAJ\nqAauD482I2kkcBexNV+eBL5jTfxAhYWFtnChz5HpXHNUVlbuniV34MCBlJaW1mqC2rVrF2+//TYv\nvvgiL730Ei+++CJr1qxp9Jx9+/alvLw81VV3EZH0upkVNnVcY1Pfz24so5md15KKEVvw68tmtlVS\nR+Bfkp4CLgLmmtl0SVOBqcCNkoYBlwLDiU0T86ykoWZWTewBgmuAV4kFlbOAp1pYL+dcAxJHtJeX\nl3PFFVdwzTXX7A4gr776KpWVlQD079+fL37xi/zgBz/gpJNOYtGiRXz3u9/d/Tn4I8H7tIYWWgEq\ngDeAHwKnAF9K3JJZrKWpjViz2hvEBlS+B/QL6f2A98L+NGBaQp5ngBPDMe8mpF8G/KmpMn2RLre/\nWrJkiQ0fPtyWLFmSdJ6amhq79dZbdy8QVXfLysqyL3zhC3bdddfZ/fffbx988EG958mERaZc69Da\nRbqAfwNOD1/WlwNzgPvNLLkG1UZIygZeB44A/mBmr0rKM7OPwiFlQF7Y7w+8kpB9bUjbGfbrpju3\nzyktLd29lsbw4cObnT9xQsbx48fXar6qrq5mzZo1rFixot5ty5Yt9Z6zV69erFmzhm7dmn5upz0t\nMuVap8GOejOrNrOnzexKYDSwHJgfngBrlXDufGAAMErSiDqfx/8SioSkyZIWSlpYUVER1WndfqS5\nndRR5q/7OG5iM1Iyampq+NrXvrb76at169aRn5/P2WefzdChQ+nSpQuDBw/mtNNO4xvf+Aa33nor\nS5Ys4eCDD+aqq67i0ksvpXPnzrXOmZuby29+85ukAgr4I8H7lcZuY4h1pl8EzAJeA/4b6J/MLVCy\nG/AT4AdkePNXS5oOPP++kb+165u3Nn/dpqNLLrnEPv74Y1u6dKk9//zzNmvWLLvjjjvsZz/7mU2Z\nMsUuueQSO/XUU2348OHWt2/feputADvkkEPs4osvthtvvNFmzJhhc+fOtdWrV9uuXbuarIM3X+1/\nSLL5q7Ev+3uI9Xf8DzAimZMlVSAcBPQK+12A/wXOAX4NTA3pU4Ffhf3hwFshwA0GVgLZ4bMFxO6i\nRKyDflxT5bckqKT7S8XzZ9aXenO/UOvmP++88+zdd9+1BQsW2Ny5c+3RRx+1e++91+644w775S9/\naTfddJNdf/31dvXVV9vIkSMtOzu7wcBQd+vdu7cNHTrUTj75ZLvwwgtt8uTJlpubW++xffv2bbNr\n6Nq/KIJKDbH1U7YAmxO2LcDmZE7ewHmPBRYBbwNLgJ+E9AOBucAy4FnggIQ8NwErwt3M2QnpheEc\nK4DfEx6RbmxrSVCJ+kvF87c+f01Nje3atcu2b99u27Zts82bN9unn35qH3/8sZWXl9v69ettzZo1\ntmrVKhs3bpx17tzZAOvcubOdfvrp9tJLL9nzzz9vzz77rD311FP22GOP2cMPP2z333+/3XPPPTZz\n5ky788477bLLLrNOnTrV+jLu0KGDjRkzxiZOnGiXX365XXTRRTZu3DgbM2aMnXjiiVZQUGBHH320\nDR482Hr06JF0QIhvWVlZ1rNnTxswYECDAaVbt27297//3UpKSuytt96y9evX244dO+q9fjNnzrSu\nXbvWyp+bm2tFRUXN+ndo7d2ia9+SDSptPk4l3Zo7TqWoqIjrr7++Vjt2dnY2xx9/PEcccQT1Xb/E\ntBUrVrBw4UKqq6tr5R85ciSHHXZYrWPj+4lpK1eu5M0339wr/3HHHcehhx5ab77E1w8//JAlS5bs\nlf/oo49mwIAB9QX9Wu/Xr1/PsmXLqKmp2Z0/KyuLwYMH07dvX8yMmpqaWnkS32/cuJF169bV+pkk\n0adPH7p3705NTU2j27Zt22qto5EpJHHwwQfTpUsXcnJy9tri6bNmzaq3/j169OC+++6je/fudO/e\nnR49euzez83N3b1OSH2/f7m5ufz+979v1oSKEyZMYPbs2VRVVZGTk8P555/PAw880PoL4fYbyY5T\n8aDShLy8PDZs2LBXelZWFoMGDQKod6GgeNqqVatqfaHHZWdnc/jhh++VP74ff33//fcbzH/00Ufv\ndXzd1yVLlrBr16698nfo0IGCggIk7bXF80vi1VdfZefOnXvl79ixI1/60pfIysqqlbfu+6effpod\nO3bslb9z58589atfJSsrq9Ft5syZfP7553vl79q1Kz/60Y/Izs4mOzubDh067N5PTPve975X79NL\nvXr14uGHH6Zjx4507NiRTp067d5PTJs1axbTpk2rtZ5Hc77UowgKUQSEpgYvOteUZINKJP0k7Wlr\nbvNXa5sOPH/7zm+W/ua/qPozvPnKtQat7VPZVzfvU/H8bf2lHkVQ8IDg0s2DSoRBJd1fKp6//X+p\ne1Bw7Z0HlQiDiln6v1Q8v3+pO5dOyQYV76h3zjnXpGQ76pNZT8U555xLigcV55xzkfGg4pxzLjIe\nVJxzzkVmv+uol1QBfJCi0/cEPkthvqaOa+zz+j5LJi3xfR/g4yTq2VItvX7Nydvca9TYZ5l2/Rqq\nU5T50nn9IHN/B1P9fzgTrt+hZnZQk0cl84iYb8ltwIxU5mvquMY+r++zZNIS35PkI4Vtff2iuobt\n/fq15hq2h+vXFtcwU/8Pt5frZ2be/BWxx1Ocr6njGvu8vs+SSWvpz9QSrSkrimvY3q9fa8rz69e6\n8lL9f7i9XL/9r/nLtZykhZbMhHKuXn79Ws+vYeu0xfXzOxXXHDPSXYF2zq9f6/k1bJ2UXz+/U3HO\nORcZv1NxzjkXGQ8qzjnnIuNBxTnnXGQ8qLhISLpA0p8lPSjpjHTXp72RdJikmZIeTndd2gtJXSXd\nHX7vvpbu+rQ3qfqd86DikFQkaYOkJXXSz5L0nqTlkqY2dg4ze9TMrgG+CUxIZX0zTUTXb6WZTUpt\nTTNfM6/lRcDD4ffuvDavbAZqzvVL1e+cBxUHcBdwVmKCpGzgD8DZwDDgMknDJB0j6Yk6W9+ErP8V\n8u1P7iK667e/u4skryUwAFgTDqtuwzpmsrtI/vqlRIdUndi1H2b2gqRBdZJHAcvNbCWApAeA883s\nFuCcuueQJGA68JSZvZHaGmeWKK6fi2nOtQTWEgssb+J/IAPNvn5LU1EH/4dwDenPnr8CIfYfuH8j\nx38bOA24WNI3U1mxdqJZ10/SgZL+CBRImpbqyrUzDV3LR4CvSLqTNExH0o7Ue/1S9TvndyouEmZ2\nG3BbuuvRXpnZRmL9US5JZlYJXJ3uerRXqfqd8zsV15B1wCEJ7weENJccv37R8WvZOm16/TyouIa8\nBgyRNFhSJ+BSYHaa69Se+PWLjl/L1mnT6+dBxSHpfuBl4EhJayVNMrNdwHXAM8A7wENmVprOemYq\nv37R8WvZOplw/XxCSeecc5HxOxXnnHOR8aDinHMuMh5UnHPORcaDinPOuch4UHHOORcZDyrOOeci\n40HFtZqkmySVSnpb0puSTgjpN0jKTWG5B0exFoSkr0p6R9JzrThHdfjZ41ujU91LukrS71taXmtJ\nOlXSSQnv75J0cSvO12B+SbdKOqWBOjzR0jLrOd9Bkp6O6nyuZXzuL9cqkk4kNuvuF8xsu6Q+QKfw\n8Q3AfcC2VJRtZuuBFn8RJpgEXGNm/0rmYEkdwoCyRJ+bWX4EdWlOma1xKrAVeCnCc+5F0oHAaDO7\nIZXlAJhZhaSPJH3RzF5MdXmufn6n4lqrH/CxmW0HMLOPzWy9pOuBg4Hn4ncAks6Q9LKkNyTNktQt\npK+W9CtJiyUtkHRE3UIkfSnhLmCRpO6SBsUXI5L0l4TPKyTdHNJ/KOm1cBf1s3rO+xPgZGCmpF9L\nypH011CXRZLGhOOukjRb0jxgbrIXJ/xsfcJ+oaT59RxzkKR/hHq+JumLIf2nku6V9CJwb508p0p6\nXtJjklZKmi7pa+H6LZZ0eEPnVmxq9G8C3w3X69/DaU+R9FI438Uhv8J1WRLOOyEh/feKLfz0LNDQ\nmjBfAXbfPSi2WNS7kt4gtshWPH1U+N1YFOpwZEh/QVJ+wnH/knRcfb8P4ZBHAV8FMp3MzDffWrwB\n3YitZ/E+cAfwpYTPVgN9wn4f4AWga3h/I/CThONuCvtXAE/UU87jwBcTyuwADAKW1DnuUGJTURwK\nnAHMAETsD6gngFPqOfd8oDDsfx8oCvtHAR8COcBVxKYMP6CB61AdrkN8m1DPNSgE5of9q4Dfh/2/\nAyeH/YHAO2H/p8DrQJd6yjsV2EQsqHcmNkHgz8Jn3wFuTeLcP0g4313ArHCdhhFbfwNiQaEEyAby\nwvXoRywgxNMPDnW5uJ563g2cG/ZziE3BPiT8mzwU/7cGegAdwv5pwD/C/pUJP8tQYGFDvw9hvz+w\nON3/L/bnzZu/XKuY2VZJI4F/B8YAD0qaamZ31Tl0NLEvqxclQayJ7OWEz+9PeP1tPUW9CPxG0t+A\nR8xsbTjPbpJyiH0xftvMPpD0bWKBZVE4pBuxL7QXGvmRTgZuDz/bu5I+IPZlBlBiZp80kK81zV+n\nAcMSfp4e8bs4YLaZfd5AvtfM7CMASSuA4pC+mNi/RVPnrutRM6sBlkrKC2knA/ebWTVQLul54Hjg\nlIT09eEOrj79gIqwfxSwysyWhTrfB0wOn/UE7pY0BDCgY0ifBfy3pB8CE4kFP6jn9yGkbyAW5Fya\neFBxrRa+WOYD8yUtJvbX5V11DhOxL+XLGjpNA/vxMqZLmgOMIxaYzgSq6hz2R2JfMM8mlHmLmf2p\nGT9OYypbkGcXe5qZcxo4JotYv0OtnycEgsbK3J6wX5PwvoY9/7cbO3dj56v3gBb4nIZ/7kQ/B54z\nswtD89x8ADPbJqmE2EqFlwAjQ/pevw9m9m4oq6Eg7NqA96m4VpF0ZPjrMi4f+CDsbwHibd2vAF+M\n95dI6ippaEK+CQmviXcw8XION7PFZvZLYlN5H1Xn8ylAdzObnpD8DDBRe/pu+qvp9eD/l9AmH+o3\nEHiviTyNWU34IiTWlFSfYmIrZxLKjbLDv6FzJ/7bNOZ/gQmSsiUdROwOZQGxu714ej/23BnV9Q4Q\n7yN7FxgU7+8BEv/A6MmeNT6uqnOOvxBbAO41M/s0/BwN/T4MBZYk8XO5FPGg4lqrG7Fmi6WS3ibW\nxPXT8NkM4GlJz5lZBbEvi/vDcS9TOzD0DunfAb5bTzk3hM7it4GdwFN1Pv8BcExC5+03zayYWJ/C\ny+EO6mGa/iK9A8gKxz8IXGXhIYQmdFHtR4rjwe1nwO8kLSTW71Kf64FCxR4mWEq0q/E1dO7HgQvr\ndNTX55/A28BbwDzgR2ZWFtKXEVvn/B7q+UMgmEOs/4dwtzQZmBM66jckHPcr4BZJi6jTgmJmrwOb\ngb8mJDf0+zAmlOnSxKe+d2knaTWxjvKP010XFz1J/wLOMbNNLcx/MLHmsKNCn09jx74AnB+/o3Ft\nz+9UnHOp9n1izYjNJukK4FViTwc2FVAOAn7jASW9/E7FOedcZPxOxTnnXGQ8qDjnnIuMBxXnnHOR\n8aDinHMuMh5UnHPORcaDinPOucj8f5F93EjhljuZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e0650b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogx(step_sizes, maximums, 'kd-')\n",
    "plt.ylabel(\"Maximum number of foxes in 600 days\")\n",
    "plt.xlabel(\"Step size for Euler method (days)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smaller step sizes give more accurate results, and we see it converges towards the \"true\" solution as we go left on the graph. Let's assume the smallest we tried is in fact \"true\", and calculate the errors relative to this. Then plot those against step size on a log-log plot. The slope should give us the order of convergence of the method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "maximums = np.array(maximums) # convert list into a numpy array\n",
    "errors = maximums - maximums[-1] # calculate error relative to the last value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEQCAYAAACqduMIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X94VOWd9/H3N0DUwAJFJcsjEvQBpVDrSrh07Q8LSHdp\nnyISkR+NiBRItVW2ax9BpdTWlta1ovUHW4pAZdsAtRRtUKSsCvaCsOuPtio+iFIXEFgXAfkRIIEk\n3+cPEjuESTKTmck5J/N5Xddc5twz3PPJXGO+1zn3ue/b3B0REZGGcoIOICIi4aQCISIicalAiIhI\nXCoQIiISlwqEiIjEpQIhIiJxqUCIiEhcKhAiIhKXCoSIiMQVqgJhZh3N7FUz+0rQWUREsl37THZu\nZouArwB73P1TMe3DgYeBdsACd7+v7qkZwJOJ9n/OOed47969U8p45MgROnbsmFIfQVH2YCh7MJQ9\nfV577bW97n5usy9094w9gKuAgcCmmLZ2wF+AC4Fc4HWgP/BFYBxwE/CVRPovLCz0VK1duzblPoKi\n7MFQ9mAoe/oAr3oCf2Mzegbh7n8ws94Nmi8Htrr7ewBmtgwYCXQCOtYVi2NmtsrdazOZT0REGmee\n4dVc6wrEM153icnMRgPD3X1K3fEE4Ap3v7Xu+CZgr7s/00h/JUAJQH5+fuGyZctSyldRUUGnTp1S\n6iMoyh4MZQ+GsqfPkCFDXnP3Qc29LqNnEC3h7k808/x8YD7AoEGDfPDgwSm937p160i1j6AoezCU\nPRjK3vqCuItpF3B+zHHPuraEmdkIM5t/8ODBtAYTEZG/CqJAvAL0NbMLzCyXkwPTZcl04O4r3b2k\nS5cuGQkoIuFRWlpK7969GTp0KL1796a0tDToSFkjowXCzJYCG4GLzWynmU1292rgVuD3wGbgSXd/\nK8l+dQYhkgVKS0spKSlh+/btuDvbt2+npKRERaKVZLRAuPt4d+/h7h3cvae7L6xrX+XuF7n7/3b3\n2S3oV2cQIllg5syZHD169JS2o0ePMnPmzIASZZdQzaQWEYm1Y8eOpNolvSJZIHSJSSQ79OrVK6l2\nSa9IFghdYhLJDrNnzyYvL++Utry8PGbPTvrKtLRAJAuEiGSH4uJi7r777o+PCwoKmD9/PsXFxQGm\nyh6hmyiXCDMbAYzo06dP0FFEJMPqrxT8+te/ZsyYMQGnyS6RPIPQJSaR7LFhwwZ69uxJ9+7dg46S\ndSJZIEQke5SXl/OZz3wm6BhZSQVCREJr586d7Nixg89+9rNBR8lKkSwQus1VJDuUl5cD6AwiIJEs\nEBqDEAm3+vWTcnJyUlo/qby8nLy8PC699NI0J5RERPIuJhEJr/r1k+qXyKhfPwlI+vbU8vJyLr/8\ncjp06JD2nNK8SJ5BiEh4pWv9pKNHj/KnP/1Jl5cCFMkCoTEIkfBK1/pJr7zyCtXV1SoQAYpkgdAY\nhEh4pWv9pPoB6iuvvDLlTNIykSwQIhJeP/zhD8nJOfVPS0vWT9qwYQOf/OQn6datWzrjSRJUIEQk\nrY4fP05tbS1nn332x2133HFHUgPUtbW1bNy4UZeXAqYCISJps2/fPqZPn87nPvc5PvzwQw4cOECH\nDh04cuRIUv2888477N+/XwUiYCoQIpI2d999NwcOHOBf//VfMTO6dOnC1VdfzYoVK3D3hPupH3/Q\nDOpgRbJA6C4mkfB5+eWXefzxx5k2bRqXXHLJx+1FRUW89957vPHGGwn3tWHDBrp168ZFF12UiaiS\noEgWCN3FJBIuNTU13HLLLfzt3/4t3/ve9055buTIkZgZK1asSLi/+gX6zCzNSSUZkSwQIhIuP//5\nz/njH//Igw8+SOfOnU95rnv37nz+85/nqaeeSqivffv28fbbb+vyUgioQIhISvbs2cPMmTMZOnQo\nY8eOjfuaUaNG8eabb/Luu+8229/GjRsBLdAXBioQIpKS6dOnc+TIEebOndvoJaFRo0YBJHQWUV5e\nTvv27Rk0aFBac0ryVCBEpMXWr1/P4sWL+fa3v02/fv0afV1BQQGFhYUJjUOUl5dz2WWXkZeXl86o\n0gIqECLSItXV1XzjG9+gV69efOc732n29UVFRfznf/4nu3btavQ1J06c4OWXX9b4Q0ioQIhIUur3\neujQoQNvvvkmRUVFdOzYsdl/V1RUBMDTTz/d6Gv+/Oc/c+zYMY0/hIQKhIgkrH6vh+3bt3/cNn/+\n/IQ2BOrXrx/9+vVr8jKTdpALl0gWCE2UEwlGqns9FBUV8dJLL7F37964z5eXl1NQUMB5552XclZJ\nXSQLhCbKiQQj1b0eioqKqKmpYeXKlac95+5s2LBBZw8hEskCISLB6NGjR9z2RPd6GDhwIL169Yp7\nmen9999n165dKhAhogIhIgmpqqqKuzd0Mns9mBlFRUX8+7//O4cPHz7lOS3QFz4qECKSkLvvvpvt\n27dz++23U1BQgJlRUFDA/Pnzk9rroaioiKqqKp577rlT2jds2EDHjh1PWehPgqUCISLNWr16NQ8+\n+CDf/OY3mTNnDtu2baO2tpZt27YlVRzg5B1K55577mmXmcrLy7niiito3759OqNLClQgRKRJH3zw\nARMnTuSSSy7hJz/5Scr9tWvXjmuvvZZnn32WyspKACoqKnj99dc1/hAyKhAi0qja2lomTpzIoUOH\nWLZsGWeddVZa+i0qKqKiooIXXngBOLmXRE1NjcYfQkYFQkQa9dBDD7FmzRp++tOf0r9//7T1O3To\nUDp37vzxZab6Aeq///u/T9t7SOpUIEQkrtdee4277rqLoqIiSkpK0tp3bm4uI0aM4He/+x3V1dWU\nl5czYMAAunbtmtb3kdSEpkCY2SfNbJ6ZLTezW4LOI5LNDh8+zLhx48jPz+fxxx/PyM5uo0aNYt++\nfbz00kts3LhRl5dCKKMFwswWmdkeM9vUoH24mW0xs61mdieAu29295uBMYC+KSIBqF+Ir3Pnzmzd\nupWbbrqJbt26ZeS9hg8fzplnnsns2bM5cOCABqhDKNNnEE8Aw2MbzKwdMBf4EtAfGG9m/eueuwZ4\nFliV4Vwi0kC8hfgefPDBhBbia4mOHTsyYMAA1q5dC5ycZ5Gp95KWyWiBcPc/APsbNF8ObHX399z9\nOLAMGFn3+jJ3/xKQ3I3VIpKyVBfiS1ZpaSlvvvnmx8e7d++mpKRERSJEzN0z+wZmvYFn3P1Tdcej\ngeHuPqXueAJwBbAcKALOAN5w97mN9FcClADk5+cXLlu2LKV8FRUVdOrUKaU+gqLswWir2YcOHUq8\nvwdmxosvvpj2LOPGjeN//ud/TmvPz88n3v/XbfVzD8KQIUNec/fm93R194w+gN7Appjj0cCCmOMJ\nwGMt6buwsNBTtXbt2pT7CIqyB6MtZj9+/LifeeaZDpz2KCgoyEgWM4v7fmaWVPYoCFt24FVP4G9s\nEHcx7QLOjznuWdeWMO0HIZI+7s60adOorKwkNzf3lOeSWYgvWY2tAJvoyrCSeUEUiFeAvmZ2gZnl\nAuOAsmQ6cO0HIZI2c+fOZd68ecyYMYNFixaltBBfMmbPnk1eXt4pbZksSJK8jK6KZWZLgcHAOWa2\nE7jH3Rea2a3A74F2wCJ3fyvJfkcAI/r06ZPuyCJZZc2aNfzTP/0TI0eO5Ec/+hE5OTkZKwgN1b/P\nzJkz2bFjB7169WL27Nmt9v7SvIwWCHcf30j7KlK4ldXdVwIrBw0aNLWlfYhku82bNzNmzBguueQS\nfvWrX5GT0/oXFIqLi1UQQiw0M6lFpPXs27ePESNGcMYZZ1BWVhaqO2wkPCK58LouMYm03PHjx7nu\nuuvYuXMna9eu1aCwNCqSZxAapBZJTv0SGkOHDuXss8/mpZdeYuHChVx55ZVBR5MQi+QZhIgkrn4J\njfpZ0hUVFdq1TRISyTMIzYMQSVy8JTSqq6sztoSGtB2RLBC6xCSSuB07diTVLlIvkgVCRBLXo0eP\nuO0anJbmqECItGF79+6ltrb2tHbNWJZERLJAaAxCpHnHjh3jmmuu4aOPPmLWrFmttoSGtB2RLBAa\ngxBpWk1NDcXFxfzHf/wHpaWl3HvvvWzbto0XX3yRbdu2qThIQnSvm0gb4+5861vf4qmnnuLhhx/m\nuuuuCzqSRFQkzyBEpHEPPPAAjz32GN/+9reZNm1a0HEkwlQgRNqQJUuWMH36dMaOHcv9998fdByJ\nuEgWCA1Si/xV/TIa9Ut19+vXj8WLFweyOqu0LZH8BmmQWuSk+mU0tm/f/vF+0tu3b2f58uUBJ5O2\nIJIFQkROireMxrFjx7SMhqSFCoRIhGkZDckkFQiRiDp48GCjq7JqGQ1Jh0gWCA1SS7Y7evQoI0aM\noKamhjPOOOOU57SMhqRLJAuEBqklm504cYLrr7+e9evXs2TJEhYuXKhlNCQjNJNaJEJqa2u56aab\nWLVqFfPmzWPs2LEAKgiSEZE8gxDJRu7OtGnTWLJkCT/+8Y/5+te/HnQkaeNUIEQi4p577mHu3Lnc\ncccdzJgxI+g4kgVUIERCKnaGdLdu3fjBD37A5MmT+Zd/+RfMLOh4kgU0BiESQvUzpOsnwX300Ue0\na9eOL3zhCyoO0mp0BiESQvFmSNfU1DBr1qyAEkk2UoEQCSHNkJYwiGSB0EQ5aeu6d+8et10zpKU1\nRbJAaKKctGUvv/wyBw4cOG2sQTOkpbVFskCItFVvvPEGw4cP57zzzuORRx7RDGkJlO5iEgmJLVu2\n8MUvfpG8vDxeeOEFevfuza233hp0LMliOoMQCYFt27YxbNgwgI+Lg0jQdAYhErDdu3dz9dVXc+TI\nEdatW8fFF18cdCQRQGcQIoGInSVdUFDArl27WL16NZ/+9KeDjibyMRUIkVbWcB/p6upqAN59992A\nk4mcSgVCpJXFmyVdVVWlfaQldFQgRFqZZklLVISqQJjZtWb2uJn92sz+Ieg8IulWWVl52hah9TRL\nWsIm4wXCzBaZ2R4z29SgfbiZbTGzrWZ2J4C7P+3uU4GbgbGZzibSmo4fP86YMWOorKwkNzf3lOc0\nS1rCqDXOIJ4Ahsc2mFk7YC7wJaA/MN7M+se85Dt1z4u0CdXV1RQXF7Ny5Urmzp3LokWLNEtaQi/j\n8yDc/Q9m1rtB8+XAVnd/D8DMlgEjzWwzcB/wnLv/MdPZRFpDbW0tkyZNYvny5cyZM4dvfOMbgPaR\nlvAzd8/8m5wsEM+4+6fqjkcDw919St3xBOAK4B1gIvAK8Gd3nxenrxKgBCA/P79w2bJlKWWrqKig\nU6dOKfURFGUPRjLZ3Z05c+bw7LPPMnnyZG644YYMp2tatnzuYRO27EOGDHnN3Qc1+0J3b/IBtAMe\naO51zfTRG9gUczwaWBBzPAF4LNl+CwsLPVVr165NuY+gKHswEs1eW1vrt912mwM+c+bMzIZKUDZ8\n7mEUtuzAq57A39hmxyDcvQb4XIvKVON2AefHHPesa0uI9oOQsIqdId21a1ceffRRbr/9dn7wgx8E\nHU0kaYmOQfzJzMqA3wBH6hvdfUUL3/cVoK+ZXcDJwjAO+Gqi/9jdVwIrBw0aNLWF7y+Sdg33kT50\n6BDt27dn4MCB2kdaIinRu5jOBPYBQ4ERdY+vJPIPzWwpsBG42Mx2mtlkd68GbgV+D2wGnnT3t5IN\nLxIm8WZIV1dXa4a0RFZCZxDuPqmlb+Du4xtpXwWsakmfZjYCGNGnT5+WxhJJO82QlrYmoTMIM+tp\nZk/VTXjbY2a/NbOemQ7XGNeWoxJCn/jEJ+K2a4a0RFWil5h+AZQB/6vusbKuLRAapJawWbRoEfv3\n76ddu3antGuGtERZogXiXHf/hbtX1z2eAM7NYK4m6QxCwqS0tJQpU6bwj//4jyxcuFAzpKXNSPQu\npn1mdgOwtO54PCcHrUWy2m9+8xtuvPFGBg8ezIoVK8jLy2PixIlBxxJJi0TPIL4GjAE+AP6bkxPd\nWjxwnSpdYpIwKCsr46tf/SpXXnklZWVl5OXlBR1JJK2aPYOoW1ivyN2vaYU8CdE8CAlCaWkpM2fO\nZMeOHZx77rns27ePwsJCVq1aFaplFETSJdGZ1HFvVRXJFg23Cd2zZw+1tbVMnjyZzp07Bx1PJCMS\nvcS0wcweM7PPm9nA+kdGk4mESLxJcO7Oj370o4ASiWReooPUf1f333tj2pyTM6tbnSbKSWvTJDjJ\nRs2eQZhZDvAzdx/S4BFIcQDd5iqtLz8/P267JsFJW5bIGEQtML0VsoiE0htvvMHhw4dPW3BPk+Ck\nrUt0DOJ5M/u/Zna+mXWrf2Q0mUgIbN68mWHDhtG1a1fmzJmjSXCSVRIdgxhb999vxrQ5cGF644iE\nx7vvvsvVV19NTk4OL774IhdddBH//M//zLp16xg8eHDQ8UQyLtHVXC/IdJBkaJBaMu2//uu/GDp0\nKCdOnGDdunVcdNFFQUcSaXVNXmIys+kxP1/f4LnA7u/TILVk0vvvv8/QoUM5cuQIzz//PAMGDAg6\nkkggmhuDGBfz810Nnhue5iwigYndKvTCCy/kgw8+YM2aNVx66aVBRxMJTHMFwhr5Od6xSCQ1nCVd\nXV2Nu7Nly5ago4kEqrkC4Y38HO9YJJLizZKuqqrSVqGS9ZobpL7UzA5x8mzhrLqfqTs+M6PJRFqJ\nZkmLxNdkgXD3dk09HxTdxSTpUlFRQW5uLlVVVac9p1nSku0SnSgXKrqLSdLh2LFjXHPNNVRVVZGb\nm3vKc5olLRLRAiGSqqqqKoqKili3bh2//OUvWbRokWZJizSQ6ExqkTbjxIkTjB8/ntWrVzN//nxu\nuOEGABUEkQZ0BiFZpaamhokTJ/LUU0/x8MMPM3WqNiUUaYwKhGSN2tpapk6dytKlS7nvvvuYNm1a\n0JFEQk0FQtq02BnSXbp04Re/+AXf/e53mTFjRtDRREJPYxDSZtXPkK6fBFdRUUH79u218J5IgiJ5\nBmFmI8xs/sGDB4OOIiEWb4Z0dXW1ZkiLJCiSBULzICQRmiEtkppIFgiRRHziE5+I264Z0iKJUYGQ\nNmnp0qXs37+fnJxTv+KaIS2SOBUIaXPKysqYMGECX/jCF1iwYIFmSIu0kO5ikjblhRdeYMyYMQwc\nOJCysjI6d+7MpEmTgo4lEkk6g5A2Y+PGjYwcOZK+ffvy3HPP0blz56AjiUSaCoS0Ca+//jpf/vKX\n6dGjB2vWrOHss88OOpJI5KlASGTFzpIeOHAgZsbzzz9Pjx49go4m0iaoQEgkNdxHura2lmPHjrF+\n/fqgo4m0GaEpEGZ2oZktNLPlQWeR8Is3S7qyslKzpEXSKKMFwswWmdkeM9vUoH24mW0xs61mdieA\nu7/n7pMzmUfaDs2SFsm8TJ9BPAEMj20ws3bAXOBLQH9gvJn1z3AOaUMqKytP2yK0nmZJi6RPRguE\nu/8B2N+g+XJga90Zw3FgGTAykzmk7aiurmbcuHHaR1qkFZi7Z/YNzHoDz7j7p+qORwPD3X1K3fEE\n4ArgHmA28EVggbv/uJH+SoASgPz8/MJly5allK+iooJOnTql1EdQsi27u/OTn/yE5557jmnTpvE3\nf/M3LFiwgD179tC9e3emTJnCsGHDMpT4r7Ltcw8LZU+fIUOGvObug5p9obtn9AH0BjbFHI/mZAGo\nP54APNaSvgsLCz1Va9euTbmPoGRb9unTpzvg3/3ud9MfKAnZ9rmHhbKnD/CqJ/A3Noi7mHYB58cc\n96xrS5j2g8g+DzzwAPfffz+33HIL3/ve94KOI5IVgigQrwB9zewCM8sFxgFlyXTg2g8iqyxevJg7\n7riDMWPG8Oijj2JmQUcSyQqZvs11KbARuNjMdprZZHevBm4Ffg9sBp5097cymUOiJXaGdPfu3Zk0\naRLDhg3j3/7t32jXrl3Q8USyRkZXc3X38Y20rwJWtbRfMxsBjOjTp09Lu5CQariP9IcffkhOTg7j\nx4/njDPOCDidSHYJzUzqZOgSU9sVb4Z0bW0t9957b0CJRLJXJAuEBqnbLs2QFgmPSBYInUG0XT17\n9ozbrhnSIq0vkgVC2qYTJ04Qr+hrhrRIMFQgJBTcnZtvvplNmzYxdepU7SMtEgKR3JNadzG1PbNn\nz2bRokXMmjVLA9IiIRHJMwiNQbQtv/rVr5g1axYTJkzg+9//ftBxRKROJAuEtB1r167la1/7GkOG\nDGHBggWaJS0SIpEsELrNtW146623GDVqFH379mXFihWN7vEgIsGIZIHQJaboql9GY+jQoVx66aW4\nO6tWraJr165BRxORBiJZICSa6pfR2L59O+5OTU0Nx48fZ/369UFHE5E4VCCk1cRbRqOyspKZM2cG\nlEhEmqICIa1Gy2iIREskC4QGqaPpnHPOiduuZTREwimSBUKD1NGzefNmDh06RE7OqV85LaMhEl6R\nLBASLQcPHmTUqFF06dKFhx56SMtoiEREJJfakOiora3lxhtv5C9/+QsvvPACV111FdOmTWPdunUM\nHjw46Hgi0gQVCMmoH/7wh5SVlfHII49w1VVXBR1HRJKgS0ySMc888wz33HMPN954I7feemvQcUQk\nSZEsELqLKfzeeecdiouLGThwIPPmzdMaSyIRFMkCobuYwu3w4cNce+215ObmsmLFCs4666ygI4lI\nC0SyQEj41K+xlJOTQ35+Pm+//TZPPvkkBQUFQUcTkRbSILWkrH6NpfplNI4dO0aHDh3YvXt3wMlE\nJBU6g5CUxVtj6cSJE1pjSSTiVCAkZVpjSaRtUoGQlDW2lpLWWBKJNhUISdno0aNPa9MaSyLRF8kC\noXkQ4bF//35KS0vp2bMnvXr10hpLIm1IJO9icveVwMpBgwZNDTpLtrvtttvYu3cvL7/8MpdddlnQ\ncUQkjSJ5BiHh8Nvf/pYlS5Ywa9YsFQeRNkgFQlpkz5493HzzzRQWFnLXXXcFHUdEMkAFQpLm7nz9\n61/n0KFDLF68mA4dOgQdSUQyIJJjEBKs0tJSnn76ae6//34GDBgQdBwRyRCdQUhSdu3axW233cZn\nPvMZbr/99qDjiEgGqUBIwtydKVOmUFVVxRNPPEG7du2CjiQiGaQCIc2KXal19erVXH/99fTt2zfo\nWCKSYSoQ0qT6lVq3b9/+cdvy5cspLS0NMJWItIbQFAgz62hmi83scTPTFNyQiLdS69GjR7VSq0gW\nyGiBMLNFZrbHzDY1aB9uZlvMbKuZ3VnXXAQsd/epwDWZzCWJ00qtItkr02cQTwDDYxvMrB0wF/gS\n0B8Yb2b9gZ7A+3Uvq8lwLkmQVmoVyV4ZLRDu/gdgf4Pmy4Gt7v6eux8HlgEjgZ2cLBIZzyWJu+WW\nW05r00qtItkhiIly5/HXMwU4WRiuAB4BHjOz/wOsbOwfm1kJUAKQn5/PunXrUgpTUVGRch9BaY3s\na9asoUOHDnTt2pW9e/fSvXt3pkyZwnnnnZfSe+tzD4ayByOy2d09ow+gN7Ap5ng0sCDmeALwWEv6\nLiws9FStXbs25T6CkunsW7Zs8ZycHJ8+fXra+9bnHgxlD0bYsgOvegJ/Y4O4lLMLOD/muGddW8K0\nH0TruO+++8jNzdWMaZEsFUSBeAXoa2YXmFkuMA4oS6YDd1/p7iVdunTJSECBbdu28ctf/pKSkhLy\n8/ODjiMiAcj0ba5LgY3AxWa208wmu3s1cCvwe2Az8KS7v5XJHJK8+++/HzPjjjvuCDqKiAQko4PU\n7j6+kfZVwKqW9mtmI4ARffr0aWkX0oTdu3ezcOFCJk2aRM+ePZv/ByLSJkXydlJdYsqsBx54gJqa\nGmbMmBF0FBEJUCQLhGTOhx9+yLx587jhhhu48MILg44jIgGKZIHQXUyZ89BDD1FZWaltREUkmgVC\nl5gy46OPPuKxxx5jzJgxXHzxxUHHEZGARbJASGY8+uijHD58mLvvvjvoKCISApEsELrElH6HDh3i\npz/9KSNHjuTTn/500HFEJAQiWSB0iSn9fvazn/HRRx9pnwcR+VgkC4SkT2lpKb169eLOO+/kzDPP\n5J133gk6koiERBCruaZME+XSo3470fod4yorKykpKQGguFib+olku0ieQegSU3poO1ERaUokC4Sk\nh7YTFZGmqEBkMW0nKiJNUYHIYrNnzyYvL++UNm0nKiL1IlkgNA8iPYqLi5k/fz4FBQWYGQUFBcyf\nP18D1CICRPQuJndfCawcNGjQ1KCzRF1xcbEKgojEFckzCBERyTwVCBERiUsFQkRE4lKBEBGRuCJZ\nIHQXk4hI5pm7B52hxczsQ2A70AWIrRaxx/F+jm07B9jbgrdv+J6JPh+vvbm8jf0c5uyxx1HI3vA4\niOyJ5ozXluh3vrWzN5e1sbyxP2dz9qayNfd8U1kL3P3cZt/Z3SP/AOY3dhzv5wZtr6bjPRN9Pl57\nc3mb+D1Cm72JzzuU2RP8rmQ0e6I5E80e7/do7ezNZVX21vvOJ/I5NHxE8hJTHCubOI73c8PXp+M9\nE30+XntzeZv6uSVaI3vscRSyNzwOIntjr2lp9tifg8oer13Zk5Ou73zSWSJ9iSkdzOxVdx8UdI6W\nUPZgKHswlL31tZUziFTMDzpACpQ9GMoeDGVvZVl/BiEiIvHpDEJEROJSgRARkbhUIEREJC4ViCaY\nWUcze9XMvhJ0lmSY2SfNbJ6ZLTezW4LOkywzu9bMHjezX5vZPwSdJxlmdqGZLTSz5UFnSUTdd3xx\n3ecdqXXfo/ZZx4rKd7xNFggzW2Rme8xsU4P24Wa2xcy2mtmdCXQ1A3gyMynjS0d2d9/s7jcDY4DP\nZjJvQ2nK/7S7TwVuBsZmMm+sNGV/z90nZzZp05L8PYqA5XWf9zWtHraBZLKH4bOOlWT2QL7jSWvp\n7L4wP4CrgIHAppi2dsBfgAuBXOB1oD9wCfBMg0d34IvAOOAm4CtRyl73b64BngO+GrXPPubfzQEG\nRjT78oh8/+8C/q7uNUuCytyS7GH4rNOQvVW/48k+IrmjXHPc/Q9m1rtB8+XAVnd/D8DMlgEj3f3H\nwGmXkMxsMNCRk/8THTOzVe5em8nckJ7sdf2UAWVm9iywJHOJT3vfdHz2BtwHPOfuf8xs4r9K12cf\ntGR+D2Cwj/Q2AAACYklEQVQn0BP4MyG4opBk9v/Xuumalkx2M9tMAN/xZAX+hWhF5wHvxxzvrGuL\ny91nuvu3OPnH9fHWKA5NSCq7mQ02s0fM7OfAqkyHS0BS+YHbgGHAaDO7OZPBEpDsZ3+2mc0DLjOz\nuzIdLgmN/R4rgOvM7GekZ1mITIibPcSfdazGPvcwfccb1SbPINLJ3Z8IOkOy3H0dsC7gGC3m7o8A\njwSdoyXcfR8nrytHgrsfASYFnaMlovZZx4rKdzybziB2AefHHPesa4uCKGeHaOePcvZYUf49lD0g\n2VQgXgH6mtkFZpbLyQHosoAzJSrK2SHa+aOcPVaUfw9lD0rQo+QZuptgKfDfwAlOXvObXNf+ZeAd\nTt5VMDPonG0te9TzRzl7W/k9lD1cDy3WJyIicWXTJSYREUmCCoSIiMSlAiEiInGpQIiISFwqECIi\nEpcKhIiIxKUCISIicalAiIhIXCoQImlSt0uYm1m/mLaeZhbeDWFEmqACIZI+44FX6/5b72pObiIj\nEjlaakMkDcysEyfX2vkH4DfufpGZfQ74HXAAOAwUed3GMSJRoDMIkfQYCTzv7q8DFWZW6O7rObma\n50h3/zsVB4kaFQiR9BgPPFn385P89TLTxcDbgSQSSZEKhEiKzKwbcAWwuq7pSWCsmZ0LHHT36sDC\niaRABUIkdaOBVe5eBVB3Kem/gU8Cu4MMJpIKDVKLpMjM1gKXAodims8GlgMDgDygxN3LA4gn0mIq\nECIiEpcuMYmISFwqECIiEpcKhIiIxKUCISIicalAiIhIXCoQIiISlwqEiIjEpQIhIiJx/X/ijhRd\ne99axgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e1a2748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(step_sizes[:-1], errors[:-1], 'ko-')\n",
    "plt.axis('equal') # make a square a square (so a diagonal has slope)\n",
    "plt.grid(True) # Turn the grid on\n",
    "plt.xlabel(\"$\\Delta t$\")\n",
    "plt.ylabel(\"Error\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The straight line slope of 1 shows that Simple Euler is a first order algorithm. \n",
    "We only get below an error of $10^0$ (1.0) in this quantity of interest with the second smallest step size, $\\sim$0.002.\n",
    "The Euler method is not very good - it requires tiny steps and converges slowly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Homework\n",
    "\n",
    "## 1. Runke-Kutta (RK4)\n",
    "Implement the [Ruge-Kutta](https://en.wikipedia.org/wiki/Runge–Kutta_methods) RK4 method.\n",
    "First set up the right hand side of the ODE as a function, because you'll need to call it repeatedly within a single iteration.\n",
    "Determine the order of convergence (plot a `loglog` plot of *error* vs. $\\Delta t$ like above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Batteries included!\n",
    "Let's solve this with an even better ODE integrator. One that we don't have to write!\n",
    "Use `odeint` from the `scipy.integrate` package. Check its documentation (online) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:py3]",
   "language": "python",
   "name": "conda-env-py3-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}