New ODE solver for scipy

solve_ivp_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of usage of solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Restricted three body problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formulation is taken from an introductory chapter to Runge-Kutta methods in \"Solving Ordinary Differential Equations I\". The system of differential equations describes the \n",
    "motion of a body with a neglible mass in a gravitational field of two bodies with masses $\\mu$ and $\\nu = 1 - \\mu$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "y_1'' = y_1 + 2 y_2' - \\nu \\frac{y_1 + \\mu}{D_1} - \\mu \\frac{y_1 - \\mu'}{D_2} \\\\\n",
    "y_2'' = y_2 - 2 y_1' - \\nu \\frac{y_2}{D_1} - \\mu \\frac{y_2}{D_2} \\\\\n",
    "D_1 = ((y_1 + \\mu)^2 + y_2^2)^{3/2}, \\quad D_2 = ((y_1 - \\nu)^2 + y_2^2)^{3/2} \\\\\n",
    "\\mu = 0.012277471\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial values are picked such that the solution is periodic with period $x_{end}$:\n",
    "$$\n",
    "y_1(0) = 0.994, y_1′(0) = 0, y_2(0) = 0 \\\\\n",
    "y_2′(0) = −2.00158510637908252240537862224 \\\\\n",
    "x_{end} = 17.0652165601579625588917206249\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is straightforward to rewrite this systems as another first order system. We define the constants and evaluation of the system right-hand side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 0.0122771\n",
    "nu = 1 - mu\n",
    "\n",
    "a = 0\n",
    "b = 17.0652165601579625588917206249\n",
    "ya = [0.994, 0, 0, -2.00158510637908252240537862224]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(x, y):\n",
    "    D1 = ((y[0] + mu)**2 + y[1]**2) ** 1.5\n",
    "    D2 = ((y[0] - nu)**2 + y[1]**2) ** 1.5\n",
    "\n",
    "    return [\n",
    "        y[2],\n",
    "        y[3],\n",
    "        y[0] + 2 * y[3] - nu * (y[0] + mu) / D1 - mu * (y[0] - nu) / D2,\n",
    "        y[1] - 2 * y[2] - nu * y[1] / D1 - mu * y[1] / D2\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ODE is integrated by a single call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, [a, b], ya, rtol=1e-3, method='RK45')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `solve_ivp` provides a continuous solution with the *same accuracy* as underlying method. Here we use is to produce a smooth plot with the compued trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xp = np.linspace(a, b, 500)\n",
    "yp = res.sol(xp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x107b6b0f0>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAHfCAYAAAD3H2TtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VOXWxuHfDoQAGpASUKQoSMsEwYZBDQSxFxQVFazY\nUMSuKOpR+PRYsHeOigUVFVEUC1IUTMAIiihmQu+d0AcIIWV/fywxgJQQJrOnPPd17SuTMMwsIEye\nect6Hdd1EREREZHQifO6ABEREZFYowAmIiIiEmIKYCIiIiIhpgAmIiIiEmIKYCIiIiIhpgAmIiIi\nEmJBCWCO4wxyHGel4zjT9vDrHRzHWe84zu9/Xw8H43lFREREIlHFID3Ou8ArwOC93CfDdd3OQXo+\nERERkYgVlBEw13UnAOv2cTcnGM8lIiIiEulCuQasneM4fziO863jOMkhfF4RERGRsBKsKch9mQI0\ndF13i+M4ZwNfAs12d0fHcXQ2koiIiEQM13X3e5YvJCNgrutucl13y9+3RwLxjuPU3Mv9dUXg9eij\nj3pegy79+8XipX+7yL707xfZV1kFM4A57GGdl+M4dXe43RZwXNddG8TnFhEREYkYQZmCdBxnCJAO\n1HIcZxHwKFAJcF3XfRO4xHGcW4ACIA+4LBjPKyIiIhKJghLAXNftvo9ffw14LRjPJeErPT3d6xLk\nAOjfL3Lp3y6y6d8vNjkHMn9ZHhzHccOtJhEREZHdcRwHN1wX4YuIiIhICQUwERERkRBTABMREREJ\nMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMR\nEREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBT\nABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERER\nkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUw\nERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJ\nMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJMQUwERERkRBTABMREREJsYpeFyAiUlaFhZCf\nb9e2bXZtv11QAPHxUKkSJCTYx+23ExKg4j5e/QL5AbJXZZNSJ4XEhMTQ/IFEJGYogIlIWNm8GRYu\nLLmWLoU1a0qutWtLbufllQSqHYPW9oC1PaDtGMy2B7aDD4ZataB27Z0/1q0LdRsGeHJlGgs2+0lO\n8jGhR6ZCmIgEleO6rtc17MRxHDfcahKR4Nq6FWbMgOxsu2bPhgULLHBt3gwNG0KjRnbVr2/BaMer\nZk37mJgIjrP/z19cDBs3wurVFuR2/LhiBUxZmcWPjdpDhUIojKdJZga+6qm0bAlt2sAxx8BRR0GF\nCkH/qxGRCOM4Dq7r7vcrkQKYiJSr5cvhl19g6lTw+y1wLVoETZpASgr4fNC8uYWtI46AOnXKFqqC\nKZAfIO3dNHJyc2heM5mBbTNZsSiRnBz44w/7s+TmQqtWFsiOPx46dIDGjb2vXURCSwFMRDyXn28B\nJSvLQldWFmzaBKmpcOyxFrhSUqBpU5sqDGeB/AD+XD++JN9upx/Xr4c//7Q/76RJMH68rTnr2LHk\natgw9HWLSGgpgIlIyBUVwZQpMHq0XVOmQLNm0K6dha527WyqLhZGhVwXZs2CceNKrkMOgS5d4JJL\nbJQsFv4eRGKNApiIhMTChSWB68cfoV49OOMMu04+2Ra3iwWyqVPh889h2DBb93bxxRbGUlMhTk2A\nRKKCApiIlAvXhZwcCxHDhsGqVXD66Ra4TjvNApjsneva+rdhw+Czz2yq9uaboUcP20wgIpFLAUxE\ngsZ1Ydq0ktC1aZON3HTtqtGbA+W6tmbs9ddhxAi44ALo1QvattUUpUgkUgATkQO2dCm8/z689541\nMr3kErtOOEGhqzysXm1/12+8YX3IHnoIzj9fQUwkkiiAiUiZFBTAt9/CoEEwcaKNcl13nUZkQqm4\nGIYPh8ces7/zJ5+EM8/U379IJFAAE5H9Mn8+DBwIgwfbTsXrr7fwddBBXlcWu1zXgljfvtCgAbzw\ngvUaE5HwVdYApkkFkRjz669w6aU2rVhUZO0SMjPh2msVvrzmOHDRRdas9qKLoFMnuPdeW4MnItFF\nAUwkBhQXw9dfW7f2Sy6Bk06yEbBnn4UWLbyuTnYVH28L87OzYeVKaN0aJkzwuioRCSZNQYpEsaIi\nGDIEnngCqla10ZSuXe2gaokcX31lbStuugkeeURnUIqEE60BE5F/uK6NeD34IFSvDv/3f3DqqVrU\nHclWrIDLL4fKlS1U16zpdUUiAgpgIjEvkB8ge1U262el8Nh/EgkEbOTrvPMUvKJFYSHcdx+MGgUj\nR9oB5iLiLQUwkRgWyA9wwhtpzFrnp+J6H68dl8l1VyRqqipKvfQSPPMM/PADNG/udTUisc3TXZCO\n4wxyHGel4zjT9nKflx3Hme04zh+O47QJxvOKCOTlwe2PZzNzrR83rhBq59DqVL/CVxS74w54/HE7\nCmr2bK+rEZGyCNYuyHeBM/f0i47jnA00cV23KdATGBik5xWJaaNGQUoKrJ+dQnKSj/i4eJKTkvEl\n+bwuTcrZtdfCww/D2WdbR30RiSxB2Qvluu4Ex3H2thrhAmDw3/ed5DhOdcdx6rquuzIYzy8Sa1as\nsFGQ336DV1+Fs89OJJCfiT/Xjy/JR2JCotclSgj07Anz5sHFF9t0pHa3ikSOUPUBOxxYvMPnS//+\nmojsp+++g2OOgcaNrU/U2Wfb1xMTEkmtn6rwFWOefBISEmxKUkQiR1i+X+rXr98/t9PT00lPT/es\nFpFwsXUr9OljPaGGDoW0NK8rknAQF2cHqB9zDHTuDMce63VFItFt/PjxjB8//oAfJ2i7IP+egvza\ndd2jd/NrA4Fxrut++vfnM4AOu5uC1C5IkX/Lzobu3aFlSzu/sUYNryuScDNoELz9th2oHqczTkRC\nJhzOgnT+vnZnBHA1gOM4qcB6rf8SKZ3Bg6FjR7jrLvjkE4Uv2b0ePSCvKMCAIVkE8gNelyMi+xCU\nETDHcYYA6UAtYCXwKFAJcF3XffPv+7wKnAVsBnq4rvv7Hh5LI2AilDTd/OYb+PJL8Gljo+xFID9A\nqxfTWLjZz9H1fEzokan1gCIhUNYRsGDtguxeivv0DsZzicSCTZvg0ksthE2apGNnZN+yV2WzdJsf\nKhSSsyoHf66f1PqpXpclInuglQIiYWbFCujQAerVg2+/VfiS0kmpk4IvyUecG0+NQvWCEwl3CmAi\nYWTxYmjfHs4/H956C+Ljva5IIkViQiKZPTJ566QMao3Q9KNIuNNZkCJhYt486NQJbr/dFtyLlEVh\nIdSuDXPm2EcRKV/hsAtSRMpo8WILX/fdp/AlB6ZiRWtXMn2615WIyN4ogIl4bNUqO1S5d2/o1cvr\naiQaNGkCc+d6XYWI7I0CmIiH8vLgggvgkkvgnnu8rkaixSGHQECtwETCmgKYiEdcF665xs501Dl+\nEkwJCZCf73UVIrI3YXkWpEgseOYZWLgQMjLA2e/lmyJ7FghAs2ZeVyEie6MAJuKBjAx44QWYPNlG\nK0SCKTcXatXyugoR2RtNQYqE2IYNcPXVdnhygwZeVyPRaMYMaN7c6ypEZG/UB0wkxHr0gMqV4Y03\nvK5EotHmzZCUBOvWaXRVJBQ8PQtSRErnp59g7FjIyfG6EolWGRlwwgkKXyLhTlOQIiFSWGh9vl58\nERJ1SoyUk++/h9NP97oKEdkXjYCJhMjgwXY0zEUXeV2JRKuCAvj0UxsFE5HwpgAmEgLbtkH//jBk\niFpOSPkZOdL6yqkFhUj40xSkSAh89pkdD3PyyV5XItHsuefgttu8rkJESkMBTCQEXnxRh2xL+Zo4\n0Q5179rV60pEpDQUwETKWXY2LF8O55zjdSUSrVwX+vSB//wHKmphiUhEUAATKWdDhkD37lChgteV\nSLQaNgy2bLEGvyISGRTARMrZN99Aly5eVyHRav16m95+6SWFfJFIok74IuVoxQpIToZVqzQ1JOXj\n+uvtZIXXXvO6EpHYpE74ImFo8mQ48USFLykfX34JP/4I06Z5XYmI7C/9WBApR3/8AW3aeF2FRKMF\nC6BnTxgxQicriEQirQETKUfz5sFRR3ldhUSbvDxrN9Gnj42wikjkUQATKUfLlsHhh3tdhUQT14Ue\nPaBpU7j7bq+rEZGy0hSkSDnavFnTQxJc/fvb9OO4cTrWSiSSKYCJlKPCQrUGkOAZOBA+/BAmTIAq\nVbyuRkQOhKYgRcpRpcQAvy7PIpAf8LoUiXBDh8Ljj8Po0XDooV5XIyIHSgFMpJwE8gP8dXwad01r\nT9q7aQphUmbDh9sh2999B40be12NiASDAphIOclelU2gsp8iCsnJzcGf6/e6JIlAX3wBt9wCI0fC\n0Ud7XY2IBIsCmEg5SamTQv0EH05xPMlJyfiSfF6XJBEkkB/g8feyuPmOACNHwrHHel2RiASTjiIS\nKUd/zQrQ4RI/8yf5qF5F2yGldAL5AVo+m8bSfD/Navj4rVcmiQn6/hEJR2U9ikgjYCLlqFWzRGrl\npTJ3un54Sum4Ltz7TDZL8/1QoZD5mzR9LRKNFMBEytlll1nrAJF9KSyEm2+GicNTSE7yER+n6WuR\naKUpSJFyNns2nHKKNc9U7ybZk40b4dJLrbnqp5+CkxDAn+vHl+TT9KNIGNMUpEiYatoU2rWDt9/2\nuhIJV/PnQ1oaHHEEfP01VKsGiQmJpNZPVfgSiVIaARMJgd9+gwsugBkzdDSR7GzMGLjqKujbF26/\nXccLiUSaso6AKYCJhEiPHlCzJjz3nNeVSDhwXXj2WXj+efj4Y0hP97oiESkLBTCRMLdqlTXS/Pxz\nOPlkr6sRL23eDNdfD3PmWKPVhg29rkhEykprwETCXJ06tg7siitgzRqvqxGvzJsHJ50ECQmQmanw\nJRKrFMBEQui886wtRZcukJ/vdTUSaiNHWvi68UZ47z3tihWJZZqCFAmx4mJrN1ChAnz0EVSs6HVF\nUt4KCuChh2yt15AhtuNRRKKDpiBFIkRcHHzwAWzYAFdeac03JXotWGCBKycHpk5V+BIRowAm4oEq\nVeDLLyEQsPYUgYDXFUl5+PxzaNvWRjxHjIDatb2uSETChaYgRTxUUAC9e8Mvv8BXX1kjTol8W7fC\n3XfD99/DJ59YCBOR6KQpSJEIFB8PAwfCddfBiSdaSwKJbDNnQmoqrF5tU44KXyKyOwpgIh5zHLjj\nDjuC5t574YYbYP16r6uS/eW68Oab1uPtllvsPMfq1b2uSkTClQKYSJho2xb++AMqVQKfz9YPaTY+\nMqxaZWv53ngDMjKgZ08dKSQie6cAJhJGqlWD11+3dgX9+sGpp1ook/D17bfQpg0kJ8OkSfZRRGRf\ntAhfJEwVFlrn/H79oFMn+M9/oEULr6uS7bZssSnj776DwYOhfXuvKxIRL2gRvkiUqVgRbr4ZZs+G\nlBT7Ad+tG/z+OwTyA2QtziKQr/4VXpgyBY491tqH/PmnwpeI7D+NgIlEiEDAdky+NDDA+i5pbE30\n46vjY8J1mSQmJHpdXkwoKoKnn4YXX4SXX4bLL/e6IhHxWllHwBTARCJM5oIsOr7fniIKoSieK7Zl\n8PC1qZqeLGfz58NVV9kmifffhwYNvK5IRMKBpiBFYkSbw1JIqesjPi6e5rWSqe366NgRjj8ennkG\nFi70usLo4rq2xqttWztEfexYhS8ROXAaAROJQIH8AP5cP74kH4kJiRQWwvjxMHSoNXNt3BjOP9+u\n1q3VEqGs1q61dXg5OXZweuvWXlckIuFGU5AiAtjxRhkZ1h7h668hLw/OPNN2Up56Khx6qNcVRoax\nY6FHD7jkEnjySahc2euKRCQcKYCJyG7NnAljxsAPP8BPP8Fhh0GHDnDSSXYdeaRGyHaUlwd9+1oj\n3HffhdNO87oiEQlnCmAisk9FRXY+4YQJ8PPPMHGifa1tW2sm2rq1fTzySIiLwRWif/4JV1xhzVQH\nDoSaNb2uSETCnQKYiOw317VF+1OmWMf9P/+0j+vXw9FHWxjbHsx8Pqha1euKy0dRETz/PAwYYB+v\nvFKjgiJSOgpgIhI0a9bAtGkWxrYHsxkzoHZtaNoUjjpq549NmkRuOFu0CK6+GoqLbbfjEUd4XZGI\nRBIFMBEpV0VFsGSJdeafPRvmzCn5OH8+1KpVEsYaNNj5OqRugPmbskmpkxJWTWOHDIE774S774b7\n7oMKFbyuSEQijQKYiHimqAgWL7YwNneu3V682ALbwuUB5qWn4db2Uzngo/2cTI48PJEGDaB+fdsU\nULeu7c6sXTs0IWjdOujVy0b2PvzQjhUSESkLBTARCUtZi7No/157CosLqejE80TTDKquTWXJEgtp\nK1bAypX2cf16W/h+6KEloaxu3Z1vb/9Yq1bZwtqPP8K118KFF9qxQlWqBP2PLCIxRAFMRMJSID9A\n2rtp5OTmkJyUTGaPPZ9dWVgIubkloWx7MNvx4/bbGzbYiNmuQW13V+3a9tgPPQQffwzvvGO90URE\nDpQCmIiErV079wdDQcHOYW3HgLbrtXp1ye875hhrM7GnsJaUBPHxQSlRRGKApwHMcZyzgBexsyUH\nua779C6/3gH4Cpj395e+cF338T08lgKYiASF68Krr8Ijj0CfPnD22XsOaduvNWugevWdQ1m9erZe\nbcfrsMOgYsV/P2cgP0D2qvDbcCAi5aOsAWw3Lx/7/cRxwKtAJ2AZ8KvjOF+5rjtjl7tmuK7b+UCf\nT0SkNNasgeuug6VL4ddfrWVGaRQV2e/dMZQtWwYLFlgD2yVL7MrNtdGyHUNZ7XoB3i5OY1mBn6Oq\n+8jskUntagphIvJvBxzAgLbAbNd1FwI4jvMJcAGwawBTW0MRCYmMDGum2rUrfPYZVKpU+t9boQLU\nqWNXq1Z7vl9BgU17Ll1aEsomL8tmcRU/blwhM9bkcFhrPzU2p9KwoZ0u0Lixfdx+NWoECQkH/ucV\nkcgTjAB2OLB4h8+XYKFsV+0cx/kDWArc57puThCeW0QEsKm/P1dk8917KbzzRiLvvAPnnFN+zxcf\nX9LnrKSGFHLe9dmGg7rJ/JTjI2+DnTYwf75dU6bAsGF2e8kSC3qNG1sPtWbNSq4mTRTORKJZMAJY\naUwBGrquu8VxnLOBL4Fme7pzv379/rmdnp5Oenp6edcnIhEskB/gxDfTmJ7r56A8H7//kkmzI0I/\n9ZeYkEhmj8ydNhxUr2K7NE888d/3Lyy0EDZ3rjW1nTXLRu9mzbLQdvjhJYGsRQtISbEjoWrVCvkf\nTUT+Nn78eMaPH3/Aj3PAi/Adx0kF+rmue9bfnz8AuLsuxN/l98wHjnNdd+1ufk2L8EVkvwz4OIv7\np7eHCoXEx8WT0SOD1PqpXpd1QAoKbN3ZrFl25eSA329X1aolYSwlxa7kZKhWzeuqRWKPZ7sgHcep\nAMzEFuEvByYD3VzXnb7Dfeq6rrvy79ttgaGu6x6xh8dTABORUtm61XY3fvldgIRb0li4Zd+9xiKd\n69qoWXa2hbHsbLumT7fpzGOP3fmqW9frikWiWzi0oXiJkjYUTzmO0xMbCXvTcZxbgVuAAiAPuMt1\n3Ul7eCwFMBHZp9mz4dJLba3UW29BxarB7zUWSYqKbCrz9993vipX3jmQHX+87dgUkeBQI1YRiRnD\nhsEtt0D//vbR0R7r3XJdWLRo50A2ebKFsnbtIDXVPh57rBb8i5SVApiIRL2CArj/fhg+HIYOhRNO\n8LqiyOO6MG8eZGXZ9csvMGOGtdzYHspOOcU2AIjIvimAScRzXdsVVlho0ynx8da/SaMbAtZv69JL\noUYNGDzYDu2W4Ni8GX77zcJYVpY1nK1RA9LTSy4FMpHdUwCTsLZhA8ycabu5Zs+2RcTLlsHy5bBu\nHaxfD4GANcGsWBHi4my0o6DApktq1rSt93XqQMOG1sCycWPb/dWiBVSp4vWfUMrTjz/CFVfAbbfB\nAw/Y94eUn+JiW+A/frxdP/2kQCayJwpgEjY2bYJJk+Dnn23NyR9/2GHI2/sZNW1qzSvr1bOrZk07\ney8x0QLYjoqLbafb2rUlx8MsWmQ9kmbPth8Sc+ZYGGvXzq4OHUp/7IyEN9eFF16AAQPgo4+gUyev\nK4pNuwtktWrBWWfZ1aGDtcYQiUUKYOKZbdts2mL0aBgzxvoVtWkDJ51kO66OOcZ2qpXXqEVBgW3D\nz8qy0PfjjzYidsYZcOGFcOqpNp0pkWXLFrjpJvt+Gj7cRj0lPBQXw59/wvff2/X773DyySWBrHlz\nLR2Q2KEAJiETyA8waUE2S39PYdTXiYwcaaNaZ5wBp59uHb8rV/auPte1d+sjR8IXX9hI2YUXwrXX\n2g8J/WAIfwsWQJcu1mj0zTc1uhLuNmywNz7ff2//7+LiLIidc469JmiJgEQzBTApd4WF8OXIADf+\nnMb6eD8Hb/Xx38aZXHphIoce6nV1e7ZoEXz6KQwaZJ/fdBPccIO6hoerceOgWzdrsHrXXQrMkcZ1\nbVflyJHw9dc2OnbaaRaozz3X1pKJRBMFMCk3Cxdao8t33oGarbOYntqeYiLvyBfXhYkT4dVXbar0\nuuvsB3y9el5XJtv973/wyCMwZIjWe0WL1avhm29sGnncOBsh79LFRqXr1bMR9exV2aTUSYnJBroS\n+coawLSXSHbLdW0repcu1qQxELDQkvVlCq3q+oiPiyc5KRlfks/rUkvNcay/0SefwJQptnYsJQVu\nvzfAyOwsAvkBr0uMWYWFcOed8Pzz9n2n8BU9ate26f+vvrJdzzffbOs1U1LghJMDNH86jfbvtSft\n3TT9H5SYohEw2Ynr2rvV//7X3rnedZe9eB50UMl9AvnRc+TLzPkBTh6Uxpo4P4dV9PHXXZnUSozs\nP1Ok2bDBphwLCqy5qqaoYkNBAbz+dRZ3T2tPsVOIUxRPnzoZ9Omeqh5vElE0AiYHxHVhxAg47jh4\n+GG4917r23XrrTuHL4DEhERS66dGfPgCWFsxmw0JfqhQyIqiHI4508/o0V5XFTvmzbPdskceCd99\np/AVS+Lj4bpzU2h1qI2oNzoomZkTfBx5JHTubCPVmzd7XaVI+dEImJCVZQue16+Hxx6zF79YaXQZ\nyA+Q9m4aObk5JCcl81C9TO67PZFTT7XpsEMO8brC6DV5MlxwATz0EPTu7XU14pVdR9Q3brTpyiFD\n7LXpvPNsFP7UU2PndUkiixbhy35bssRGuiZOtOB11VX/boQaC3b9ARAIWCD95ht4/3174ZfgGjEC\nrr8e3n3XfsCK7M6qVSU7mDdsgB497GrQwOvKREoogEmpFRTAc8/Bs89Cr152tIv6LP3b6NH2zvuG\nG2xnXsWKXlcUHQYOhP79bZSjbVuvq5FI4LrWzmLQIAtkJ5xgAb5zZ0hI8Lo6iXUKYFIqU6ZY+4V6\n9awdQ5MmXlcU3lasgCuvtB8AQ4fa8StSNsXFNt34+efWI0rfe1IWeXn2PTRokDVcvuIK6+3XsqXX\nlUms0iJ82att2+yH3znnwH332YJn/QDct0MPhVGjbHNC27b2gi/7r6AArrnGzhGcOFHfe1J2VarY\nm6Jx42yNWNWq0LGjtS4ZPtxamohEAo2AxYAZM+wF67DD4O23oW5dryuKTIMH25q5Tz7RurD9kZcH\nXbva7aFDNd0twZefb6Nir75qa1tvvtmWDtSp43VlEgs0Aia79cEHkJYGN95oC58Vvsru6qstQFx+\nOQwb5nU1kWHDBjsTsHp1G51Q+JLykJAA3bvDzz/Dl1/C3LnQrJn9n5082evqRHZPI2BRautWuOMO\nm/L5/HPrOi3B8ccfdqbd//2fLQSW3cvNtfCVmgqvvKIWAhJaa9bYLtvXX7eRsHvvtZM9YnGnt5Qv\nLcKXfyxfbuesNWxo5zeqsXvwzZ5t05D9+imE7c7ixXD66XDxxfD44zpQW7xTVGQ7bp95xtpa3H33\nv0/3EDkQmoIUwLZqn3ginH++TZcpfJWPpk3hxx8tgL37rtfVhJf586F9e1uD89//KnyJtypUgIsu\nsgX7gwfD2LF28sIjj8DKlV5XJ7FMASyKjBplUz4vvGDHCekHX/lq2hR++MF2l375pdfVhIe5cyE9\n3aZ77r3X62pEdnbyybYWccIEGw1r0cJaWMyZ43VlEosUwKLEBx/YgtMvv7RpHwmNZs3g66/tRfyn\nn7yuxltz5lg7gL597QxRkXDVrJk1BJ4501rNpKbaTvGcHK8rk1iiABYFBg60UZhx4+xgYwmt446D\njz+2VgvTp3tdjTdmzbLw9Z//WAsAkUhQp45tppk7F3w++x7u2tU22oiUNwWwCPfSS/D007bbMTnZ\n62piV6dOMGCArb1bs8brakJr5kzbkNC/v7U7EYk01avbyO28eTYads45dsyRWlhIedIuyAj2+ut2\nnuO4cdCokdfVCNgh3pMmwZgxUKmS19WUv3nzoEMHO8z92mu9rkYkOPLy7Kijp5+2Fj6PPQbHH+91\nVRKutAsyxgweDE8+aYvAFb7Cx5NP2s7TPn28rqT8LVkCp51mGz4UviSaVKkCvXvbusbzz4cLLrCd\nlNnZXlcm0UQBLAJ99539gB892rZTS/ioUME2RHz9NXz6qdfVlJ+VKy183Xor9OzpdTUi5SMhAXr1\nsiB28sm21KB7d+sDKHKgFMAizJQpNtrw5ZfQsqXX1cju1Khhpw/07h2di/LXroUzzoBu3eCee7yu\nRqT8Vali3+tz5thi/ZNOsgbMCxd6XZlEMgWwCLJ0qS0M/d//bKGohK82bawJabdudlBwtAgE4Oyz\nLYA98ojX1YiEVmKi7TifNQsOOwyOPdb63a1d63VlEokUwCLE1q22BqF3bzvPTMLfjTdCkya2uyoa\nbNtm34MvWHB8AAAgAElEQVStW9uOTzX6lVhVo4YdsZWdDZs2QfPmdtTR1q1eVyaRRLsgI8T118PG\njXa8kH7wRY61a2007K234Mwzva6m7IqLrdFvIGDTqxUrel2RSPiYMcPeaP3+u+2YvOIKHfodS3QY\ndxT74AObzvrtNzj4YK+rkf01dqwF6L/+gmrVvK6mbO6/HzIz7c9StarX1YiEp4kT4b77YPNmaxF0\n+uleVyShoAAWpWbNst03P/wARx/tdTVSVjfcAPHx8MYbXley/15+2eqeMAFq1fK6GpHw5rp23mSf\nPrZR6tlnbYpSopf6gEWhwkI7n6xfP4WvSPfss9aaYvx4ryvZP599Zuu9vv9e4UukNBzH1kr6/dak\n+JRT4I47tFBf/k0BLIw9/bQdkdGrl9eVyIE65BAbRbrhBpueiASTJ9v33jffqNmvyP5KSLAdkjk5\nUFAALVrYaPLaTQGyFmcRyA94XaJ4TFOQYcrvh/R06/vVsKHX1UiwdO8Ohx4Kzz/vdSV7t3gxtGtn\nx1117ux1NSKRLzsb7rgvwMTmaRTW8JNS10dmj0wSExK9Lk0OkNaARRHXtaHryy6zTuMSPVavhlat\n4MPPAlRtlE1KnZSwewHevNmmTbp3twXFIhIcPy/Kov277SmiEKc4nuGdM7jgODV1jHRaAxZFBg+2\nw2BvvtnrSiTYateGJ54LcO7nabR/rz1p76aF1VTEhrwA59yURcqxAe691+tqRKJLq7oppNT1ER8X\nTx0nmevO8/HUU9ZjT2KPAliY2bTJ+sm8/rr6yESr5qdks62an8LiQnJyc/Dn+r0uCYBAfoBmT6WR\n0aQ9fx6XxqZt4RMMRaJBYkIimT0yyeiRwey+mUyekMjEiTYqPnq019VJqCmAhZkBA+DUU+GEE7yu\nRMpLq7opNKvhg6J4mtZIxpfk87okAF4Zms0q1w8VCpmxJnyCoUg0SUxIJLV+KokJiTRpYrujn3sO\nbrkFLr5Y50vGEgWwMLJiBbz2GjzxhNeVSHlKTEjk11syuWxLBsdPC49FuDNmwHP3p9D0EJseSU4K\nn2AoEu3OO882XrVubedLPv64jjWKBVqEH0buvtuOfHnxRa8rkVDYsMG2po8Y4e2IZyAAJ55o33+X\nXRXAn+vHl+QLi2AoEmvmz7f/i9nZ8NJLcM45Xlck+6JdkBFuxQpITrZ3QYcd5nU1EirvvWf9wbKy\nIM6D8WjXtd221arB22+H/vlFZPe+/x5uv9266b/4Ihx5pNcVyZ5oF2SEe+klO8BV4Su2XH21Ba/3\n3/fm+V94AebNg1df9eb5RWT3zjrLzo898UQbIe/f33bHS/TQCFgY2LQJjjjCOo83bux1NRJqU6bY\nGpCZM0N7WPfkyfa8kyfb95+IhKdFi2xacupUGw07/3yvK5IdaQQsgr3/PrRvr/AVq447zt7tPv10\n6J5zwwbo1s2mPxW+RMJbw4YwbBgMHGjNkc87D+bM8boqOVAaAfOY69rOlxdegE6dvK5GvLJ0qX0f\nTJlS/ucuuq6Frxo1LICJSOTYts1GwQYMsNYVfftC1apeVxXbNAIWoX791Y5+6djR60rES4cfDr17\nw4MPlv9zvfOObfYI9/MoReTfKlWCPn3gjz9g9mzbvDV8uL2xksiiETCP9eplP3wfesjrSsRrmzdD\ns2b2Ytq2bfk8x/TpNt3900/2wi0ike3HH+G226BBA3j5ZXsNkdDSCFgEKiy0ef1u3byuRMLBQQdZ\nA8a77y6fd7MFBXDVVfDYYwpfItHi1FNtNOz00+Gkk2wUffNmr6uS0lAA89C4cbYAWovvZburr7Zd\nsV98EfzHfuwxqFsXevYM/mOLiHfi4+Gee2DaNDvKqGVLe3MfQ5NJEUlTkB669Vbb3XL//V5XIuHk\nhx/gppsgJwcSEoLzmJMmQefO9k5ZveZEoltGhq0prVvXFuz7dKpYudIUZIRxXfj2Wzj3XK8rkXDT\nqZO9YL72WnAeb8sWG1l79VWFL5FY0L49/P679Qvr2NE66q9d63VVsisFMI9M+StAXu0sGh4V8LoU\nCUMDBsCTT8Lq1Qf+WPffb520u3Y98McSkchQsaIFr5wcKCqyacnXX7e1xxIeNAXpgUB+gJbPpLGs\n0M/Rh/rI7JGpg4/lX267DRzHdjaVVUYGdO9uR5rUqBG82kQksvz1F9xxB+Tm2rSk+k4Gj6YgI0j2\nqmyWFfhxnUJycnPw5/q9LknC0KOPwscf2xFFZZGXBzfcYFOZCl8isa1VK1tf2r8/3HgjXHSRnQMr\n3lEA80BKnRQqrvNRMS6e5KRkfElaISn/Vru2TR/26VO23//oo3DssXDBBcGtS0Qik+NY8MrJsWUJ\nbdta24qAVsJ4QlOQHlizBo5sHuD73/20quvT9KPsUX6+rd0YNGj/Tkv49VdbgPvXX5CUVH71iUjk\nWrYMHngAxoyxZuA33WSd9mX/aAoygkydCm2SEzmpYarCl+xVQoId0n3PPVBcXLrfU1AA119vRw0p\nfInIntSrB4MHw8iR8N130KIFfPRR6V9r5MAogHkgO9vm40VK45JLoEoV+OCD0t3/xRfteCudsCAi\npdGmjQWwd9+1djXHHGOfR/lklOc0BemB3r2haVPbkSJSGr/8YkFs5kw7smhPFi+2F89ffoGjjgpd\nfSISHVwXvvrK1obVrg1PPWVHHMmeaQoygsyZYwFMpLRSUyEtDZ59du/3u+suC/gKXyJSFo4DF15o\nxxr16AGXX24bef76y+vKoo8CmAcWL7YjiET2x5NPwiuv2MLZ3fn+e1tfqKOtRORAVaxoAWzWLOus\nf/rptoNy6lTrZZm1OItAvrZPHghNQXqgVi37pq5Vy+tKJNL07QsrV8I77+z89a1bISXFmraec443\ntYlI9NqyBd58E55+McDmy9LIO8iPr44aiYPHU5CO45zlOM4Mx3FmOY6z2/ffjuO87DjObMdx/nAc\np00wnjcSbdtmPVfUGFPKom9f27H0xx87f33AADj6aIUvESkfVavCnXfCxz9ks7mqn0K3kGnLchg4\n3K9dk2V0wAHMcZw44FXgTMAHdHMcp8Uu9zkbaOK6blOgJzDwQJ83Um3cCNWqQZwmf6UMqlWzBqv3\n3FOyQ2nuXHjpJdv9KCJSno6rn0Kruj7i4+JpUDWZj573kZwMb71lp29I6QUjBrQFZruuu9B13QLg\nE2DX3tsXAIMBXNedBFR3HKduEJ474mzatPddbCL7csMNsGIFfPuthbDbb7du+VpXKCLlLTEhkcwe\nmWT0yCD7rkymTkrkjTds52SDBvbmcM6c8q1hd2vQdv3ajp/v6/ayjcsYO28sY+eN/ef2iJkjGDtv\nbLmuc6sYhMc4HFi8w+dLsFC2t/ss/ftrK4Pw/BFlyxYbyhUpq4oVbTfkXXfZO87582H4cK+rEpFY\nkZiQSGr91H8+79jRrnnz4H//g3bt7Bi0m26C886zhtL7K5AfIHtVNil1UnZaY7Zxa4B2b6cxa62f\nwyv5uKFCJqtWwoeV0lgf76fKZh+Hj/+OhSefw7bqfiqsb4HrQnGNGVRY3wLHgcLqM6gUaEGFCpB3\n0HQcKuI6+eBABSpRxLZ/nq9pjZb07/AIizctIqlKEksCS8CF+ArxjJ43hv4d+pX57zEYAUz2g+tq\n+lEO3Fln2a7ISy+1A3Z1fIiIeK1xYzu5o39/GDrUdm3ffLO1srjmGjjuOGtzsS+B/ABp76bhz/XT\nsLKP6+MymZuTSHY2/LUum7zL/VChkEV5OWRv9lOrlsuGbX5cCtlWLYduj3zLE3/4wS2E2tNxcP65\n7f59u+iQ6RTiAEW4FP3z3EXuNtihxtlrp9N9+J67Wo9b8GOZ/76CEcCWAjtOftT/+2u73qfBPu7z\nj379+v1zOz09nfT09AOtUSSqOA7Ur2+3jznG21pERHZUuTJcfbVd8+fbcUeXXmojYV27WlPpVq12\nDmOuCzNmwE8/wYip2fxZ10LW/EAOUwN+TmubSo8e0LBpCp2H+8jJzSG5bjJv9fABMPHdv7+WlEzP\n9HP5apl93rxWcwBmrpm529szVs+gQlwF8ovyAahUoRLbiktGwNhdYJwPLDjwv6cDbkPhOE4FYCbQ\nCVgOTAa6ua47fYf7nAPc6rruuY7jpAIvuq6buofHi+o2FNOnQ5cu9o0mUlYrV4LPZ+8oGzeGN97w\nuiIRkT0rLobJk2HYMLsSEqzDflwcbN4MGRk2kt+xIxxzYoDX8tKYv8kC1a6tLgL5Afy5fnxJvn++\nvuvXdvwc2OvthtUaMn21RZaWtVsyffV0thRsAeDe0fcya+0se+Lt0WTXUNaPMrWhCEofMMdxzgJe\nwhb1D3Jd9ynHcXoCruu6b/59n1eBs4DNQA/XdX/fw2NFdQBbvNjmx5cs8boSiWQ9e8LBB8N//gPJ\nyfDll9B215WXIiJhYPt6Ll9SCovmJPLVV/DIIzsf+n300dCrF3ToAM2bw6ZtgX/C0cINC/+1FixU\nNVeLr8Up/zuT9UULqFqhNltYvdP9jjzkSObfOd+7ABZM0R7A1q2DI46ADRu8rkQi1V9/QadOdi5k\njRrw4Yfw/PP27rKiVnWKSBgJ5Ac44Y00Zq/3U3G9jzpfZ9Ll3EQuuMCOV6tUyUb0x4yBUaNgwgTr\nlXnyyXD8SQHei0tj0VYbsSqPpq+7LvZ3XfjdH+DCEWks2+YnbuMRFCbOh7gi4uPiaVS9EQs3LKRB\ntQa8eNaLpB+RTrXK1RTAIkFRkc2Pb9kC8fFeVyORxnXhjDOgc2e47baSr3XqZOe33X67t/WJiIC1\nyvngA3hzZBZzTmkPFQqp6MST0SODdg12uwLpH0uXwsSJ8NkvWQw72H4vRfGcmJNBav1UjjySf64G\nDaB69d0v7t/TTkqA/HyYPi/Axd+ksWCLn9rFPo6ZmsmvExOp1DiLlee0x42zmhsd0ohFGxaSnJTM\nd92/Y9HGRTtNf5a1E74CmAcOOwx++w0OP9zrSiTSfPcd3H23jYLtGOCnT7fz2qZNs+8vEZFQKyiw\nkzoGDbI1XRddBBd3C9B3dhrTV+9+PdfebN8NmZObw1HVk3m0YSZL5iYyfz7/XEuWWDueGjXseL+a\nNW19mVspwG9Hp7G5ip+D83wkT8okb30iGzfaCFsgALXbZLH8LAtaccTzeOMMrj0tlYNrljzvnkLX\njhTAIsixx9qZWscf73UlEkkKCmydxDPPWG+dXT34oL0gffxx6GsTkdi1dKn1/3r7bRuVuv562/V4\n8MH267tbNF9apfm927bZ8p61a+3Kz4e/1mdxz1/tKaKQisTzynE28paYCImJFtS2FO4ctHYMh/tT\nswJYBLnoIujWzbbjipTWa69Zw9UxY3Y/3L5li+2MfPNNOP300NcnIrHDdSEzE159FcaOtZ9pvXrZ\na1A42HH0bG8jbwcSDrdTAIsg999vc9YPPuh1JRIpNmyAZs1g9Gho3XrP9/vmG5uinDbN1hqKiART\nQQF89pmdxrFlC/Tubf2+qlXzurJ/C0a4Ko2yBjD1ZPdAs2Ywa5bXVUgkeeopOPfcvYcvsKlJnw8G\nDAhNXSISGwIBeOEFOOooG2X/v/+DnBwLYOEYvqDkyKRQtq/YHxoB88DkydbHaepUryuRSLB4MbRp\nA3/+WdL9fm8WLbJ1hhMmQIsW5V+fiESvdevg5ZdtqrFjR7j3XvUc3JVGwCJIq1bWCT8/3+tKJBI8\n8oidp1aa8AXQsKGdxXbdddb2RERkf+XmQt++NuK1cKG1hRg6VOErmBTAPFClCjRtaiMaInszbZq1\nnujTZ/9+3y23WFPWV14pn7pEJDqtXQsPPGDd6DdsgClT4J13bOmMBJcCmEfS0mwHicje9OkDDz9s\nmzb2R1yc9eJ5/HGYO7d8ahOR6BEIwGOPWdBat87e/L3+up3cIuVDAcwj7dtbozqRPRkzxsJTz55l\n+/1Nm9pO2+uv3/nMNRGR7fLy7Cizo46yzWGTJllPr9IueZCyUwDzSHq6BbCCAq8rkXBUXGyjX08+\naWelldUdd1iTwldfDV5tIhL5CgosaDVtaht2fvjBjg5q0sTrymKHAphH6ta1dxw//+x1JRKOhgyx\n4zQuvvjAHqdCBRg82KYWpk8PTm0iErmKiuDDD22H9BdfWHPnL76AlBSvK4s9CmAeOucca5wpsqPc\nDQHufTGL/k8Gdtvxfn8ddRT897/Q7ZoAGfOyCOQHDvxBRSTijBljLWreeMMW1o8aBSec4HVVsUt9\nwDz0559w4YUwb97uj5aR2BPID9D86TRWFPk5+jDffh1cuzcbtwao/6gdTNvq0OA9roiEP78f7rsP\nZs+Gp5+GLl30MyeY1AcsAh19tE0zTZ7sdSUSLibOzmZ5oR83rpCc3Bz8uf6gPK4/N5u8g/0UO4X4\nVwXvcUUkfK1caT0EO3aEM86wIHbRRQpf4UIBzEOOYweYfvSR15VIuPj23RRqFfuIj4snOSkZX1Jw\nTrZNqZOCL8lHReJxVidzeKUwOTFXRIIuLw+eeMKOJata1Rp/33nngW3okeDTFKTHFiyA44+HJUt0\neHKsmz/fvhcmTQ2wOi74B8huP5j23QE+1ixL5LPP9E5YJJoUF9sGnocesrVdTz1la0ClfJV1ClIB\nLAyceaadJn/FFV5XIl7q3t26Tz/6aPk+z9atcNJJcMMN0KtX+T6XiIRGZibcfbc1YX7uOTjlFK8r\nih0KYBFs+HAYMACysryuRLzy22/QubM1Qjz44PJ/vtmzLYSNGWMHfYtIZJo9G+6/H37/3foGXnaZ\nhTAJHS3Cj2CdO8OqVeoJFqtc13Yo9esXmvAF1nzxpZfsxTqgrhQiEWfNGlvX1a6dHZA9fbqtKVb4\nihz6pwoDFSrYf6Rnn/W6EvHC55/bi+l114X2ebt3hw4d7HljbNBZJGLl59vRQS1a2CkXOTl2eHaV\nKl5XJvtLU5BhYvNmOwJi9GhrTyGxYcsWSE6G996z46lCbetWO5e0a1cbhROR8BPID/DXymzmZqXQ\n78FEWraEZ56Bli29rkxAa8CiwnPP2TqwYcO8rkRCpX9/ewf76afe1bB4sU1hfPghdOrkXR0i8m+B\n/ADHvpbGnPV+Kgd8DD0rk/PPVBPlcKI1YFHgllssgP36q9eVSCgsXAivvGLvZL3UoIFtXb/iCqtJ\nRMLDggVw0c3ZzFnvhwqFFNXMIcmnJsrRQgEsjFStaiMi99yjNTmx4N574Y47oGFDryuxTtn33WeH\nf+fleV2NSGxbt87+Px53HJxwRAqtDg1+c2bxnqYgw0xRERxzjPWCuvhir6uR8vLjj3D99Tb9GC6L\nZ13XRsFc10bE1KRVJLS2bbODsv/7XzsnuH9/OOywkibKwW7OLMGhNWBRZPx4a8yakxO6tgQSOgUF\nFrIfe8wOxQ0neXk2Gnb22eXfEFZEjOvabugHHoBmzawvZEqK11VJaSmARZmrroK6ddWaIho98YR1\nrf7uu/AcZVqxAlJT7RiTyy/3uhqR6JaVZctOtmyx1/vTTvO6ItlfCmBRZuVKa0fxzTd2ppdEhxkz\n7IiQKVOgUSOvq9mzadPsB8GIERbGRCS45s6Fvn0tgD3+OFx5pfWElMijXZBRpm5d61R+zTXWq0ki\nX3Gxrfvq3z+8wxdY+H/nHbjoIu2MFAmmNWvgrrvgxBOhdWuYOdNe5xW+Yo8CWBi77DLw+eDBB72u\nRILhtddsyvGWW7yupHTOO8/WpJx1lv3QEJGyCwRspKt5c3tT7ffDQw/Z7neJTZqCDHNr1tiC7ddf\ntx+IEpkWLIDjj4eJE+0FOJL06QMTJsDYsfphIbK/tm6FgQNtTWWnTnbma9OmXlclwaQ1YFFswgRr\nSfHbb9Y0UyKL69ooUseONqIUaYqLbYpk40bbqVWxotcViYS/ggJ4/334v/8r2fWsY+aik9aARbFT\nTrFdMl26qElmJHrrLcjNtX/DSBQXB4MG2Tv5Xr3UJFhkb4qL4ZNPbPnIkCEwdCh89ZXCl/ybRsAi\nhOtCt25QqZK9qwrH9gXyb9t3PWZmRv7BuYGAjeKdc469qxeREq5ru9YffhgqV7Z2MzpbNTZoCjIG\nbN4MaWnQtattX5bwlp9vLRx69oSbb/a6muBYtQrat4cbbrCjlEQExo2zzVKbNtlC+86d9SY5lpQ1\ngGk1RwQ56CB7h9WunZ0feMUVXlcke3PffdZuomdPrysJnjp1bDF++/a2IL9XL68rEvHOL7/Af/4D\n8+bZqPDll6udhJSeAliEqVfPOqh37Ai1atnibgk/n3wC335rGyei7Z1w/foWwjp0sDcF11zjdUUi\noZWZaYvqZ8ywVhLXXQfx8V5XJZFGU5AR6uef4YILbFda+/ZeVyM7ysmxcDJmDLRp43U15WfGDDj1\nVGsY3LWr19WIlC/XtanGxx6z5sQPPmhn9laq5HVl4jXtgowxJ51koyyXXGK9pSQ8rFljwfiZZ6I7\nfAG0aAEjR0Lv3rbLSyQauS6MGmXrb2++Ga691rrX33CDwpccGI2ARbjRo+0MsaFDIT3d62pi27Zt\ncMYZdnbnM894XU3oTJliOyPfeMOOLhKJBkVF8MUXMGCAHZT98MNw6aVa4yX/pl2QMWzcOHtheOOd\nAIe3ySalTgqJCYlelxVTXNfeEefmwvDhsfciPXUqnH02vPKKpiMlsm3daq1+nn0WkpLg/vvh/POt\nH57I7iiAxbgfJgQ469M0imv7aVXXR2aPTIWwEHroIRuN/PFHSIzRv/Y//7RNIS++aOeYikSS9ett\nFPfll+3YsD59rIdftG2ikeDTGrAYV7VRNiT5KaaQv1bkkL3K73VJMeOll2DYMNudGqvhC6B1awuh\nd95pHcBFIsGCBXZKRZMmtrFkzBj4+mtb86XwJeVJASxKpNRJwVfHR3xcPJUDyTzTx8eWLV5XFf3+\n9z947jkLHklJXlfjvVatrEXFfffZEUwi4ch1ISPD1iwef7xNL06dalOPKSleVyexQlOQUSSQH8Cf\n66dJoo97bkskO9tGZho39rqy6BPID/D4m9kMeTGF8aMTadLE64rCy+zZtiGhVy8LYyLhID/fdo+/\n+KKdq3v77dZK4uCDva5MIpnWgMlOXBdee8161rz5prVGkOAI5AdoMSCNZQV+WtTyMflmrbfbnSVL\nLIRdeCH897+azhHvLFtmr4MDB1p7mDvugDPP1MJ6CQ6tAZOdOI71Zxoxwt7l3XabveOTA+O60Kt/\nNssK/FChkLkbc/Dnar3d7tSvb9M8Y8bYSFhRkdcVSSwpLralARddZNOKK1fajvHvv7cduwpf4jV9\nC0a5E0+EP/6w9gjHHWfrHKRstm61Y3dm/FSy3i45KRlfks/r0sJW7drwww+2uPnKK20KSKQ8rVoF\nTz8NTZtaC4kzz7TO9W+8AS1bel2dSAlNQcYI14WPPoK774Ybb7QDZCtX9rqqyLF8ub2Trl8f3nsP\niivaejtfkk/Tj6WwdSt0725b/b/4Ag45xOuKJFoE8gP8tTKbdTNT+GBQIt9/DxdfDD17WlNkTX1L\nedMaMCmV5cttajInx94Rqnv+vmVmWni46Sbrhq0X9LIpKoK77rJpoO++gwYNvK5IIt20mQHO+DiN\nlcV+Ejb6eKxxJjdenaiALyGlNWBSKocdZgd4P/GETad16wZLl3pdVXgqKoLHH7fO7gMH2qihwlfZ\nVahgPdOuvRZOPhn++svriiQSBQLwzjt24H37S7LJxdZjFtfKIe0iv8KXRAwFsBjVpQtMnw5HHQVH\nHw2PPmovbGJmz4aOHa2n1ZQpcO65XlcUHRzHml4OGACdOtn6MJF9KSqy/4tXXmkjpyNGWMPfuT+n\n0OpQrceUyKQpSGHhQptaGzsW+va1NWJVqnhdlTcKCuwokiefhEcegVtvjb1zHUPlp5/sDNMnnoDr\nr/e6Ggk3xcXw88/w6afWz7BevZJR+x2bHm/vf6j1mOIVrQGTAzZ1KvTrB7/9Zueg3XgjVK3qdVWh\n88MP1h+oXj1bH6fmquVv5kw47zzrFfbUUwq7sc51YfJkC12ffWabNS67zIJ6s2ZeVyeyewpgEjS/\n/25rnyZMsJ1EvXtD3bpeV1V+pk2z9V1//QXPP29Na7XWK3TWrIFLLrFzND/6KLbP04xFrmvT/EOH\n2lW5soWuyy6D5GSvqxPZNy3Cl6A59lhrFTBhgv1wbNHCdgFmZNiLZbT480+4/HLr1t6xo+0MvfBC\nha9Qq1ULRo2ykH/KKeCfEyBrcRaBfC1KjFb5+dYQtVcvW9PVvTtUrAhffWVrU/v3V/iS6KcRMNmn\ndetg8GDbCQhw1VVwxRXQqJG3dZVFUZH9sH/+eXuhv+MO+yGgs+C857rw5PMBHl2QBkl+fHV8ZPbQ\nMU/RYt06az/y1VfWod7ng86dbcS5RQuvqxMpO01BSrlzXcjKgg8+sPUZLVtac9IuXeCII7yubu/m\nzLEQ+d57NtJy++02xVGpkteVyY6yFmeR9k57iiikAvFkXpdBuwapXpclZeC6Nqo8ahR8842tLe3Y\n0QLXuedG97IGiS0KYBJS+fm2a3L4cNsSfuihduTHGWfYNJLXuyhd19Z2ffON7aBavtwC1/XXW9sN\nCU+B/ABp76bhz82hwppkLtuSyZuvJJKQ4HVlUhpr1tjrwqhRNspVqZK9Lpx9Npx2Wmxt6pHYoQAm\nnikqgl9/tRfdUaNsbdXRR1sQa9vW1pQdeWT5Hn5bWGiL6H/5xdaujR0L1arZC/8ll1jjT+2wiwzb\n2wo0rOKj942JLF9uzYPr1fO6MtlVQQFMmlTyf3/mTGuQesYZFryOOkprKiX6KYBJ2Ni82baST5xo\nwWzqVNiwwRbVNmtmV6NGcPjhdiUlWVjaXUAK5AfIXpWNLykFNz+RNWusc/+iRbBggR3ynJNjL/yN\nGkFqKrRrZ++2w31aVPatuNj6hL3xBnzyCaSleV1RbMvPt//bP/1k16RJduj19sB10kma1pfYowAm\nYZ5LjxIAABn8SURBVG3NGgtLM2fCrFkWoJYssTC1Zo114T/4YJu6jI+3q8AJsOKcNApr+CHXx8FD\nM6mVmMjhh0PDhna1aGGLeVu2VPuCaDZyJPToYYfJ33tv+Y6mSom8PBtV3h64fvvN/s916GDXKadA\njRpeVyniLQUwiWhFRbBxI2zdatMaBQXw55osLhvVnsLiQuLj4snokUFqfS3IjlWLFtk6vtq14f33\noWZNryuKLq5rp2L88kvJlZ0NrVqVBK6TT7bRahEpoQAmUWf7guyc3BySk5LVkkDYtg3uvx++/NK6\npbdt63VFkSsQsBGtHQNXXJxN4aem2nXccXDQQV5XKhLePAlgjuPUAD4FGgELgEtd192wm/stADYA\nxUCB67p7fNlUAJMd6Zw32Z0vvoCbb7YTDHr31kLvfVmzxtZi7ngtWgRt2ljQOvFE+9iggf4uRfaX\nVwHsaWCN67oDHMe5H6jhuu4Du7nfPOA413XXleIxFcBEZJ/mzoWuXW2n3dtva2oMbCp//nzw+3cO\nWxs2QOvWcMwxJVdysq21FJED41UAmwF0cF13peM4hwLjXdf9V09jx3HmA8e7rrumFI+pACYipbJ1\nK9x5J/z4ozUHbt3a64rKx/bdwCl1UkhMSKS4uCRo5eTYR7/fNrnUrm0bU9q0KQlbjRtr44JIefEq\ngK11Xbfmnj7f4evzgPVAEfCm67pv7eUxFcBEZL98+CHcdRc89RRcd110TKO5rk0d/jE9wHWZaSzd\n5icx30ejHzKZ40/8J2j5fDaapd3AIt4oawCrWIoHHgPseGiEA7jAw7u5+56S08mu6y53HCcJGOM4\nznTXdSfs6Tn79ev3z+309HTS09P3VaaIxLArr7SGv127wrhx8PrrkTElWVho7Vjmz7e+dvPn29Tq\n7Nl2ARx2QjZL2vlx4wrZVDmHOx7307VdqoKWiEfGjx/P+PHjD/hxDnQEbDqQvsMU5DjXdVvu4/c8\nCgRc131+D7+uETARKZMtW2wkbOxY+OgjW1jupc2bYfFiu5YssY/bg9aCBbBsGdSpYydFHHGEfWzc\n2JqbNmsGtWrBpm3aDSwSzrxchL/Wdd2n97QI33GcqkCc67qbHMc5CBgN9Hddd/QeHlMBTEQOyBdf\nwC23wB13WNuKYB9DtXUrrFwJK1bYOaMrVpTc3h60Fi+2+zVoAPXrl3zcHrSOOMKaCZemc7x2A4uE\nL68CWE1gKNAAWIi1oVjvOM5hwFuu657nOM6RwHBserIi8JHruk/t5TEVwETkgC1eDFddZevBPvjA\nws/eFBbamqvtwWrXa8egtWUL1K0Lhx1mB9HveO0YuGrWjI71aCKyZ2rEKiKygy1bLEzdcw8MHw4X\nX2yd3HNz7Vq9eufbGzbYsTp7ClY7fq1GDQUrETEKYCISFXZtuQB2KPf69f8OTXu7XVRkB70nJcGc\nOdb5HeDBB23qLynJWjZsv0+NGsGfqhSR6KcAJiIRZfNmm9bbfq1cCUtXB3i7KI21FfxU3eKjwZhM\n1i5PZO1aOxJne1jaMTjt7nbt2na4+46jVBs2QK9e8Mcf8PHHcPTR3v3ZRSR6KICJSFgoKLCF6AsX\n2i6/5ctLPu54bdtm03rbr7p1YWtSFu857SmmkIpOPO92yKBTs1Rq1w5O13bXtZ5hd98NDz8Mt92m\nBqUicmAUwEQkJLYHrAULdn+tWGGBqmFDOPzwnUPW9qtePahe/d/rqEJ1APvcudY7rFo1eO89q0lE\npCwUwEQkqDZvhhkzSo672X4tXmwL0Y84YvdX/foHNloVqpYLhYXw+OMwcKBdF15Ybk8lIlFMAUxE\nyqSoCKZPhylTIDu7JGitXGnNQJOTd74aNy5d76pIkZVlo2GnngovvGBrx0RESksBTERKZelSmDjR\ngsevv9qi9Hr14LjjbGH69nMFjzwydnYFBgJw++329/Lhh9C2rdcViUikUAATkX/s2Mph/apExo6F\nH3+ECRMsbJx0kvXEOuEEO0PxkEO8rjg8DBsGt95qi/P79o2dACoiZacAJiIArN4YoO3/0li42U/8\neh8HDc3k9PaJdOoEaWnQvLmaiO7NkiVwzTWQn28d9I880uuKRCSclTWAaQO2SBTYtAmGDIFLLoFG\nbbOZH/BT7BRSVDOHr3/x88kncOON0KKFwte+1K8PY8ZAly42FfnBB9a+QkQkmBTARCJUQQF89RVc\neqm1e/joIzj3XMj+IYXW9XzEx8Xjq5NMq7o+r0uNOHFxdoTR2LHw1FPQrRusW+d1VSISTTQFKRJh\n5s6Ft9+2/lVNmth02UUXQa1aJfcJVSuHWJCXB/ffb2H3/fchPd3rikQknGgNmEgUc13IzITnnoOf\nf4arr4YbboCWLb2uLHaMHGl/51deCY89Fl2tOESk7BTARKKQ68KoUfDoo3YY9V13WfiqWtXrymJT\nbq6FsMWLbcpXAVhEFMBEokxGBjzwgB0i3a8fXHyxzi0MB64Lb75pZ0n27w+33KKNDSKxTAFMJEos\nXAh9+sAvv8ATT8Dll6sfVTiaOROuuMKOZRo0yA4TF5H/b+/eo6ys6z2Ov3/AgBcG8jJcFMVLqbGJ\nPFAIKYhZKUpcXOUF83ggPVpSapaXann8Q1taq5Iycx0vpSneUgnvmslFxEuSyQwimEpeEUFwq4XA\n/M4fv+FIyQwDM/t59t7zfq31rL0ZHmd/WY/PzGf9fr/n++t4bEMhVbjGRpg6NXWkHzAgbQ90/PGG\nr3K1775pPd6gQbD//nDrjCLzXp5HcU0x79IkVQBHwKQysHRpWtwdI1xzTdqDUZXjvj8VGTt9BOt3\nbKDQq8Dcr8/x6VOpg3AETKpQd9+dGn6OHQuzZhm+KlHPT9QT61Lz2wVvLOSePzfkXZKkMmcAk3IS\nI1xyCZx6Ktx+O3zve043VqqBvQZSqEvNb/t1G8A3v1pg2rS8q5JUzpyClHLQ2AhTpsDcuXDPPamT\nvSrbxs1v//ZsLUcfnZq2Tp0K226bd3WSSsWnIKUK0diY9mV8/nm4807o0SPvilQKxSKccgrU18Mt\nt6R9OCVVH9eASRXi9NNh8eK09svwVb1qa1Oz1m99C0aMSO8laQNHwKQM/epXcPnlqX1Bz555V6Os\nPPMMfPWr8PnPpylJtzGSqodTkFKZmzcPxo9Pr3vtlXc1yto776RWI6tWwa232rhVqhZOQUplrFhM\nv3yvuMLw1VH16AHTp6eF+UOHwlNP5V2RpDw5AiZl4Oyz4Y034Lrr8q5E5eC221L7kalTYeLEvKuR\n1BZOQUplaskSGD48PQ3Xp0/e1ahcLFgA48alTdYvvtgecFKlMoBJZWryZNhjDzj//LwrUblZsQKO\nOQZqauDmm30qVqpEBjCpDL3+OhQKqefXjjvmXY3K0bp1qSnvY4+lpry77JJ3RZK2hIvwpTJ09fVF\nDjxmHjXbF/MuRWWqSxf49a/h2GM/nKqWVP0cAZNKpLimSO/zRvBBzwYG9i4wZ9IcarvV5l2Wyti0\naXDmmfCHP8CwYXlXI6k1HAGTyszcJfX8o3sD61nHwuULaVjekHdJKnMTJ8Jvfwtjx8If/5h3NZJK\nyQAmlciqJQPp/s8CNZ1qGFA3gEJdIe+SVAFGj05tKiZOTHuFSqpOTkFKJfKjH8Gbq4oc++0GCnUF\npx+1RZ58Eo48Eq65BsaMybsaSc1xClIqM0uWwMB9ahnWb5jhS1vss59NI2CTJ6enIyVVFwOYVCIr\nV0JdXd5VqJIdcEBakH/iiTB3bt7VSGpPBjCpRIpF6N497ypU6YYPh+uvh6OOggaf45CqhgFMKpEu\nXVKTTamtDjsMfv5zOOIIWLYs72oktQcDmFQi228P776bdxWqFhMnwqRJMGECrFmTdzWS2soAJpXI\nrrvCK6/kXYWqyfnnp62KpkzJuxJJbWUAk0pkzz3hb3/LuwpVk06dUqPW2bPhxhvzrkZSWxjApBIZ\nPBjmz8+7ClWb7t3hppvg9NMN+FIlsxGrVCKrV0O/frB8OWyzTd7VqNr89Kdw993w0EMQtrgFpKT2\nYiNWqcz07AmDBsGsWXlXomp0xhnpIY+rr867EklbwwAmldCYMTB9et5VqBp17gxXXQXf/z6sWpV3\nNZK2lFOQUgktXQpDhsCrr0K3bnlXo2p08smw/Q5FjvlWPQN7DXTbKyljTkFKZah/f9h/f7j11rwr\nUbX63g+K/PL9EYz87UhG/GYExTXFvEuS1AoGMKnEzjwzLZh2YFelsKJzPezcwLrGdSxcvpCG5e5X\nJFUCA5hUYqNHw9q16Yk1qb0N7DWQfXYowPoa9ttpAIW6Qt4lSWoFA5hUYp06wUUXwXnnwfr1eVej\nalPbrZYnTp3D4L/O5uyd57gGTKoQBjApA2PHwsc+BldemXclqka13Wo5bdww7rjZ8CVVCp+ClDLy\nzDNw6KFQXw+9e+ddjarNsmWw336p8W+XLnlXI3UcPgUplblBg+Ckk+DUU12Qr/bXu3faeeHpp/Ou\nRFJrGMCkDF1wAbz4ot3LVRqDB6eRVknlz4FqKUPdusENN8CoUTB0aBoVk9rLJz8JCxfmXYWk1nAE\nTMpYoQBTp8KECfD223lXo2rSp09aAyap/BnApBxMnAjjx8NRR8GaNXlXo2rRsyesXp13FZJawwAm\n5eTHP4YddoDJk6GxMe9qVA3WroWamryrkNQaBjApJ507p/Vgf/87TJnik5Fqu2IRunfPuwpJrWEA\nk3K07bZpi6I//xnOOssQprZ58UXYY4+8q5DUGgYwKWc9esB998Ejj8Bppzkdqa3X0JCasUoqf3bC\nl8rEO+/AmDHQv3/qE9a1a94VqZI0NkJdXeoDtuuueVcjdRx2wpcq3IaRsGIRRo+GVavyrkiV5Ikn\nUjd8w5dUGQxgUhnZbju47TYYOBCGD4dFi/KuSJXi2mvh+OPzrkJSa7UpgIUQvhJCqA8hrA8hDG7h\nvMNDCItCCItDCOe05TOlate5c2rU+t3vwogRcMsteVekcrdyZfr/5MQT865EUmu1dQRsATABmNXc\nCSGETsBlwGFAATguhOAyUWkzvv51uP9+OPdcOPPM1ONJ2pRLL007K/Trl3clklqrTQEsxvhcjHEJ\n0NLis6HAkhjj0hjjWuAmYFxbPlfqKAYPTi0qliyBQw6Bl1/OuyKVm5degssvhx/+MO9KJG2JLNaA\n7Qps/GvjlaavSWqFHXeEGTPSE5KDB8NVV9kvTEmM8M1vwhln2P9LqjRdNndCCOFBoPfGXwIi8IMY\n452lKOqCCy74//ejRo1i1KhRpfgYqWJ06pSmIo88EiZNSut9rroKdt8978qUp1/8At56C84+O+9K\npI5j5syZzJw5s83fp136gIUQHgbOijHO38TfDQMuiDEe3vTnc4EYY7ykme9lHzCpBevWwU9+Aj/7\nGVx0EZx8MoQt7kCjSvfww3DMMfDYY7DXXnlXI3Vc5dAHrLkPfxL4eAihfwihK3AsMKMdP1fqULp0\ngfPOg5kzU8PWz30OHn8cimuKzHt5HsU1xbxLVAkV1xS5YfY8jj6hyE03Gb6kStWmEbAQwnjgl8DO\nwCrg6Rjj6BBCX+DKGOOYpvMOB6aSAt/VMcaLW/iejoBJrdTYCNddB+ddUOSfx43g3e0aKNQVmDNp\nDrXdavMuT+2suKbIZ389gudWNrD7dgXqz/Q6S3nb2hEwtyKSqsBDi+fxpWkjaQzr6EQNDx0/m1Ef\nH5Z3WWpnV943j/9+dCR0XkdNpxpmT5rNsH5eZylP5TAFKSknQ/sP5FN9CnTpVEPtPwdw3BcKXHop\nvP9+3pWpvfz+93Du5IHs2b1ATacaBtQNoFBXyLssSVvJETCpShTXFGlYnqYgn19Yy4UXwqOPwne+\nA9/4BnTvnneF2hpr16YnYG+7DW6/HT5R+PA6O/0o5c8pSEkfsWABXHhhWrB/xhmpZ1TPnnlXpdZ6\n9tm0vVBdXVrrt9NOeVck6d85BSnpIz71Kbj55tSyoL4e9twTTjsNFi7MuzK15IMP4OKLYeTItCXV\nXXcZvqRqYwCTOoABA+CGG9KI2M47w6GHpmP69NRXTOXjgQfg05+GRx6BJ56AU06xz5tUjZyClDqg\nDz5Ia4ouuwxeeQVOOgm+9rU0QqZ8zJ+f1nq99FJqtDt2rMFLqgROQUpqta5d4bjjYO5cuOMOWLYM\nhg6Fgw6CK66AFSvyrrDjmDcv7fM5ZgxMmAANDTBunOFLqnaOgEkC0tN2998P118P994LhxwCxx8P\no0f7BGV7W7cuBd+pU9MI5DnnpD0+t9km78okbSmfgpTUblavTi0PbrwxjdAcdBB8+cvp2G23vKur\nXAsWwO9+l9bj7b03fPvbMH582l5KUmUygEkqidWr08jYnXemkbHddktB7Igj4DOfMTxszmuvwbRp\nKXi9/XYaVTzhhPRghKTKZwCTVHLr1qURsRkzUihbujRtBn7wwTBqFAwZAjU1eVeZr/Xr4amnUli9\n91547jk46qgUukaOhE6uvJWqigFMUubeegvmzEmNXmfNghdegOHD4cADUxgbMgT69Mm7ytKKMf27\nH30U7rsvtZHo1SutnRs9Ok3fduuWd5WSSsUAJil3K1akQPbYY2kUaP78FD42hLEhQ2DQoDSNWc4j\nQcU1RerfrGdgr4Ef2e5n9Wp48sn0b3z88fTatWsKnl/8Ihx+OPTvn1PhkjJnAJNUdmJM05QbwthT\nT6U2CytWpEXo++zzr8eee0Lv3tC5c341F9cUOeg3I1i4vIHdty1w2jZzeGlxLYsWwaJFqfYhQ+CA\nA2DYsPTar19+9UrKlwFMUsV47z14/nlYvPhfjxdegJUr096Hu+zy4dG3bwpmtbXQo0d63fj99tun\nvlmbOmJMn/fuux893n4bXn89Ha+9ll5fWDuPNw4fCZ3XEdbXMHblbA7eexj77gv77ZdGt/IMiJLK\niwFMUlVYuzY1ht0QijYcb74JxSK880563fj9e++loLWpI4TUx2xTR8+eKdxtCHl9+0KPuiIn/GkE\ni1YsZEDdAOZMmvORaUhJ2sAAJkntpLimSMPyBgp1BcOXpBYZwCRJkjLmXpCSJEkVwgAmSZKUMQOY\nJElSxgxgkiRJGTOASZIkZcwAJkmSlDEDmCRJUsYMYJIkSRkzgEmSJGXMACZJkpQxA5gkSVLGDGCS\nJEkZM4BJkiRlzAAmSZKUMQOYJElSxgxgkiRJGTOASZIkZcwAJkmSlDEDmCRJUsYMYJIkSRkzgEmS\nJGXMACZJkpQxA5gkSVLGDGCSJEkZM4BJkiRlzAAmSZKUMQOYJElSxgxgkiRJGTOASZIkZcwAJkmS\nlDEDmCRJUsYMYJIkSRkzgEmSJGXMACZJkpQxA5gkSVLGDGCSJEkZM4BJkiRlzAAmSZKUMQOYJElS\nxgxgkiRJGTOASZIkZcwAJkmSlDEDmCRJUsYMYJIkSRkzgEmSJGXMACZJkpQxA5gkSVLGDGCSJEkZ\na1MACyF8JYRQH0JYH0IY3MJ5L4UQ/hpC+EsI4Ym2fKbK18yZM/MuQW3g9atcXrvK5vXrmNo6ArYA\nmADM2sx5jcCoGON/xBiHtvEzVab8IVLZvH6Vy2tX2bx+HVOXtvzHMcbnAEIIYTOnBpzulCRJArIL\nRRF4MITwZAjh5Iw+U5IkqSyFGGPLJ4TwINB74y+RAtUPYox3Np3zMHBWjHF+M9+jb4zx9RBCHfAg\nMCXG+Egz57ZckCRJUhmJMW5uJvAjNjsFGWP84taV8y/f4/Wm1+UhhDuAocAmA9jW/CMkSZIqSXtO\nQW4yOIUQtgshdG96vz3wJaC+HT9XkiSporS1DcX4EMLLwDDgrhDCvU1f7xtCuKvptN7AIyGEvwCP\nAXfGGB9oy+dKkiRVss2uAZMkSVL7yrU1hI1cK9sWXL/DQwiLQgiLQwjnZFmjmhdC2CGE8EAI4bkQ\nwv0hhJ7NnOf9VyZacy+FEH4RQlgSQng6hLB/1jWqeZu7fiGEg0MIq0II85uOH+ZRpz4qhHB1CGFZ\nCOGZFs7Zonsv795cNnKtbJu9fiGETsBlwGFAATguhLBfNuVpM84F/hhj3Bf4E3BeM+d5/5WB1txL\nIYTRwN4xxk8ApwBXZF6oNmkLfhbOjjEObjouzLRIteQ3pGu3SVtz7+UawGKMz8UYl9DMAv6N2Mi1\nDLXy+g0FlsQYl8YY1wI3AeMyKVCbMw64tun9tcD4Zs7z/isPrbmXxgHXAcQYHwd6hhB6o3LQ2p+F\ndgIoQ02ts95u4ZQtvvcq5YeqjVwr167Ayxv9+ZWmryl/vWKMywBijG8AvZo5z/uvPLTmXvr3c17d\nxDnKR2t/Fg5vmsK6O4QwIJvS1A62+N5r01ZErdGaRq6tcODGjVxDCM8218hV7audrp9y0sL129Ta\nkuaeyPH+k7LxFLB7jPH9pimt6cA+OdekEil5AMu6kavaVztcv1eB3Tf6c7+mrykDLV2/pgWlvWOM\ny0IIfYA3m/ke3n/loTX30qvAbps5R/nY7PWLMb670ft7QwiXhxB2jDGuzKhGbb0tvvfKaQrSRq6V\nrbl1C08CHw8h9A8hdAWOBWZkV5ZaMAP4r6b3JwJ/+PcTvP/KSmvupRnAfwKEEIYBqzZMMyt3m71+\nG68ZCiEMJbWKMnyVj0Dzv+u2+N4r+QhYS0II44FfAjuTGrk+HWMcHULoC1wZYxxDmj65o2mPyC7A\nDTZyLQ+tuX4xxvUhhCnAA6TAf3WM8dkcy9aHLgFuCSFMBpYCR0NqpIz3X9lp7l4KIZyS/jr+b4zx\nnhDCESGE54H3gEl51qwPteb6AV8JIXwDWAv8Azgmv4q1sRDCNGAUsFMI4e/A/wBdacO9ZyNWSZKk\njJXTFKQkSVKHYACTJEnKmAFMkiQpYwYwSZKkjBnAJEmSMmYAkyRJypgBTJIkKWP/By8TkWNHbjaW\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107b65048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "plt.plot(yp[0], yp[1], '-')\n",
    "plt.plot(res.y[0], res.y[1], '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that at some regions the solver takes quite large steps, but an accurate solution is available at all points nevertheless."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Bouncing ball"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example illustrates event detection capabilities. The problem is very simple: a ball bounces from the ground and loses some fraction of its vertical velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g = 9.81  # Gravity.\n",
    "alpha = 0.8  # Ratio of preserved velocity.\n",
    "v_hor = 1  # Horizontal velocity.\n",
    "n_bounces = 10  # Number of bounces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the system right-hand side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(x, y):\n",
    "    return [v_hor, y[2], -g]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the bounce event:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bounce(x, y):\n",
    "    return y[1]\n",
    "\n",
    "bounce.terminate = True  # This event is terminating.\n",
    "bounce.direction = -1  # Detected when y[1] is decreasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver runs until the bounce event is detected and then gets restarted with a new value of the vertical velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ya = [0, 0, 10]  # Start at (0, 0) and with the vertical velocity of 10.\n",
    "\n",
    "sol_points = []\n",
    "sol_smooth = []\n",
    "\n",
    "for i in range(n_bounces):\n",
    "    # Use a 3-rd order Runge-Kutta method.\n",
    "    # Integration interval is not relevant. \n",
    "    res = solve_ivp(fun, [0, 100], ya, events=bounce, method='RK23')\n",
    "    sol_points.append(res.y)\n",
    "    t = np.linspace(res.x[0], res.x[-1], 100)\n",
    "    sol_smooth.append(res.sol(t))\n",
    "    ya = res.y[:, -1].copy()\n",
    "    ya[2] *= -alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sol_points = np.hstack(sol_points)\n",
    "sol_smooth = np.hstack(sol_smooth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x10814f940>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAFwCAYAAACRj46qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX6P/DPGULL0Jv0oDSpUkVxxQCiiAUrirprdNlV\n1x/oCrrqqkl21+/qipVVdq1BFxtYwQaoASkCSo00ERI6Kj0RKZnn98czgYAkmXJm7r0zn/frlRcw\nTO595gYynzzn3HOMiICIiIiIouNzugAiIiKiRMBQRURERGQBQxURERGRBQxVRERERBYwVBERERFZ\nwFBFREREZIGVUGWMqW2MmWSMWWmM+dYY08fGcYmIiIi8IsXScZ4C8JGIXGWMSQGQaum4RERERJ5g\nol380xhTC8BiEWltpyQiIiIi77Ex/HcygJ+MMS8bYxYZY54zxlS3cFwiIiIiz7ARqlIA9ADwjIj0\nAPAzgHssHJeIiIjIM2zMqdoEYKOIfB3882QAfzn+ScYYbjJIREREniEiJpznRx2qRGS7MWajMaad\niKwBMBDAijKeG+3pyCFZWVnIyspyugyKAL923savn3fxa+dtxoSVpwDYu/tvFICJxpjKANYBuNHS\ncYmIiIg8wUqoEpGlAHrbOBYRERGRF3FFdQpJenq60yVQhPi18zZ+/byLX7vkE/U6VSGfyBjhnCoi\nIiLyAmNM2BPV2akiIiIisoChioiIiMgChioiIiIiCxiqiIiIiCxgqCIiIiKygKGKiIiIyAKGKiIi\nIiILGKqIiIiILGCoIiIiIrKAoYqIiIjIAoYqIiIiIgsYqoiIiIgsYKgiIiIisoChioiIiMgChioi\nIiIiCxiqiIiIiCxgqCIiIiKygKGKiIiIyAKGKiIiIiILGKqIiIiILGCoIiIiIrKAoYqIiIjIAoYq\nIiIiIgsYqoiIiIgsYKgiIiIisoChioiIiMgChioiIiIiCxiqiIiIiCxgqCIiIiKygKGKiIiIyAKG\nKiIiIiILGKqIiIiILGCoIiIiIrKAoYqIiIjIAoYqIiIiIgsYqoiIiIgsSHG6AK8LBIA1a4ANG4DK\nlYGOHYGTTnK6KiIiIoo3K6HKGJMPYA+AAIBDInK6jeO62e7dwBNPAC+8AFStCrRuDRw4AOTlAe3a\nAaNHA1dcAfjYCyQiIkoKRkSiP4gx6wD0FJFd5TxHbJzLDT76CBgxAjj/fOCuu7Q7VeLwYeDjj4Gs\nLKB2beCVV4DmzR0rlYiIiCJgjIGImHA+x1YfxVg8lqs9/jjwxz8Cr78OvPzysYEKAFJSgIsvBhYs\nAAYOBPr0Ab75xplaiYiIKH5sdqp2AygG8JyIPH+C53i+U/Xww8CECcC0aUCLFqF9zjvvALfcAnzy\nCdCjR2zrIyIiIjsi6VTZmqh+lohsNcY0BDDdGLNSRGZbOrYrTJgAPP888OWXQNOmoX/e5Zfrrxdd\nBHz1FdCyZWzqIyIiImdZCVUisjX464/GmHcBnA7gV6EqKyvryO/T09ORnp5u4/Qxt2CBzp2aOTO8\nQFXi8suBggJg6FBg7lygenX7NRIREVHkcnNzkZubG9Uxoh7+M8akAvCJSKExxg9gGoBsEZl23PM8\nOfy3dy/QvTvw6KNHu06REAGuuQZo1AgYN85efURERGRfJMN/NkLVyQDeBSDQztdEEXn4BM/zZKi6\n5RaguFiH/qK1ezfQtSuQkwMMGBD98YiIiCg2HAlVIZ/Ig6Fq/nzg0kuBlSuBOnXsHPO994B77wWW\nLgWqVLFzTCIiIrLLySUVEk5xMfCnP+mwn61ABei8qtatdWkGIiIiShzsVJXhmWeASZOAL74ATFg5\ntWLr1gG9ewOLFgFpaXaPTURERNHj8J8lhYXaTZoxA+jSJTbnePBB3S8wJyc2xyciIqLIMVRZ8vDD\nOufp9ddjd47du4E2bXTtqjZtYnceIiIiCh9DlQX79mmXKjf311vQ2Pa3vwHff68LixIREZF7MFRZ\n8M9/AsuXA6+9Fvtz7dmjXaq5c4G2bWN/PiIiIgoNQ1WUioqAk0/WldM7dIjPOf/+d+1WcW4VERGR\nezi5919CePVVoG/f+AUqALjtNh1u3L4dOOmk+J2XiIiI7OI6VUEiun3MqFHxPW+9esBVVwHPPRff\n8xIREZFdDFVBn3+u61H17x//c48cCYwfDxw8GP9zExERkR0MVUHjxmm4sb3QZyi6dAHatwfeeSf+\n5yYiIiI7OFEdwPr1usJ5QQHg9ztTw7vv6pY4c+c6c34iIiI6inv/Reill4Drr3cuUAHAxRfrCut5\nec7VQERERJFL+lAVCACvvAJkZDhbR0oK8NvfciFQIiIir0r6UJWbC9StC3Tr5nQlwA03ABMnAocP\nO10JERERhSvpQ9WECRpm3ODUU4GWLYHp052uhIiIiMKV1BPV9+0DWrQAVq92z8Kb48dr9+zNN52u\nhIiIKHlxonqY3n4b6NfPPYEKAK65Bvj0U2DXLqcrISIionAkdaiaOFEnh7tJ3brAoEFcs4qIiMhr\nkjZU/fQTsGABcOGFTlfya1ddBUya5HQVREREFI6kDVXvvQecdx6Qmup0Jb82ZAgwbx6wc6fTlRAR\nEVGokjZUTZqkHSE3qlEDOPdc4P33na6EiIiIQpWUoWrHDu0EDRnidCVlu/JKYPJkp6sgIiKiUCVl\nqHr/fZ0MXqOG05WU7aKLgC+/BHbvdroSIiIiCkVShqrJk7UT5GY1awIDBgAffOB0JURERBSKpAtV\ne/dqB+iii5yupGJXXMGlFYiIiLwi6ULV9OlA377aCXK7Cy4AvvgC+OUXpyshIiKiiiRdqJo61Rtd\nKgBo0ADo0gWYOdPpSoiIiKgiSRWqAgHgo4/cueBnWS68EPjwQ6erICIioookVaj6+mvt/pxyitOV\nhO6ii7S75rK9qImIiOg4SRWqvDT0V6JzZ6C4GFi50ulKiIiIqDwMVS5nzNFuFREREblX0oSqzZuB\nggLgzDOdriR8DFVERETulzShato03U8vJcXpSsKXng4sXqxrbBEREZE7JU2omj5dt6bxourVgT59\ngNxcpyshIiKisiRFqAoEgBkzvBuqAK19+nSnqyAiIqKyJEWoWr4cqFMHSEtzupLIDRqkwZCIiIjc\nKSlC1fTpOp/Ky7p1A376Cdi0yelKiIiI6ESSJlR5eegPAHw+YOBADgESERG5lbVQZYzxGWMWGWM+\nsHVMG375BZg7F+jf3+lKosd5VURERO5ls1N1O4AVFo9nxdy5uip5nTpOVxK9c88FPvtMJ94TERGR\nu1gJVcaY5gCGAHjBxvFsmjFDh80SQVoaULu2TrwnIiIid7HVqXoCwF0AXLftb25uYgz9lUhPB2bO\ndLoKIiIiOl7UocoYcyGA7SKyBIAJfrhCURGwdKk3t6YpyznnMFQRERG5kY1NW84CcIkxZgiA6gBq\nGmNeEZHfHf/ErKysI79PT09Henq6hdOXbd48XYogNTWmp4mrc84B7rgDENHNlomIiCh6ubm5yI1y\n6xIjYm/EzhhzDoDRInLJCf5ObJ4rFA88ABQXA//3f3E9bcy1bg188AHQqZPTlRARESUmYwxEJKz2\nRUKvUzVrlnZ2Ek2/fvraiIiIyD2shioRmXmiLpUTfvkF+OYboG9fpyuxj/OqiIiI3CdhO1Xz5+vw\nWM2aTldiX0moivNoKhEREZUjYUPVzJmJOfQHAK1aAZUrA99953QlREREVIKhyoOM4RAgERGR2yRk\nqDp8GFiwIDHnU5U4+2xgzhynqyAiIqISCRmqli3TLV3q1nW6ktjp21f3NSQiIiJ3SMhQNXduYnep\nAKBjR+CHH4Aff3S6EiIiIgIYqjzL5wPOOENXjSciIiLnMVR5GIcAiYiI3CPhQtXmzUBhIdC2rdOV\nxB5DFRERkXskXKiaN0/DRjJsNnz66cCiRcDBg05XQkRERAkXqpJl6A8AatUC2rQBlixxuhIiIiJi\nqPI4DgESERG5Q0KFqv37geXLgV69nK4kfvr25SKgREREbpBQoWrxYuDUU4HUVKcriZ8zztDNo4mI\niMhZCRWqFiwA+vRxuor4at1a73bcts3pSoiIiJJbQoWqhQv1jrhkYoy+5oULna6EiIgouSVUqFqw\nAOjd2+kq4q93b33tRERE5JyECVU7dwLbt+ucqmRz+ukMVURERE5LmFC1cCHQsydQqZLTlcRf7976\n+kWcroSIiCh5JVSoSrb5VCUaNwZq1gTWrnW6EiIiouSVMKEqWedTlSjpVhEREZEzEiJUiWioStZO\nFcB5VURERE5LiFC1caP+2qKFs3U4iaGKiIjIWQkRqkq6VMY4XYlzevYEli4FDh1yuhIiIqLklBCh\n6ptvNFQks5o1gbQ0YMUKpyshIiJKTgkRqhYvBnr0cLoK5/XoASxa5HQVREREycnzoUpEg0T37k5X\n4rwePbRrR0RERPHn+VC1ZYsGq2bNnK7EeT17slNFRETkFM+HqpKhv2SepF6iWzdg2TKguNjpSoiI\niJJPQoQqDv2p2rWBJk2A1audroSIiCj5pDhdQLQWLQKuucbpKtyjZLJ6x45OV+KsPXuAGTOAvDzg\nwAGgeXMgPZ3XhYiIYichOlW88++oZL8DcNcu4M9/Blq1Al54QdftSk3Va3LeeUDfvsCXXzpdJRER\nJSJPd6p27tSP1q2drsQ9evQAHnrI6SriqKAAyMkBAgFs2urD1R9loPOFaVixQodCSysuBt56Cxg+\nHLjxRiArC6hUyYmiiYgoEXm6U7V4MXDaaYDP06/Cru7d9boEAk5XEgcFBcC4ccCYMZjcJRu/eXcM\ncnqOw3/vK/hVoAI0QA0frstOzJ0LXHKJDg0SERHZ4N04UlCAao9k474DmUB2tr7BEho0AOrUAb7/\n3ulK4iAnB8jOxpTP/bjtNuD9GX60fS1bHy/HSScBn36qw4JXX82tfYiIyA5vhqpgh+LFumOw/U/Z\nwJgx2rFgsAKQRPOqAgEs+NaP3/8emDpVu5bw+0Nq06WkABMnAocPA7fcEvtSiYgo8XkzVAU7FF8t\n9+tyCn6/dqsq6FAki27ddHPlRPfzLz7ccGUR/vtfoHfv4INFRSGPB1epArzxhg4Fvvxy7OokIqLk\n4M1QFQigCH7k55e6RT7EDkUy6NpVFwFNZCLA//s6A882ysRl5xXpg0VFQGYmkJER8nFq1AAmTwbu\nvluXXyAiIoqUN+/+8/mQN78IHTr4Ubly8LEwOhSJLhlC1cSJwKIdafjv5JHA2LEaqH0+YORIIC0t\nrGN16gT84x/AiBHAnDm8I5CIiCJjRCS6AxhTFcAsAFWgIW2yiGSf4HkS7bmOKCjAkt+Pw/PNs/FM\njv9ohyKCN9REFAjo6uobNgB16zpdjX07dmgQ+uAD4PTT7RwzEAD69dO7A2+7zc4xiYjIu4wxEJGw\nNsGLOlQFT5wqIj8bYyoBmANglIgsOO459kIVgLuvLsCVhTk4vVewQ5GRwUBVSt++wMMPa1BINCNG\n6J17Tz9t97grVuj1Wr7812tcERFRcokkVFkZ/hORn4O/rRo8pr30VIbP1qbhsn9nAmfG+kze1LWr\nTlZPtFC1fDkwZQrw3Xf2j92xo2bz7GzgP/+xf3wiIkpsViYhGWN8xpjFALYBmC4iC20ctyyHDwMr\nVwJdusTyLN6WqPOq7rsPuPdeoFat2Bz/3nt14nosQhsRESU2W52qAIDuxphaAN4zxnQUkRXHPy8r\nK+vI79PT05Genh7R+dau1eGZGjUiqzcZdO0KTJjgdBV2zZ6tnarJk2N3jvr1gTvvBO6/H3jzzdid\nh4iI3CU3Nxe5ublRHcPKnKpjDmjMAwCKROTx4x63Nqdq0iS9++u996wcLiHt2QM0a6a/JsrdbAMG\nAL/7XVgrJkSkqAho0waYNo3dUCKiZBXJnKqoh/+MMQ2MMbWDv68OYBCAVdEetzx5eXyzq0jt2kDD\nhsC6dU5XYsfChbr1znXXxf5cfj8wahTw6KOxPxcRESUOG3OqmgD4whizBMB8AJ+KyEcWjlum5cuB\nzp1jeYbEUDJZPRE88ggwejSOrksWY7feCnz4IXc+IiKi0EUdqkRkuYj0EJFuItJVRB6yUVh5li9n\npyoUiTJZffVqYNYs4Pe/j98569TR8z3+eMXPJSIiAjy4TU1REbB5M9C2rdOVuF+ihKqnn9ZNj/3+\n+J73jjuAV17ReWlEREQV8VyoWrkSaNcufsNAXpYIoWrfPuD114Gbb47/uZs2Bc47D3j11fifm4iI\nvMdzoYpDf6Fr0wbYvh3Yu9fpSiI3cSLQv7/eyeiEP/0JePZZ3cCZiIioPJ4MVZykHppKlXSPvOXL\nna4kMiIaaG691bka+vUDjNE5XUREROXxZKhipyp0Xh4CnDsX+OUXXZ/KKcZoqBs/3rkaiIjIGzwX\nqrhGVXi8HKpycvQOPJ/D/0qvuw745BNg1y5n6yAiInfzVKj66SftXDRv7nQl3uHVtar27wfefhu4\n/nqnKwHq1tUJ62+95XQlRETkZp4KVSXzqUxYi8Ynt86dgW+/9d5E6/feA3r3dm6C+vFuuCHx9lIk\nIiK7PBmqKHQNGgDVqgFbtjhdSXgmTNB9/tzivPN0y581a5yuhIiI3MpToYrzqSLTsSOwYoXTVYRu\n61Zg/nzgssucruSoypWBa6/lmlVERFQ2T4Uq3vkXGa+FqrffBi6+GEhNdbqSY117rc6r8tpQKhER\nxYdnQlUgoHODOPwXPq+FqrfeAoYNc7qKX+vZEzh0yLt3UxIRUWx5JlQVFAC1a+udWBSejh01kHrB\n5s06zDtokNOV/JoxGvZ4FyAREZ2IZ0JVXh67VJHq1Ek7VV4Ytnr7beCSS4CqVZ2u5MRKQpUXriUR\nEcWXZ0LVihXacaHwNWyoC2hu3+50JRVz69Bfie7dNVAtWeJ0JURE5DaeCVUrVwIdOjhdhTcZ4415\nVdu26TDluec6XUnZjAGuvFI7akRERKUxVCUJL4SqDz/U9aCqVHG6kvINHQq8/77TVRARkdt4IlSJ\nAKtWMVRFo1Mn909WnzJFl1Jwuz59gB9/1MVAiYiISngiVG3ZoquC16vndCXe5fZO1S+/AF98AVxw\ngdOVVMznAy66CPjgA6crISIiN/FEqOLQX/TcHqo+/xw47TSgfn2nKwkNhwCJiOh4DFVJonFj4PBh\nHbZyI68M/ZUYOBD45htg506nKyEiIrdgqEoSbr4DUASYOtVboSo1FUhPBz791OlKiIjILRiqkohb\nJ6svWaJz5tq3d7qS8Jx/PkMVEREd5ZlQxYU/o+fWTlXJ0J8xTlcSnvPPB6ZN4+rqRESkXB+qdu4E\n9u8HmjZ1uhLvc2uo+vBDvZvOa9q0AapXB5Yvd7oSIiJyA9eHqpUrgVNP9V4Xw43cGKp27dKazjrL\n6UoiwyFAIiIq4YlQxflUdjRrpl2/HTucruSoL77QQOXWDZQrMngwQxURESmGqiTixjsAZ8wABg1y\nuorI9e8PzJ8PFBU5XQkRETmNoSrJtG8PrF7tdBVHTZ/u7g2UK1KzJtCzJ5Cb63QlRETkNIaqJOOm\nUJWfD+zdC3Tp4nQl0eG8KiIiAlweqn7+Gdi2DTj5ZKcrSRxuClWffaYrk/tc/a+wYgxVREQEuDxU\nrV6tt62npDhdSeJwU6iaPt3b86lKdOumk/83bXK6EiIicpKrQxWH/uxr0wYoKAAOHXK2jkBAO1Ve\nnk9VwucDzjkHmDnT6UqIiMhJDFVJpmpVXVph3Tpn61i6FKhfH2jRwtk6bElP52R1IqJkx1CVhNww\nBDhjhs6nShQMVURE5OpQtXq1rqZOdrkhVM2apWs8JYpOnXR1+M2bna6EiIic4tpQFQgA338PtG3r\ndCWJx+lQFQgAs2cDZ5/tXA22+XxAv36cV0VElMxcG6o2bgTq1QP8fqcrSTxOh6q8PKBRI+Ckk5yr\nIRY4BEhElNxcu1jBmjVAu3ZOV5GYnA5Vs2YlVpeqRHo68O9/O12FXSK6rdHKlfrnZs2A3r25zAkR\n0YlE/a3RGNMcwCsATgIQAPC8iDwd7XEZqmKnSRPgl190DlDduvE//6xZwEUXxf+8sda5s65XtWUL\n0LSp09VE59Ah4NlnNSQWFwNdu2qQWrtWl+T4wx+Au+4CGjZ0ulIiIvewMfx3GMCdItIJwJkAbjPG\nRD29nKEqdozRa+tEt0oE+PLLxOxUJcq8qm+/BXr0AD78EHj1VZ3b+N57wOTJwJIluhxGUZFuL/T+\n+05XS0TkHlGHKhHZJiJLgr8vBLASQLNoj8tQFVvt2+s1jre1a7Xj0apV/M8dD16fV/Xpp3pX5p13\n6u/POENDeGktWwLPPAO88w4wahTw8MPO1EpE5DZWZ0YYY1oB6AZgfrTHYqiKLafmVc2apd2c49+o\nE0W/fsB//uN0FZGZMQP47W+Bd98Fzjqr4uf37QvMnatbDR04AGRmxr5GIiI3sxaqjDE1AEwGcHuw\nY/UrWVlZR36fnp6O9PT0Ex7rwAFd74cbKcdOu3Y6nBNviTr0V6JLF/23u3On3r3qFUuWAMOHa/cp\nlEBVolkz7cydeaZ2sG68MWYlEhHFVG5uLnKjHGowIhJ1IcaYFABTAXwsIk+V8RwJ9VwrVwJDhzoz\nPJUsFi8Gfvc7YPny+J73lFOAqVOBjh3je954GjgQGD0aGDLE6UpCs2sX0KsX8NBDwDXXRHaMlSt1\n/8MpU4A+fezWR0TkBGMMRCSscRVb61S9BGBFWYEqXBz6i7127XR+U3Fx/M65aROwd2/ibz3Uty8w\nb57TVYRGBBgxArj44sgDFaBf02efBa69Vr/GRETJKOpQZYw5C8B1AAYYYxYbYxYZYwZHc0yGqtjz\n+4EGDYANG+J3zrlzdWgpUedTlTjzTH2tXvDGGzq37pFHoj/WlVcCAwboJHciomQU9ZwqEZkDoJKF\nWo5Yswbo2dPmEelESiarx2vu2ldfaeBIdGecASxcCBw+7O5FMn/6CbjjDl06oWpVO8d87DFgULsC\nFNyUg7QWAV1nIiMDSEuzcwIiIhdz5TY1333HTlU8xPsOwK++0sCR6OrVA5o31+143Oz++3XIr1cv\ne8estasAr54+DpfPHYND92cDY8YA48bpiqFERAnOlaGKw3/xEc9QdeCALhpp8w3czUqWG3CrJUt0\n6YRSN+TakZODtq9lo14LP55/HjrOnJ0N5ORYPhERkfu4LlTt2wfs2eP9bT68IJ4LgC5dCrRtC9So\nEZ/zOc3t86ruvVfXlbK+TVEgAFPDj3/9C/jb3/T/M/x+IBCwfCIiIvdxXaj67jt98/W5rrLE07at\n3gEYD8ky9FfCzZ2qr77STZJHjIjBwX0+oKgI3bvroqCPPw7d04b/oYkoCbjuO92aNfpmT7HXsiWw\nbZsOzcVasoWq9u2B3bv1+rpNVhZw331AlSoxOHhGhrbAiorwwAPAS+OKcPDeTH2ciCjBuTJUcT5V\nfKSkaLBavz7250q2UOXz6RCg29armjdPF+qM2crnaWnAyJHA2LFoNzETYxuPRU7Nkbz7j4iSAkNV\nkmvTJvZDgNu366rdyfZ1dWOo+vvfdT5VTLpUJdLStFuVnY02r2YiOyctLt1QIiKnMVQluXiEqvnz\ndeuSZJtW07u3rlflFqtXA19/Hd+RuO7dgdNOAyZMiN85iYic4qq3ORGGqniLR6hKtqG/Er16AYsW\nuefGt3//G/jDH4Bq1eJ73r/8RSesW9hmlIjI1VwVqn78EahUCahf3+lKkgdDVezUr69bAcVzgdWy\n7N0L/O9/wK23xv/c/frp/L0oN38nInI9V4WqkuUUKH5iHaqKi3XI6fTTY3cON+vVS1+/03JygPPO\n05Xe480YDXPjx8f/3ERE8eSqULV2rb7JU/y0agVs3AgcOhSb469eDZx0km7dkozcMK9KRIf+Ro50\nrobf/haYPh3YutW5GoiIYs1Voer774HWrZ2uIrlUqQI0awbk58fm+N98k9ybY7shVH35pX6dzzrL\nuRpq1QKGDQNeeMG5GoiIYo2himI6BPj118kdqnr0AJYti10nMBQ5OXrHnzHO1QDoEODzz7tn4j4R\nkW0MVRTTUPXNN8mzifKJ1KypQ6x5ec6cv7BQN06+7jpnzl9at246eX/mTKcrISKKDVeFKs6pckas\nQlVxsW6k3KOH/WN7iZNDgG+/rcN+TZo4c/7jXX+93oVIRJSIXBOq9uzRPegaNXK6kuQTq1C1ahXQ\nuDFQu7b9Y3uJk3cA5uTEcEuaCAwfrp2z/fudroSIyD7XhKrvvwdOOcX5eR/JKFahKtmH/ko41akq\nKACWLwcuuij+5y5L06b6b2LKFKcrISKyz1WhivOpnHHKKfoGfPiw3eMm+yT1EqedpktLxLs7M2kS\ncPnlQNWq8T1vRTgESESJiqGKUK2aDrtu3Gj3uMm+nEKJatWADh2AJUvie9633gKuuiq+5wzFZZcB\ns2YBP/3kdCVERHa5JlStXctQ5STbQ4CHD3OSemk9egCLF8fvfPn5wPr1QP/+8TtnqGrWBAYNAj74\nwOlKiIjsck2o+v573vnnJNuhatUqnT+T7JPUS3TvHt9QNWmSdoRSUuJ3znBcfjnwzjtOV0FEZJer\nQhU7Vc6xHao49Hesbt3iH6qGDYvf+cJ14YU6BLh3r9OVEBHZ44pQdeAA8MMPQIsWTleSvBiqYqtr\nV2DFivisrL5+vQ7/pafH/lyRqlUL6NcP+PBDpyshIrLHFaFq/XoNVG4dqkgGDFWxVaMG0LKlDovG\n2rvvAkOHuv//0+WX6+KkRESJwhWhipPUnde6NbBunZ192QIB3e+uW7foj5VIunWLzx2AU6YAl1wS\n+/NE65JLgOnTgZ9/droSIiI7XBGqOEndeX4/UK8esGlT9Mdat06PVbdu9MdKJPGYrL5rl3YJBw6M\n7XlsaNBAFwKdNs3pSiq2dat21caP1zW28vIAEaerIiK3cU2oYqfKeaecooEoWkuW6IKXdKx4TFb/\n5BPgnHOA1NTYnseWiy9297yqhQuBCy4AOnUCXn5Z/21/9JHW3aEDMGGCne4uESUGhio6wlaoWrqU\nQ38n0r3mJhWTAAAgAElEQVS7vinHssMxZYq+4XvFkCEaUtzW9Tl0CPjLX3SI8rLLgM2bgalTgf/+\nF3jtNf1/Mn68fgwcCGzZ4nTFROQGDFV0xMkn600D0Vq6lJ2qE2nUSDtIBQWxOf6hQ9qpctNefxVp\n106vydKlTldyVFGRBtOlS3XvxD/+Eahe/djnGKMLq86ZAwwYAPTpo0OCRJTcHA9VxcV6+/cppzhd\nCXH4L/ZiOVl9zhz9GjZtGpvjx8qQIe4ZAjx4UANV48bamWrQoPznV6oEPPAA8MgjwHnn6bIZRJS8\nHA9VmzbpN67jfxKk+LPRqdq5E9i9myG5LLGcrD51qre6VCUuvFCHAJ0mAowYobsAvPhieEtSXHst\n8M9/aiDbuTN2NRKRuzkeqjj05x42OlXLlgFdugA+x/9luVMsJ6tPm6aTqr2mXz8dOtuxw9k6xo7V\ndcQmTtQOVLhuuEHnX11zje59SUTJx/G3PoYq92jSBNizR+eURGrJEk5SL0/JZHXbtm7Vrm+vXvaP\nHWvVqunq759+6lwNixcDjz6q2/tEc+fkww/r3YCPPGKvNiLyDsdDFRf+dA+fD0hL0zlukeIk9fKd\nfLIOj9oeIpo+Xe9Ci6TD4gYXXKCT7J2wfz9w3XXAE0/ov/9opKTo0gtPPsmJ60TJyPFQtW4dQ5Wb\nRDsEyFBVPp8P6NxZ7yqzado0nSjtVeeeC3z2mTNLK/zrX8Cpp+q8KBtatAD+7/+AG2/UG3GIKHk4\nHqrWr9ef3skdopmsfuiQzknp0sVuTYmmSxe7oSoQ0E7VoEH2jhlvrVsDlSvHZ2/E0vLzgaef1s6S\nMfaOO2KEDiO+9JK9YxKR+zFU0TGi6VStWqWbBntlNW+n2A5Vy5YBdeoArVrZO2a8GaPdqhkz4nve\nO+/Uj5Yt7R7XGODxx4EHHwT27bN7bCJyLyuhyhjzojFmuzFmWTift3cvcOAA0LChjSrIhlNOibxT\nxUnqobEdqrw+9Fci3qFq3jzg66+B0aNjc/yePYHzz9elFogoOdjqVL0M4PxwP2n9ev3p2mbbnaJz\n8smRd6qWLQO6drVbTyLq0sXuhrxeH/orMWAAMHNm/JYjePBB4P779e7DWPnHP4D//Af44YfYnYOI\n3MNKqBKR2QB2hft5HPpzn5I5VWG94RcUANnZ6PZuJq7Iy47dPiwJol49oGZNO5fpwAHgq690E2Wv\na9RIf8hauDD255o1S5dzufHG2J6neXNdt+rxx2N7HiJyB0fnVOXne3seSCKqXRuoWhX48ccQP6Gg\nABg3DhgzBvccyEbVv47RPzNYlatrVztDgAsX6p1rtWtHfyw3KLkLMNays3V7mcqVY3+ue+4Bnn/e\n+cVNiSj2HA1V7FS5U1iT1XNygOxs7D7kx65dQMsOfn3HysmJYYXeF/W8qmB3MPBAJh6qkjjdwYED\nYz+vavFiYPVq4PrrY3ueEi1bApdfDjz1VHzOR0TOCWN3q+hlZWUd+X16ejrWr09Heno8K6BQlExW\nP+OMEJ4cCAB+P76dA3TqFNyexu/Xx6lMXbpEsd9dSXcwOxvZs/y4684iYFwmMHJk9KtXOuyss3Ty\n+IED2jGNhcceA26/PT5dqhJjxgBnnw3cd19s53ARUeRyc3ORm5sb1TFshioT/ChT6VAF6HsAO1Uu\nU1CAEZtzUO35ALDGB2RklP9G7fMBRUXIy/Ojc+fgY0VF3PyvAl26RLGVSbA7eCDFjwULgDPf8QOD\ns3XzusxMm2XGXa1aQIcOwIIFGkJs27gR+Phj4Jln7B+7PO3b692Ar70G3HRTfM9NRKFJT09HeqlO\nT3Z2dtjHsLWkwmsA5gJoZ4zZYIypcPqnCIf/XCfYAdl49RhMOCVbf7yuaH5URgaQmYk1i4s0VBUV\n6Rt7RkacivamDh10ovSBAxF8crA7uGBBqflUCdQd7NdPJ5LHwrPPAr/7nTNz0G6/XYcAnVg1noji\nw9bdf9eKSFMRqSoiLUXk5Yo+56efgCpVEmeCbUIIdkBanOrXOVX+EOZHpaUBI0ei86djccmiTO2W\nJMAwVKxVrao/UES0gniwO5ibi6PD5wnUHTznHF1awaqCAhx+MBsNxmXir4edmYN23nnAwYOxC4xE\n5DzHvgvzzj8XCnZAjpmoHkIHRFqm4a59mUj9V7Z2qRioQhLxZPVgd/Crz4rQvz8Srjt49tm6TMSh\nQ5YOGOzAvt9mDKb2zkaDh525Q9XnA265Re8EJKLE5Fio4tCfCwU7IC1bAlu3Bt/UQuiA/PCDDmk0\nbhyfMhNFxKEqLQ0Hbx6JvnPHYsDMxOsO1q2r3xsWLbJ0wGAH9tkJftx8M0LrwMbI9dcDU6cCu8Je\n1Y+IvIChio4KdkAqHyxCkybAxlWhdUDy8oDOnbkyfriiWVZh8c40TOqYiWqPJGZ38JxzLA6TBQJY\nu9WPvDzgssuCjzk0B61+fWDwYGDixLif+hgHD2o38KWXgBdeAHJz9ecnIooOQxUdFZwfhbFjkW0y\nISF2QEpCFYWnc2dgxYrIPnfOHF1+IFH162dxXpXPhzdfKsK115ZapsHBOWgjRmiQcWLCemGhbs/T\nogVw6616jefM0aUeWrTQ/+7cUococgxVdKy0NCAzE7MGZOOzs0LrgDBURaZVK30DKywM/3PnzgX6\n9rVekmv06wfMnm2nmSQ3ZOCkZzPx28uDrRiH56ANGADs3q17ZcbT7Nn6/3TdOu0CLl4MTJgAvPyy\n/ntatgxISdEO6uTJ8a2NKFHEdfFPADo5NCcHF38dQO+PfECnjIQbukgE4ayqnpent6lTeCpV0vWL\nVq4EevcO/fNEtLswdmzsanNao0ZAgwZ6bTp1iu5Y87el4bX6I/H7z8YCMwLaoXJwDprPp/sBvv46\ncNppMTxR8HstAgEsy/NhVG4GnnklDRdeeOKnN28OPPEEcO21wBVXAGvWaAeLiEIX305V8C6cwJ1j\nMKYoG/5M7hPnViUbK1dEBPj22+jf+JJVx47hDwHm5+v8tUT/WeTMM4F586I/zsSJwIAb02CyMnWC\nugvmoA0fDrzxRgyHAEvtyflGh2xcPncMpl08Dhd2rvh7be/ewPz5ulDpP/8Zo/qIElR8Q1XwLpyt\ne/2oUwdIbch94twqLS20rLthA1CzJlCvXuxrSkSRhKqSob9EvzHARqg6dAh4803tvrhJ165Aaqqd\n0HhCwe+1Xyzw4/bbgXen+dHgmdC/1zZpAkyfDrz4og4PElFo4huqgusgHTOfKoFWgk4krVqFFqrY\npYpOJKEq0Sepl7ARqqZNA9q21eFsNzHm6BBgTAQC2LjTj2uv1Y5Tly4I+3ttkybAlCnA3XcD33wT\nozqJEkx8Q1VwHaRjQlUCrQSdSJo0AXbuBH75pfznrVihwYAiE02nKtF16aJ79UWzptMbb7ivS1Vi\n+HDgrbeAw4ftHztgfBgxvAijRgEDBwYfjOB7bYcOurXPsGGR3VBBlGzim2aC6yBtWl2koSrBVoJO\nJD6f3mK9YUP5z1u5Ur/xUmRatwa2bAF+/jm05+/dC6xdC3TvHtu63CAlBejVS+f3ROLgQeDDD0ut\nTeUybdvq/7EvvrB/7P9VysAN6zJx923R3/F41VXAb34D3HOP1RKJElJ87/4LroPU4dKxaNk8AIx1\n9i4cKl9amk6Kbteu7OesXAnccEPcSko4KSlAmzbA6tWhBaX584EePXTfzGRQMgQ4eHD4n5ubq3dX\nNm1qvSxrrr4aePttYNAge8fMzwfufCoN894YiUpPjNUhvyjveHzySe0cDh+eHEPPRJGK/5IKaWl4\nqk4m7r8DwMAKn00OqmhelQg7VTaUDAGGEqrmzk2uN7Uzz9Sb2CLx7rvu7VKVGDpUN8V+9ll7syDu\nvBO44w6g7blpwLmZVo5Zty7wyCPA7bcDCxZwxgZRWRz5r7F+PTdT9oKSTlVZtm/XtZYaNoxbSQkp\nnHlVc+Ykx3yqEmecoW/i4d7LEggA773n/lDVrh1Quzbw9dd2jjd7tk4qHz3azvFKu/ZaoHJl4JVX\n7B+bKFHEPVQdOqSb9bZsGe8zU7gq6lSxS2VHqKEqENCAceaZsa/JLRo21I9wJ/N/9ZUuHtq2bWzq\nsmnoUOD996M/jggwZgzwj38A1atHf7zjGQM89RTw179yn0CissQ9VG3aBJx0kv7EQ+5WUaeKocqO\nTp1CCw1r1uh6YA0axL4mN4lkaQUvDP2VsBWq3n5bJ+dfd130xyrL6afr8PN//hO7c/xKQcHRRVuz\ns7lYNLla3ENVQQGH/ryCnar4aNNG77I8cKD85y1cGN52Noki3FAl4q1Q1acPsGMH8P33kR9DRDtU\nf/tb7Oc73X+/bpEU6h2rUSm1Mjyys/VX7sJBLuZIqOLQnzc0a6Yb/h48eOK/Z6iyo0oVXbdtzZry\nn8dQFZo1a3R9tW7dYleTTT4fcPHF0XWrPv5Yg1VZ+/rZ1LWrfk3++9/YnkcEyM/OQSayMexGP666\nCrjn737MGpgNeTknticnipAjoYorKHhDSoouArpp04n/nqHKnlDmVSVrqOrcWRcB3bs3tOd//DFw\nwQXe2sYnmiFAEeChh4B7743fa37wQeDRRyvurkZq1iwNxR9N0V04rrgCuPJKoFo14La7/fjPswF8\n9llszk0UDYYqKldZewDu2aMfLVrEv6ZE1KGDbvlTloMHgWXLgJ4941eTW6SkaHdk8eLQnl8Sqryk\nf39g0aLQg2Nps2drR/mqq+zXVZZu3XQu4Ftv2T1uIKDh8LrrdJjx1tt8yL67CFdfrWt6ZWUBy+YV\nof9AH266SbfQKS62WwNRNOIeqjZsYKjykrImq69apQsrcr0aOzp00AVAy5KXp0OENWrEryY36dUr\ntGUHfv5Z1/I699zY12RTaqoOqX3+efifO26crh9VqZL9uspz++16N6CIneMdOqTb4cybpwH6qqsA\nc2OGTlAvOroyvMnKxKkPZ2DRIv03ceWV+rlEbsBOFZWrrMnqHPqzq317DaplSdahvxKhhqovvtBu\nXq1asa/JtsGDgU8+Ce9zNm8Gpk8Hfve72NRUniFDtFs9d270xyou1h109u/XTbCP3OEa3IUDY8dq\nuBo79sjK8PXrA59+qp97/fXsWJE7xDVUBQI6N4JDRt5RVqeKocqudu2A774re5FLhqrQQpUXh/5K\nnH++hoRwOj/PPaeLcjoRIkt2vnnqqeiPlZmp7w2TJ59gC6a0tKPLKWRmHvNTeeXKOgT5ww+6fhaR\n0+Iaqn74QYcv/P54npWiwU5VfNSsqWtQbdx44r9P9lDVvj2wbRuwa1fZzxHxdqjq2FGHsSq6C7TE\nwYMaqm67LbZ1lScjQztLP/wQ+TE++EBXaZ88ObJFS6tV02D1+uvAO+9EXgeRDXENVRz6856yOlUr\nVjBU2VbWEGBRkXaxunaNf01uUamS7o24aFHZz/nuO70brUuX+NVlkzFHu1WhmDpVO5wdO8a2rvLU\nqgVceinw6quRff727cAf/gC8+SbQqFHkdTRsqKHsllt0xw4ip8Q1VHGSuve0aAFs2QIcPnz0sQMH\ndJmFNm2cqysRtW9/4snqixfrnVZVq8a/JjepaAhw+nRg0CBvLaVwvMGDQw9VL78M3HhjbOsJxU03\nAS++GP6EdRHg1luB3//eztZLvXtrqLr5ZnuT54nCxU4VlatqVZ00umXL0cfWrtVhQW41ZNepp544\nVCX70F+JikLV558DAwfGr55YOPdc4MsvdfHS8mzbps+78sr41FWes8/WH7q++iq8z3v/fe3MZmba\nq+X++4H167VrZdOhQ8DSpcBnn+mG1fv32z0+JQ6GKqrQ8fOq1qzRYQeyq6zhv4ULdc+1ZNerl76h\nnUhxMZCbCwwYENeSrKtbV4fzKlpBfuJE3YbHDUtsGHO0WxWqAwd0x5mnnrLbga1SBXj6aV2/qqJg\nGor167Xz1aABMHy4bgV00026f+3w4bp2HFFpcQ9V3KLGe46fV7V6tQYAsqus4b+vv9ZAkezatNE9\n8nbs+PXfLV2qc3KaNo1/Xbb1769LQ5RFBMjJ0UnibnH99TpJPNQV1seN0zmZgwbZr6V/f+C004An\nn4z8GCLAY49ph7hxY/1/uWKFfl2WLtXvh716AeedB/zlL1wni45ip4oqxE5VfLRsCezcCezbd/Sx\nfft0LaJTT3WuLrfw+YAePU7crfrsM+93qUpUFKqWL9eV188+O341VaR5c71B4OOPK37ujh3AI4/o\nklOx8uijevxI7ko8cEBD4ltvAQsW6EoOjRsf+5x69YDRo7VTtXy57rm4Z4+d2snbGKqoQuxUxYfP\nB7Rte+wt9UuX6iT1lBTn6nKTsuZVJcJ8qhJ9++rNCT//fOK/f+stXXncbbsZDB+uyxpU5IkndOgy\nlt9D2rbVFdkffzy8zzt8WD9v/34dTj7llPKf36iRLgnRpg1wQ3oBDvw1++iaWidai4YSXly/VR8+\nrAmfvKVVq+DEz4ICICcHly4OoPsUH5CWwZRsWckQYM+eAAoKUPn/cvBgcQDI9ul4T5Jf71699Pb7\nIwoKcPjFHJzzeQCDT/MBPTM8f41q1NDhqzlzjhseKyiAvJyDxuMDuOwKH1CQ4arXeuWVOhRWWHiC\nuV7B7x2//BxAjXE+/PHTDACxrf2ee7SzedddQP36FT9fBPjjH/V96u23Q78RJyUFeObuAky7eByG\nzc/G25/4kXKgSMNVcPV3Sh5x/VknLc3btzsnq7Q04OB3BcC4cdh50xj8o3I2amSO0YkR/GnMqiN3\nABbo9c5pMAab/pCts3p5vdGzZ6lOVfAaLTh7DCZ1zkbqA4lzjX41BBh8rXmDx+CxWtlo+rj7XmuD\nBsBvfqN39R0jWDvGjMG//NlYd9kYtHw/9rWnpWlHrNwV3wsKjqzUvvCibGybX4BJk8K/s9lMyMHA\n2dk4kOLXld39fj1uTk7YdR88qEOK8+bpndZl7bJALiUicfkAIIMHC3nQzz+L/K1SlhTvLZR580R6\n9Qr+RWGhSFaWo7Ulmv/9T+Tqq0Wva2GhdO8uMm9e8C95vaW4WKRmTZEdO+TINcrOFhkzJviEBLlG\nM2aInHFGqQeCr/W++0Tuvjv4mAtf64QJIkOHHvdgsPbCQpEGDURWr5a41b52rUi9eiJ7957gL/Pz\nRUaPFikslK+/FklrUCg7fz9aHw/Xgw+KiMiPP4q0aCEyZcqxj59Qfr5egwcfFMnKktXT8uX66/Xf\nd/v2In36iKSliTRqJPLnP4ts3Rp+WRQdjUjhZZ24dqpu+4FjzV5UvTpQvWoA2/b5j51P5ffzxyjL\njiyrEAjgYGU/Vq0qtZI6rzd8Ph0aW7oUei38fuTmAunpwSckyDXq21e7FUduWggEIKn+I/OpALjy\ntV50kXbYjpkPFvw6vfqqdrLatUPcam/dWm9gmDDhBH+ZkwNkZ+NgZT8yMoD/e8qPuk9F1l2CzwcU\nFaFBA+B//9NlGHZtKip74lup7t2h+7Px4N4xmDF0HH7TogDff6/fA776SueyzpmjQ5OdO0dWWulu\nHN9/Yy+uoWr5lRzG8KqadXzYvKbo2Dv/isr5pkERad8+uLGy8WH1oiK0agWkpgb/ktcbANCtG7Bk\nCQCfDwd3FWHhQuCss4J/mSDXqHp1HeqcMyf4gM+HZfOKEAjoPCEArnyt9erpvLdp00o96PNBCovw\n9NPAqFHBx+JY+6hR+pbzqwwXDHuPPqp33g4fjsjDXkaGhpaiIvTrBwy7sAgLLsgse92LYKDbW+zH\noEHAwhV+XLUiGzdXzUHDhsc+tU0bndyfm6t3Nd5yi67LFpJS4Q3ZfP+Nh7jPqYpmrJmcs7hrBmo9\nlon8b4u0U1UUnIjppsVyEkDNmkCdOsDW8zJgsjLRp3OR/gWv9xHduundccjIwI7/l4lOrYpQpw4S\n7hr1769vpACAjAwUjs7ENRcX6bxUF7/Wyy4D3nuv1AMZGdhwYyZqmCLtKMa59t/8Rn8wOSboAYDP\nh7VLi/DEE8D48Th6XSMJe2lpOil97FggMxP/bDAWmT+NxLwtZUxSDwSw+5Af556ri71++CHQsFX5\nga5zZ2D+fJ1zmZERYvYLhrcDKX58+SXwv3f9mNwlG1sfznFbkzNhxPXuvyM3QbiwbU3lq35qGmam\njsRZL4xFv0YBYJWPd7bESPv2wLeFaZjXdCRGbBoLZAb0Gz2vNwANVU89BSAtDVNajcQD3yXmNTr7\n7FJbuKSl4aG9IzF+r/tf69ChQFaW3kWXkgIgLQ1/2zUS41uPhcmKf+3GaLfq6ad1b8UjMjKwYlAm\n7hmZjZYt/UfD3siRkZ0oLe3IF6wagD+dCtx5JzB37q9v0CoWH357eRF69fLjmWdCD3Q1amgAGzxY\n727817/KL6lwbwD3/9WPnBzteLVrBxw65Me2zwN49CPgz3/W/RfDXtU+eDcnAsGvJ+9MPircSViR\nfgCQDRuCs79cOMGSyvfEEyK33SZSvbrIvn1OV5PYbrlF5N//FunXT2TaNKercZ/9+0WqVdNfL7lE\n5I03nK4oNvbtE0lN1de5YYNI/foihw45XVVoevYU+eIL/f2mTSJ164oUFTlXz/79ev1Kz0GfO1ek\nT+N8OfDXo5PFI5qkXobiYpHu3UXefPPXf/fgDfky+eTRcmh3oT5QWKiT5kM8/08/ibRpozcGlOXN\nN0Ue8WfJ3bcVyqZNpf4i+P67YIHIkCEinTuLLF0a+usqPcE/ktq9BBFMVLfSqTLGDAbwJHQ48UUR\neeREz2v2QjZwzTDdJCrSnwbIEe2rFaDq2zl4KCWAGo/xJ5NY6lG/AM1eycEFSwLoO8MHtMvgtS6l\nWjXg7JYF2HF7Dn4zPYAL2viAMzIS7hrVqAEMaF2A7bfl4MdtATzXzIeUzRmeeJ2XXgrMfKUA6TNz\nsH5WADkn+5D6Y4ZjtVerBlxzjU5Yf/AGXe9r1UsBPHu6D1X+EJu6fD7d6mbECB0SrbxFuztrVgVQ\n52MfznvzMqQ8OfZotyeM7l39+sC77+oQ8Zln6mKnJd0jKQ5g5pc+PL02A8++noGuMzOBOtkAju3G\n9U4Dpk4FXn1VN/J+5ZVSnbzjOlG7hmYg54s0fPwxMHhBDv7jz0aThX4MHgxkZPjRJDv7yNBn0gs3\nhR3/AQ1Sa6EruVUGsATAqSd4nsg994iceabI7NmxDphkU36+bL1utKSiUAYMkIT+ycRx+fmy9tLR\n0q1toTRvLrzWJ5KfLx92GC33jiqUk0+WxL1G+fnyec/RMja7UC64QGRSjnde56pP8+W/tUZLYF+h\ntG0r8tVnztf+zTciZzXPl8Cdo+WzDwqlfXuRw3tiX1d6usjkx7S7s31doTRuLDL7Uzvnffppkd69\nRQ5/f7R7lJkp0qdzoRTeEjz+cUs3nOics2eLNGwY7IyX6kQdPizy70cK5d/VRssdl+XL+++L7Bz5\noKxfL/LJJyI336xdyMxMkcN/LWf5CI9CBJ0qG6HqDAAfl/rzPQD+coLnibRqpQvxcOjPW7KyZPff\nn5SdqCNFlWuJ1Kkj8uST/DrGQlaWbLtXr/W+SrzWJ5SVJdMvflJ2mTpSlJLA1ygrSxZlPCl7U+rI\nbtSS4treeZ2BYcNki6+p/FK9tuwxtSUwbJjIihWO1h4IiOTWvFAOpNaWfZVqyS+pdUTGj4/5dJSF\nT86WbZUaS6BWLSmsXFs2NuquKcbCeQMBkfs7Tpb9VWuJ1Kolv1SvLVP8w+SHhflhH//LL0VuqDlZ\nDqXqsYpr1ZbcJsPkil75suqbUscKrjtWYuNGkaHnFspzzbP0vBUEOC9xKlRdAeC5Un++HsDTJ3ie\nyLp1IiedJHLVVTG+FGTVkCESSE2VLpXy5IknRCQvTyd7DBnidGWJJ3itOyJPHnhAeK1PZMgQOVQ1\nVdojT55/XhL3Gg0ZIsXV9XUOGCDeeZ2TJ4tUry5fdhghqSiUZ+5epytYXnyxyKhRztU1frwc9qXI\nec3zpGVLkYNLgtdz/PjyF+mMxuzZEmjXTn5KOUnuuWm7nNZkuxzu3EVX95w9O/rzTp4sB+s0kCWV\nusuUNwqle711cqhe8Frn54d3/MmTZX+NBrKscnfJ/7ZQhnRcJ7urNZLDFx53rBPMqQrcOVr+e8Ns\neaHOaNmwMnHmWrk/VInoN4Y6dWJ4Gci6atVEFiyQtm1FPvww+NiCBfo42RW81m3aiLzzTvAxXutj\nVasmP05bIIDIqlXBxxLxGpX6f/fEE8HHvPA6W7US6dNHpry4XQCRzZtFf6Bu3lzHwpxSp44U3vQn\nSUWhPPZY8LG8PJHatWPXqUpPF+nbV7677C5JRaG89ZaIbN8u0revyNlnR3/eVq1Ebr9dZpxxn6Si\nUN59V/Rat2wpct994R0/eKxPe+mxRo0SCXxfxrFONKSYlSVPPlQo7dqJ/PBD8HkevynNyeG/T0r9\nuczhvxRABgFyByDgh2c+dgPyKCCpGCBAdUkN/nm3C2pLtI9jr7Wf17rca3SxAEjYa3T0dfYToK5n\nXud6QJ4A5FH4JBVDjnyNdgcfd6quXYC0BORRpEkqqhypa3/w8VicMzf4oedtI6nBx2cBssHCedcD\nkgVIS1STR3HKkePnAzInzOMfPVaKPIp2YR8r68jvRwhQ9wSPe/Mj3ExkgoEnYsaYSgBWAxgIYCuA\nBQCGi8jK454nkpUFDBwI/Pa3wPr1UZ2X4qhuXeCtt3TBlZI7Vfr21f0ydu1yurrEwmtdsWS5Rl59\nnSefrItV3Xyz1h8IAHv2AJMmAX/8o3N3iNWtC8yerbdVltzZ9tNPuq/Mnj2xOWf//rpD8rvvAvv3\n63n37dPr0qyZ7pocjZJr/dBD+lpycoAdO/RWvhtv1KXY43WskhXb/f6jjxUVefquQGMMRMRU/MxS\nou1UBUPZYGiw+g7APWU85+icqsmTI27HkQPGj9e5B3l5+ue8UnMRyC5e64olyzXy6uucPFl3Th4x\nQrUxF3oAAAWnSURBVId/1pWaU+Xk/Bonrufs2SLt2omcdpoO+23fLtKl1JyqaNm81tEeKwHXr0IE\nnSoroSqkEwE6ZstA5U3jx+tcuFrBu63c/o3dy3itK5Ys18irr3PyZJ1DVbu2fgwb5o43Vyeu5+zZ\nugpoybXo3t3uskI2r3W0xwph+QYviSRURT38FypjjMTrXERERETRiGT4z11bnBMRERF5FEMVERER\nkQUMVUREREQWMFQRERERWcBQRURERGQBQxURERGRBQxVRERERBYwVBERERFZwFBFREREZAFDFRER\nEZEFDFVEREREFjBUEREREVnAUEVERERkAUMVERERkQUMVUREREQWMFQRERERWcBQRURERGQBQxUR\nERGRBQxVRERERBYwVBERERFZwFBFREREZAFDFREREZEFDFVEREREFjBUEREREVnAUEVERERkAUMV\nERERkQUMVUREREQWMFQRERERWcBQRURERGQBQxURERGRBQxVRERERBYwVBERERFZwFBFREREZAFD\nFREREZEFDFVEREREFjBUEREREVnAUEVERERkAUMVERERkQVRhSpjzJXGmDxjTLExpoetooiIiIi8\nJtpO1XIAlwGYaaEWcrHc3FynS6AI8Wvnbfz6eRe/dsknqlAlIqtF5DsAxlI95FL85uBd/Np5G79+\n3sWvXfLhnCoiIiIiC1IqeoIxZjqAk0o/BEAA/FVEpsSqMCIiIiIvMSIS/UGM+QLAaBFZVM5zoj8R\nERERUZyISFjTmyrsVIWh3BOHWxgRERGRl0S7pMKlxpiNAM4AMNUY87GdsoiIiIi8xcrwHxEREVGy\ni/ndf8aYwcaYVcaYNcaYv8T6fGSPMaa5MeZzY8y3xpjlxphRTtdE4THG+Iwxi4wxHzhdC4XHGFPb\nGDPJGLMy+H+wj9M1UeiMMX8OLo69zBgz0RhTxemaqGzGmBeNMduNMctKPVbXGDPNGLPaGPOpMaZ2\nRceJaagyxvgA/BvA+QA6ARhujDk1luckqw4DuFNEOgE4E8Bt/Pp5zu0AVjhdBEXkKQAfiUgHAKcB\nWOlwPRQiY0xTACMB9BCRrtD5y9c4WxVV4GVoVintHgAzRKQ9gM8B3FvRQWLdqTodwHciUiAihwC8\nAWBojM9JlojINhFZEvx9IfSbejNnq6JQGWOaAxgC4AWna6HwGGNqAThbRF4GABE5LCJ7HS6LwlMJ\ngN8YkwIgFcAWh+uhcojIbAC7jnt4KIAJwd9PAHBpRceJdahqBmBjqT9vAt+UPckY0wpANwDzna2E\nwvAEgLug68qRt5wM4CdjzMvB4dvnjDHVnS6KQiMiWwA8BmADgM0AdovIDGerogg0EpHtgDYZADSq\n6BO4ojpVyBhTA8BkALcHO1bkcsaYCwFsD3YaDbiVlNekAOgB4BkR6QHgZ+hQBHmAMaYOtMuRBqAp\ngBrGmGudrYosqPAH1FiHqs0AWpb6c/PgY+QRwdb1ZACvisj7TtdDITsLwCXGmHUAXgfQ3xjzisM1\nUeg2AdgoIl8H/zwZGrLIG84FsE5EdopIMYB3APR1uCYK33ZjzEkAYIxpDOCHij4h1qFqIYA2xpi0\n4J0P1wDgXUje8hKAFSLylNOFUOhE5D4RaSkip0D/330uIr9zui4KTXDIYaMxpl3woYHgDQdesgHA\nGcaYasYYA/368UYD9zu+q/8BgIzg728AUGFjweaK6r8iIsXGmP8HYBo0wL0oIvyH5RHGmLMAXAdg\nuTFmMbT1eZ+IfOJsZURJYRSAicaYygDWAbjR4XooRCKywBgzGcBiAIeCvz7nbFVUHmPMawDSAdQ3\nxmwAkAngYQCTjDE3ASgAMKzC43DxTyIiIqLocaI6ERERkQUMVUREREQWMFQRERERWcBQRURERGQB\nQxURERGRBQxVRERERBYwVBERERFZwFBFREREZMH/B9Bd3p97HvJ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107d619b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(sol_smooth[0], sol_smooth[1], '-')\n",
    "plt.plot(sol_points[0], sol_points[1], 'ro', fillstyle='none')\n",
    "plt.hlines([0], 0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 3. Van der Pol oscillator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system for Van der Pol oscillator is given as:\n",
    "$$\n",
    "y_1' = y_2, \\\\\n",
    "y_2' = \\mu (1 - y_1^2) y_2 - y_1 \\\\\n",
    "y_1(0) = 2, \\quad y_2(0) = 0\n",
    "$$\n",
    "It becomes stiff for high values of $\\mu$, meaning that regions of rapid transition are followed by regions where the solution varies slowly. Explicit methods either diverge or make prohibitevely many steps for stiff problems, thus implicit methods should be used. Our function `solve_ivp` implements a one-step fully implicit Runge-Kutta method of Radau II A family."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve this system for $\\mu = 10^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 1e3\n",
    "\n",
    "def fun(x, y):\n",
    "    return [y[1], mu * (1 - y[0]**2) * y[1] - y[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, [0, 3000], [2, 0], method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x107d882e8>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVNWd9/HP6VXQQkFosEUQRMSuAlxAWqSw1VEJmmg0\nk1ET80iMmSySeT2ZmRiTyYiPSSaLZmIWjYNbojGaxMREQVGiTbfQIk8AoarZomERhGaVi0BDdZ/5\n41ZvyNbUvbXc+32/XqVVl+pb5/Q99etTv7OUsdYiIiLBV5TrAoiISHYo4IuIhIQCvohISCjgi4iE\nhAK+iEhIKOCLiIRESaYnMMaUA3VAWfp8v7fW3p3peUVExFvGi3n4xpie1trdxphiYC7wFWvtmxmf\nWEREPONJSsdauzt9txy3l6/VXCIiecaTgG+MKTLGLAI2Aq9Yaxd4cV4REfGOVz38VmvtucBAYJwx\npsqL84qIiHcyHrTtzFq70xjzGjAJaOz8b8YYpXlERI6BtdZ4cZ6Me/jGmL7GmBPT93sAlwPLD/Zc\na21gb3fddVfOy6D6qW6qX/BuXvKih38K8EtjTBHuH5BnrLUzPTiviIh4KOOAb61dCpznQVlERMRH\nWmnrkZqamlwXwVdBrl+Q6waqn3TwZOHVUb2QMTZbryUiEhTGGGy+DNqKiEhhUMAXEQkJBXwRkZBQ\nwBcRCQkFfBGRkFDAFxEJCQV8EZGQUMAXEQkJBXwRkZBQwBcRCQkFfBGRkFDAFxEJCQV8EZGQUMAX\nEQkJBXwRkZBQwBcRCQkFfBGRkFDAFxEJCQV8EZGQUMAXEQkJBXwRkZBQwBcRCQkFfBGRkFDAFxEJ\nCQV8EZGQUMAXEQmJjAO+MWagMeZVY0zSGLPUGPMVLwomIiLeMtbazE5gzABggLV2sTHmBOCvwDXW\n2uUHPM9m+loiImFjjMFaa7w4V8Y9fGvtRmvt4vT9XcAy4NRMzysiIt7yNIdvjDkdOAeYf7B/nzpj\nKht2bqBhXQNOs+PlS4tkTfUvqjF3G0bcP4LZ78xWW5aCkXFKp/1EbjqnFrjHWvung/y7ZRoUmxKM\ngWi/KDNvmsma99cQq4gRKY94Ug4RP1X/opr5m+aDBdIfsmP9Ysz69Cy1ZfGFlykdTwK+MaYEeAF4\n0Vp7/yGeY7m404HBRZx8xlDeZzVn940y99Z6vVEk75m7D/K+aymm574h7O2xmiHHR5n7uXr6n6S2\nLMemtraW2tra9sd333133gX8XwFbrLVfPcxzLNOgpKgEMJxcPJjN+1bTalKQKmXUwjomj6omHofx\n4+GkkzIulojnht8/nFU7VnU5NjgylHedtbSQgpZSyp+qY1Sfai66iPbbCX0cEk0JfQKQbsurQVtj\nzEXAp4BLjTGLjDELjTGTDvbcqWOmsuZf1lA/pY6Ft89h5IAopUWlxE6p4p6pUcrK4L774LTTYPRo\n+PKX4emnYfk7jvL+khfunXRvl8f31NzDvM/VE+vvtuXRlVWs/f9R7r0XKirg8ceh6hyHk78W56KH\nJ1J1b5wlKxw0YU1ywbMc/hFf6CDTMp1mh+TmJNF+0S69nv37YdEiqK+H1+Y5vHRKnJY+SXq3RPnh\n2fVMvizCKadkpdgiXcx+ZzaXP3F5x+ObZ3PZ0MsO2ZYB5q5toObxiaRsiqLWUnr/qY6e26qpqaH9\nNmQIGE/6cBI0eZfDP6oXOsZ5+A3rGpj4+ERSrSmKKSX+tzqWzKymf3+49FL3dvHFcPLJPhRa5ABO\ns8OIH47nvdQyohVnM+/WeUdM0TjNDvHH4jRubqSqXxV1t9SzcW2E2lrabyUl6A+AHFSoAv6Bb5b6\nKfX0LInw1lvw6qvu7fXXYdiwjj8A8ThElCYVn/zwwQ38efkMnrnnKip7VR7VzxzuE4C1sGoVh/wD\nMOYih109Eozsr/x/GIUq4MPh3yzgpoAWLOj4A/DmmzBqVMcfgAsvhB493PNo4Ewy4TQ7VN0bZ8P+\nJCMHRKmf4v3sss5/AF6Z4/Bcnzip3klObo3y03PrueryCL16efqSksdCF/C7a88emDev4w/A0qVw\n3oUOKyfE2WqSVPWL8vpnNQ1Uuq9hXQMTHp1IKylKi0qpm1JH9cBqX1+vLaVZRCljltTR+HI1558P\nkya5t9Gjlf4JsryapZOPevSAyy6D73wHGhpgwwb42OcSbLZJUjbFkvcaueFfkjz7LOzcmevSSiGJ\nVcSoLIlSZEup6ldFtF/U99eL9nNnAI3sX8Xsp6Js2gR33OG2609+Eior4ZZb3BltW7f6WhwpcIHs\n4R9M57GAM06s4lZTz+wZEebOhepquPpq93bGGTkrohSIY8nhZ+JIKc2334ZZs+Cll9w0UFVVR+9/\n7FgoLlY6s5AppXOMDvbG2bULZs+G55+HGTOgd2838H/0o+4CsJKSnBZZ8kw2cviZaG6GuXPd4P/S\nS7B+PdRc6bAgFue9Frft51uZ5fAU8H3S2gp//Su88IJ7+/vf3V7S1Ve7/+/TRz2lsMt2Dj9T69fD\nz//cwPc2TsQWpTAtpfxb3zr+7YZqKipyXTo5Ggr4WbJ+Pcyc6fb+a2shdr7D6kvibCZJtEI9pTDq\n6OE3MnJAVUG0gc7pzMqyKsYscdOZ55/vjgFcdx3065frUsqhKODnwJ498MDzDXytcaK7/09LKV+J\n1PG1m6o5Vbv/h8pPH3KYvSTJkz86eE49Hx2YztyzB158EX73O/f/Y8a4wf/jH1fwzzeapZMDPXrA\n56+Jte//MzRSxbblUUaOhEsugf/5H82QCIvjiiJU7KsumGAPECmPUD2wo8w9erg9+9/8xp3t88Uv\nulOYhw2Dyy932/PmzTkutHhOPfxuOrCntHev20P6zW/cmRITJ8KNN8LHPga2VPn+IJo+3V3cN316\nrkvivd273fb829+6g77jxsE//qPb8y+PqD3nglI6ecpx4Lnn3OD/+psORbfF2XWcm+/XQq/gCHLA\n7+yDDzoF/1cd7Gfj7Dk+SVVFlLlqz1mjlE6eikTg5pvdgd5nahM4xyVpwV3o9Y0fJ9m2LdclFDl6\nxx8Pn/iEG/Cfa0iw53i3PS99r5GvfDvJ2rW5LqF0lwK+TyacGWNkeo/0M0+qYktjlDPOgFtvdbd+\nFikkYwfF2vf8P6tPFcftjHLuuXDNNW4qs7U11yWUo6GUjo8OzPc3NcEjj8CDD7pf8nL77XD99VBW\nluuSSnf89CGHV5Yk+PWPwpXLPrA9f/CBm7584AF3i5IvfAGmTIGyE5Tr95Jy+AUulXLn9v/sZ9DY\nCLfdBp+a4rCtRG+SfJfvK21zwVp3TOOBB+C5mQ7Ft8VxjtNaFa8oh1/gSkrcWQ9/+Yt727jdoeo+\n9yvwxj4Yx2l2cJr1tY75KNGUYEMqSatJ0bi5keTmZK6LlHPGuLN5fvlLePq1BO+Xd2xS+Nta/X7y\niQJ+jlVVwZSvJSjqn8QWpVixrZHrpr7JBb+IM/HxicQfiyvo55Fs75ZZaDqPXVWWVvGfX4wyeTLM\nn486MXlAKZ080Hnp+1l9qjh3y7080fIRKE5RbEr47yt/xNn9zmbcqeP08TgPFOJK22zqnOsvI8Jj\nj8G3f+iw8/o4e05Qqqe7lMMPoM5vEoALp8dp3JLEpoqhpBkMxPrFjuo7VMVfYZmH76U57zRw6RPu\npnNFtpRZN9QxbkhUg7tHQTn8AOq89D1SHqHhtnoe+ugDFJemIH2pGzcv4+nE0/pILAXnvFM7Uj0n\n7a/ipqsGEb1PactsU8DPU5HyCDfEbuDsfme3H2tNlfDPz3+J8Y/oDSKFJVIeoX5KPXVT6lj9n/X8\n1wNrWLc3SapVg9/ZpICfxyLlEebdOo/ZN8/mJ5N+QklZC9akSGxs5Kk5b2oALEf2tjpsKtPvvrs6\nf4r9ZE2MUae4g9+2qYo1SwapPWeBcvgFovPAbv/is3jvPbB9lzG09xDm3DInK1+1J5qH76W2cavV\nbw3i5pcm09I3wdDeQ6hTe+5Cg7Yh1fYG2bVvF5OemEQLLQD07dmXuVPmMrzv8ByXMPgK7RuvCkHD\nugbij8VpsS1gYWjvoTz00Yc0Ky1Ng7Yh1faReNyp4xjSe0j78S27txB9MMrKLStzWLpw0Dx878Uq\nYgw5Kd2eDbyz/R0mPTmJ8Y+MV4rHYwr4BShSHmHOLXPo27Nv+7FUa4oJj01gw84NOSxZ8EXKI3yt\nbz1Xb65TOscjbe15WJ9hFKVDUottIbE5Qe3q2twWLmAU8AtUZa9K5k6ZS0lRSfuxzbs3c9FjF6lX\n5LNC/MarfFfZq5KFn1/Ijyf9uH0aMsDnn/+8OjEe8iTgG2MeMcZsMsYs8eJ8cnSG9x1O8otJ+hzX\nB9LDI6t3rObhhQ/ntmAixyBSHuGWc27h9BNPbz+28YONjJ0+Vp0Yj3jVw38MuNKjc0k3DO87nDsu\nuqNLr+irL3+VhRsW5q5QAadpmf6JlEe4f9L9XY5t2LWBZxufzVGJgsWTgG+tfR3Y7sW5pPs+PerT\n7p1Ok6DunH1nbgoTcE6zww+2xJnRTytE/XLJkEvoYXp0ac+/Wfqb3BUoQJTDD4DKXpWcdsJpXY6t\n2roqR6UJNm2P7L9IeYRBJw3q8ql1/fvrc1egACk58lO8M23atPb7NTU11NTUZPPlA+39Pe93eYNs\n3b01d4UJsLZpmRv2N2papo82Ohu7PH7XeTdHJcm+2tpaamtrfTm3ZwuvjDGDgeettaMO8e9aeOWj\n0T8fzZItHWPmo/uNZvGXFuewRMGl7ZH9p/bcIV8XXhm69DElm2L9Y10fV8QO8UzJlKZl+k/t2R9e\nTct8CpgHDDfGrDXGTPHivHL0JgyecNjHIoVE7dkfnuTwrbU3eXEeOXbOPuewj0UKidqzPzRLJyBG\nnDyi6+O+Iw7xTJH8p/bsDwX8gBhTOYYiysBCWVEZ559yfq6LJHLMxlSOoVjt2XMK+AHRuKWRVvaB\ngX2t+1i2ZVmuiyRyzBq3NNKi9uw5BXwRkZBQwA+Iqr5VlFAOFsqLyzm779lH/iGRPFXVt4pSo/bs\nNQX8gFjz/hpaaQEDrbaVtTvX5rpIIsdszftraLVqz15TwA+IWEWMvozA2GLOOvksLfn3kXbL9F+s\nIkZl+QhoVXv2kgJ+gJhO/xV/aLfM7FF79p4CfkAkmhJsZjnWpFixdYV2cfSJdsvMjkRTgvXNy6FI\n7dlLCvgBoZROduhLzLNDKR1/KOAHiD4C+09fYp49as/eU8APCKV0ske7ZfpPKR1/KOAHRKwiRgVK\nNUgwxCpiDDwuimlVe/aSAn5ARMojfKZoJle2/JyZN81U71MKWqQ8wrdOn8nod9WevaSAHxBOs8Ov\nWiczq/hLTH5qsqYLSkFzmh3uWT2ZtwaqPXtJAT8gEk0JmtB0QQmGRFOCd/cmsUVqz15SwA8I5fCz\nRytt/accvj8U8AMiUh7htuJ6bk5puqCftNI2OyLlEb57Rj3xv6k9e0kBP2AsNtdFCDSttM02tWcv\nKeAHhNPsML0lzpMl6nn6SStts8NpdvjG23Hqh6k9e0kBPyA0aJsdWmmbHRq09YcCfkBo0DZ7tNLW\nfxq09YcCfkBo0FaCRIO2/ijJdQHEO+Umwqm2mkh5rksikrkexRH67FZ79pJ6+CIiIaGALyISEgr4\nIiIhoYAfIM3W4V2jJf8SDHtaHLb1VHv2kicB3xgzyRiz3Biz0hhzhxfnlO7Rwqvs0V46/tPCK39k\nHPCNMUXAz4ArgShwozFmRKbnle7Rwqvs0F462aGFV/7wood/AbDKWrvGWrsfeBq4xoPzSjdo4VV2\naC+d7NDCK394EfBPBdZ1evxu+phkkRZeZYf20skOLbzyR1YXXk2bNq39fk1NDTU1Ndl8+VDQbpn+\nattLZ/aSJE/eGVUg8l342nNtbS21tbW+nNtYm9kv1BhTDUyz1k5KP/46YK213z/geTbT15JDc5od\nhn8vTpNNMnJAVL0iH02fDm++6f5f/OE0O4z8cZy1u5OMOiXc7dkYg7XWeHEuL1I6C4BhxpjBxpgy\n4Abgzx6cV7pBg7bZo1k6/tOgrT8yDvjW2hbgduBlIAk8ba1dlul5pXs0aJsdmqWTHRq09Ycn8/Ct\ntS9Za8+y1p5prf2eF+eU7tGgbXZolk52aNDWH1ppGyAdu2XqzeEXzdLJno7dMtWevaKAL9IN+sYr\nKWTaD1+kmzq+8SrXJRHpHvXwRURCQgFfRCQkFPBFuknz8KVQKeCLdIPm4UshU8APEH0Biv80Dz97\n9AUo3lPADwh9AUp2aB5+dugLUPyhgB8Q2ksnOzQPPzu0l44/FPADQnvpZE/HPHwFe79oLx1/KOAH\nhPbSyR7N0vGf9tLxhwJ+gGgvHf9plk72aC8d7yngi3SDZulIIVPAF+kGzdKRQqaAL9INkfIIX+kz\nk/Hbf87Mm2Yq3SAFRQFfpBucZoefbJvMvN5fYvJTk5XDl4KigB8gWmnrP+Xws0crbb2ngB8QWmmb\nHbGKGANKRmBsCWedfJZy+D7RSlt/KOAHhFbaZpFt/4/4RCtt/aGAHxBaaZsdiaYEG1uWY00LK7au\nUCDyiVba+kMBPyC00jY7Bp84mCJKwEJxUTGDeg3KdZECSStt/aGAHyBaaeu/Ne+voYX9YCDVkmLt\nzrW5LlJgaaWt9xTwRbph8ImDKaYULJQUl6iHLwVFAV+kG9TDl0KmgC/SDW4PXzl8KUwK+CLd0Lil\nkRTNYKC5pZllW5blukgiR00BX6Qbdu/ffdjHIvkso4BvjPmEMSZhjGkxxpznVaFE8tXufQr4Urgy\n7eEvBT4OzPGgLCJ5b9XWVV0er9y6MkclEem+kkx+2Fq7AsAYY7wpjkh+29m8s8tjZ6/2eJHCoRx+\ngGi3TP8t3bz0sI/FO9ot03tH7OEbY14B+nc+hLtz1Dettc9358WmTZvWfr+mpoaampru/LgcRttu\nmU0lSd56LKrl6D5Z9N6iwz4Wb7Ttlrl2WJJ4yNpzbW0ttbW1vpzbWJv5rn/GmNeAf7XWLjzMc6wX\nryUH17CugQmPTqSVFKVFpdRNqaN6YHWuixU4fb/fl617trrdHqCiZwWb/n1TbgsVQA3rGog/OpEW\ntWeMMVhrPUmbe5nSUR4/h7Rbpv+cZoe9e/d2OTa2cmyOShNs2i3TH5lOy7zWGLMOqAZeMMa86E2x\npLu0W6b/5q+fzwd80KVrc/Pom3NXoADTbpn+yHSWznPAcx6VRTLUsVtmrksSPE6zQ93bDe7oVTrg\nV0YqmXzm5JyWK8g6dsvMdUmCI6OALxIGTrPDuOnjPrSNwi+u+oV6nlJQFPBFjuC1v7/Gsq3LuqRy\nhvYeSs3pNTkrk8ix0Dz8ANE8fO+9u2MDNz31hS5fYVt5QiX1tyiv7DfNw/eeAn5AtM3Df7JkIvHH\n4nqTeGD9FocRP7iYD4re65K3X3DbAip7Vea2cAHXNg+/fpjas5cU8AMi0ZSgiSStJkXj5kZ9uXaG\n5s6FsVcl2F329/ZgP7DXQBZ8TsE+GxJNCd7dm8QWqT17SQE/IDQP3xvbdjlM+VYDH7/B4Ud3xBg1\nIEaJKWFY72HMv3W+gn2WaB6+PzxZaXtUL6SVtr77z287rN2b5Kffiiq/fAxmveZw7fNxmiNJzu4X\n5Y3b6gFIbk4S7affabY99XuHh/6Y5IVHw/27z9eVtpJjHfPww/vm6C6n2eEPCxq47kaHz/x7gv0n\nummEVdvdNEKkPEL1QP1Oc6FjHr5+915RwJfQ+ttahzO+E+f65ycyb0SceTMHE+sfpbRIaQQJJs3D\nl9BZvRp++EN44rUEH3wyCSbFtqJGNu9bS/2UeqVwJLDUww8QzcM/OKfZoWFdA39d6nDLLXD++RCJ\nwKJZMUYO6NqjVwonf2gevvc0aBsQTrPD8O/FabJJRg4I1/7hh+M0O4x5IM6qHUmKtkX5er96/nVq\nhN69O/5dPfr84zQ7jPxxnLW7k4w6JdztWYO28iGah9/V7t3w+ONw4bUJVm53B2KLKhq5+rPJ9mAP\nqEefpzQP3+X1pxsF/IDQPHxXIgFTp8Jpp8Hvfgf/8fkYo07RQGyh0Tx8N9hfOD3u6Tk1aBsQbfvh\nr92b5KdTwpWeaNrh8LPfJnj5yRjr3o5w662wcCEMHgwQ4apmDcQWmrb98B/6Y5IXvhGu67Z3L7zw\nAtz/hwTJYd5+slHAD5Aw7Ye/fz+8+io88YzD0z3jtJ6cZPBVURK319P7+K7BoS1tI4UlTPvht7a6\n23k88QQ8+yyMHg03firGjr1RErzl2eso4EvBaGmBOXPgmWfgD3+AM86AcdcnMHuT2NYU6/c1smJ7\nkurjFdwlvznNDommBGXvx3jumQhPPgk9e8LNN8PixW5KEiLc3FxPr9t7efa6CviS11pbYd48N8j/\n/vdQWQn/9E+wYAGcfjo4zTHmPBalcXNjaHO9UliWv+Nw2a/jbNifpGR7lH8urecPf4hwzjlgDpiL\n43UqSwFf8kZbryfaL8byJRGeeQZ++1s46SQ3yNfVwZlndv2ZSHlEi6Uk761c6X4q/eMfYdmuBLs+\nkYTiFKaikU9PSXJullKOCvgB4i68SuA0xwou8O3Y7TD2wTjvOEmKt0cZPLueG6+P8NJLED1Cp105\n+mByF14VTntu67DEKmKcUBZh8eKOIL9tG1x7LXznO3BedYxLn8zNp1ItvAqIQlx4tXUrzJoFM2bA\nC281sPO6iVCcosSUUjeljgtPUxAPq0JbeOU0O0x4NE7j5iQnpaL0eKqeMiJcdx18/OMwbhwUFXV9\n/tF+KtXCK/mQQlh4ZS0sWuT2csaPh6FD3ZTNxRfDG3+OMbrSnS8fragiVqFcfJgVysKrdevg4Yfh\nY7clWPJekpRNsa24ke8+kmTVKvjBD+DCC7sGe8jdgj+ldAKibeFVk82vwcudO2H2bLcX/+KLcMIJ\ncNVVcPfdMHEilLdPuVMuXjq0Lbxau7uRqv7505737nXHkmbNgpdegqYmuOIKuOHyGE27oqza7r7/\nrqmOfmgANh8opRMgd9yzgaX7ZvDwv1+Vk29mcpodlm5KULYjxpxXIsyY4c6mGT/eDfIf+ciHB11F\nDuWRZzbws5dnMOO/c9Oewf1UunKlG9xnzYLXX4dRo2DSJLjySjjvPCgudp/r175MXqZ0FPADIpc5\n/C1b4MW/OPzfpXG2FiUp2RHlM6l6rv1IhEsvheOPz0oxJEBykcNvG3Qd1CPGgtcj7UE+lXID/KRJ\ncNll7qyxbPIy4CulExAHy+H7NXNl+3b3Y+1rr7m31ashOinB9qokkML0a+S2Kf69vgTfwXL4frWn\nfftgzhsOn3ktzsaWJEVbo1z8Tj1XXx5h6lSoqvrw/PhCpYAfEH7m8HfuhPr6jgC/cqU7EHXJJfDQ\nQ+7+8ntbY8S1AEo84mcOf98+N9U4Zw7U1kJDA1RekKAp7s6NLx7QyHfvDGaHRSmdAMkkh995DnFR\nKsLcuR0BPpGAsWPh0kvdIH/BBVBWdvBzaNBVvOJVDr8twNfWurc33oDhw6Gmxr1NmAAlPR3ij8Xb\nOyz5NA00b3L4xpgfAB8FmoG3gSnW2p2HeK4Cvo8yyeG/t81hwiNxVn+Q5LhdUXiknvNiES65xA3w\nF14Ixx3ncwVEOjnWHL7T7LBwfYLmtTHefD1CbS3Mnw9nneVO/62pgXj84Hn4fO2w5FPA/wfgVWtt\nqzHme4C11t55iOcq4PuoYV0DFz0ax9JCiSmh/rP17R9JO/feI+URNm1yd+Z7/XU3VbN0RwPNN3Ys\nenr5xjouOTN4H2elcDSsa2DCo3FaD9KeD7Rjh5uWeXWuw4PNcT7okaTHriifK6rnipoIEyZkf6DV\nS3kzaGutnd3p4RvA9ZkVR47V4BMHYzFgwWIpKyqjYV0Dg3oN5h8en8zKHUlO3Belzx/r2fpehPHj\n3Y+yP/oRnDUyxhVPd+TfxwxS/l1yq3N7bqWVNTvWEO0XxVqYtTjB9hUxFja4qcfVq2HMGBg8IcHe\nE9yJA6nejdykiQMf4uWg7WeBpz08n3TDrLdnASkw0EILY6ePw7aCef90WnuthuIUO8sa+cmDSW6K\nVx+w8k+LniS/zHp7Fjbdnltp5YZnbyCy52x27y6m5aTlnLgvyp0j6rn1VneXydJSd+fUxZo4cFhH\nTOkYY14B+nc+BFjgm9ba59PP+SZwnrX2kD18pXT81ee/+rB93/aOAxYwUGxKGNL7dNbsWJN3g1Ei\nh/Kh9gwYiikyhhaborTI3W/pwB58vubhM5HVlI619vIjFOYWYDJw6ZHONW3atPb7NTU11NTUHOlH\n5Cgd+ObA0P49rjNvmsnanWsD9SaQYPtQewZGnDyc4qJiVmxdccgefBB2Tq2traW2ttaXc2c6aDsJ\nuA+YaK3deoTnqofvo8sevYxX173a/nj8KeO5b/J9CvJSkA5sz+f0O4e6W+sAAteDP5J8mqWzCigD\n2oL9G9baLx3iuQr4Pmt7k1x62qX85bN/yXVxRDKi9uzKm4DfrRdSwBcR6Tbthy8iIt2mgC8iEhIK\n+CIiIaGALyISEgr4IiIhoYAvIhISCvgiIiGhgC8iEhIK+CIiIaGALyISEgr4IiIhoYAvIhISCvgi\nIiGhgC8iEhIK+CIiIaGALyISEgr4IiIhoYAvIhISCvgiIiGhgC8iEhIK+CIiIaGALyISEgr4IiIh\noYAvIhISCvgiIiGhgC8iEhIK+CIiIaGALyISEhkFfGPM/zPGvGWMWWSMeckYM8CrgomIiLcy7eH/\nwFo72lp7LjADuMuDMhWk2traXBfBV0GuX5DrBqqfdMgo4Ftrd3V6eDzQmllxClfQG12Q6xfkuoHq\nJx1KMj2BMebbwGeAHcAlGZdIRER8ccQevjHmFWPMkk63pen/fxTAWvsf1tpBwK+BqX4XWEREjo2x\n1npzImNOA2Zaa0ce4t+9eSERkZCx1hovzpNRSscYM8xa+7f0w2uBZYd6rlcFFhGRY5NRD98Y83tg\nOO5g7RrgC9ba9zwqm4iIeMizlI6IiOQ331faGmMmGWOWG2NWGmPu8Pv1/GKMWd1pkdmb6WO9jTEv\nG2NWGGPwyo0fAAADr0lEQVRmGWNO7PT8O40xq4wxy4wxV+Su5AdnjHnEGLPJGLOk07Fu18cYc156\nEH+lMebH2a7HoRyifncZY941xixM3yZ1+reCqZ8xZqAx5lVjTDI9ieIr6eOBuH4Hqd/U9PGgXL9y\nY8z8dCxZaoy5K33c/+tnrfXthvsH5W/AYKAUWAyM8PM1fazLO0DvA459H/ha+v4dwPfS96uARbhj\nJKenfwcm13U4oOwTgHOAJZnUB5gPjE3fnwlcmeu6HaZ+dwFfPchzzy6k+gEDgHPS908AVgAjgnL9\nDlO/QFy/dFl6pv9fDLwBXJCN6+d3D/8CYJW1do21dj/wNHCNz6/pF8OHPxFdA/wyff+XuAPXAB8D\nnrbWpqy1q4FVuL+LvGGtfR3YfsDhbtUnvZVGxFq7IP28X3X6mZw6RP3AvY4HuoYCqp+1dqO1dnH6\n/i7cyRIDCcj1O0T9Tk3/c8FfPwBr7e703XLcQG7JwvXzO+CfCqzr9PhdOi5cobHAK8aYBcaYz6WP\n9bfWbgK3kQIV6eMH1ns9hVHvim7W51Tca9qmEK7v7caYxcaYhzt9ZC7Y+hljTsf9JPMG3W+PhVS/\n+elDgbh+xpgiY8wiYCPwSjpo+379tFvm0bvIWnseMBn4sjEmjvtHoLOgjYAHrT4PAEOttefgvtHu\ny3F5MmKMOQH4PfAv6Z5woNrjQeoXmOtnrW217h5kA3F761GycP38DvjrgUGdHg9MHys4Nj3d1Fq7\nGXgON0WzyRjTHyD98aop/fT1wGmdfrxQ6t3d+hRUPa21m2062QlMpyPNVnD1M8aU4AbDJ6y1f0of\nDsz1O1j9gnT92lhrdwK1wCSycP38DvgLgGHGmMHGmDLgBuDPPr+m54wxPdO9DYwxxwNXAEtx63JL\n+mn/B2h74/0ZuMEYU2aMGQIMA97MaqGPjqFrTrRb9Ul/7HzfGHOBMcbg7qn0J/JHl/qZrtt3Xwck\n0vcLsX6PAo3W2vs7HQvS9ftQ/YJy/YwxfdvSUcaYHsDluOMU/l+/LIxGT8IdZV8FfD1Xo+IZ1mEI\n7gyjRbiB/uvp432A2en6vQyc1Oln7sQdTV8GXJHrOhykTk8BG4BmYC0wBejd3foA56d/J6uA+3Nd\nryPU71fAkvS1fA43Z1pw9QMuAlo6tcmF6fdZt9tjgdUvKNdvZLpOi9P1+Wb6uO/XTwuvRERCQoO2\nIiIhoYAvIhISCvgiIiGhgC8iEhIK+CIiIaGALyISEgr4IiIhoYAvIhIS/wsQLERjBcQpjQAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107d88358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(res.x, res.y[0])\n",
    "plt.plot(res.x, res.y[0], '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the solver pass flat regions in very few number of steps. It would not be the case with any explicit method, which will have to make very small steps to overcome the stability issue. For example `method='RK45'` will make more than million steps on this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}