{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
     "grade": false,
     "grade_id": "cell-98f724ece1aacb67",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CDS: Numerical Methods Assignments\n",
    "\n",
    "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
    "\n",
    "- Solutions must be submitted via the Jupyter Hub.\n",
    "\n",
    "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
    "\n",
    "## Submission\n",
    "\n",
    "1. Name all team members in the the cell below\n",
    "2. make sure everything runs as expected\n",
    "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
    "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
    "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
    "6. submit your solutions  **in time (before the deadline)**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "team_members = \"Koen Vendrig, Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "7317367b26066eae942edeea4e7d4280",
     "grade": false,
     "grade_id": "cell-e9dc85dd9ad77a5e",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Ordinary Differential Equations - Initial Value Problems\n",
    "\n",
    "In the following you will implement your own easy-to-extend ordinary differential equation (ODE) solver, which can be used to solve first order ODE systems of the form\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = \\vec{f}(x, \\vec{y}),\n",
    "\\end{align*}\n",
    "\n",
    "for $x = x_0, x_1, \\dots, x_n$ with $x_i = i h$, step size $h$, and an initial condition of the form $\\vec{y}(x=0)$.\n",
    "The solver will be capable of using single-step as well as multi-step approaches.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "910692c2f4b93f251a3a128caf2b0c37",
     "grade": true,
     "grade_id": "cell-acb87410eaef9fe1",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "79e708efbac889db28e537bcf01e26c6",
     "grade": false,
     "grade_id": "cell-a27c54a734bf2232",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 1\n",
    "\n",
    "Implement the Euler method in a general Python function $\\text{integrator(x, y0, f, phi)}$, which takes as input a one-dimensional array $x$ of size $n+1$, the initial condition $y_0$, the *callable* function $f(x, y)$, and the integration scheme $\\Phi(x, y, f, i)$. It should return the approximated function $\\tilde{y}(x)$. \n",
    "\n",
    "The integration scheme $\\text{phi}$ is supposed to be a *callable* function as returned from another Python function $\\text{phi_euler(x, y, f, i)}$, where $i$ is the step number. In this way we will be able to easily extend the ODE solver to different methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "8396264309b023abf89a62200fa55379",
     "grade": true,
     "grade_id": "cell-1f9655333f2dc3cb",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def integrator(x, y0, f, phi):\n",
    "    \"\"\"\n",
    "    Numerically solves the initial value problem (IVP) given by ordinary differential\n",
    "    equation (ODE) f(x, y) = y' with initial value y0 using the integration scheme\n",
    "    provided by phi.\n",
    "\n",
    "    Args:\n",
    "        x:   size n + 1 numerical array of x values\n",
    "        y0:  an initial value to the function f\n",
    "        f:   a callable function with signature (x, y), with x and y the current state\n",
    "             of the system, and f such that f(x, y) = y' in the ODE\n",
    "        phi: a callable function with signature (x, y, f, i), with x and y the current state\n",
    "             of the system, i the step number, and f as above, representing the integration\n",
    "             scheme to use\n",
    "\n",
    "    Returns:\n",
    "        An n + 1 numerical array representing an approximate solution to y in y' = f(x, y)\n",
    "        given initial value y0, steps from x, and using integration scheme phi.\n",
    "    \"\"\"\n",
    "    \n",
    "    eta    = np.zeros(len(x))\n",
    "    eta[0] = y0\n",
    "    \n",
    "    for i in range(1, len(eta)):\n",
    "        h      = x[i] - x[i - 1]\n",
    "        eta[i] = eta[i - 1] + h*phi(x, eta, f, i)\n",
    "    \n",
    "    return eta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "35d6f653b74f2c8e2ff509890754df72",
     "grade": true,
     "grade_id": "cell-04ac3a0537092706",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_euler(x, y, f, i):\n",
    "    \"\"\"\n",
    "    Returns the integrator phi = f(x, y) for the Euler method to solve the ODE\n",
    "    y' = f(x, y).\n",
    "\n",
    "    Args:\n",
    "        x: size n + 1 numerical array of x values\n",
    "        y: size n + 1 numerical array of y values with estimates up to index i - 1\n",
    "        f: a callable function with signature (x, y), with x and y the current state\n",
    "           variables of the system, and f such that f(x, y) = y' in the ODE\n",
    "        i: the current step number to calculate phi for\n",
    "\n",
    "    Returns:\n",
    "        The integrator phi(x, y, f, i) = f(x, y).\n",
    "    \"\"\"\n",
    "    \n",
    "    return f(x[i - 1], y[i - 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "a300f75f62faab4732ebcac7bde5e18e",
     "grade": false,
     "grade_id": "cell-047f2dfab5489a88",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 2 \n",
    "\n",
    "Debug your implementation by applying it to the following ODE\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = y - x^2 + 1.0,\n",
    "\\end{align*}\n",
    "        \n",
    "with initial condition $y(x=0) = 0.5$. To this end, define the right-hand side of the ODE as a Python function $\\text{ODEF(x, y)}$, which in this case returns $f(x,y) = y - x^2 + 1.0$. You can then hand over the function $\\text{ODEF}$ to the argument $\\text{f}$ of your $\\text{integrator}$ function.\n",
    "\n",
    "Plot the solution you found with the Euler method together with the exact solution $y(x) = (x+1)^2 - 0.5 e^x$.\n",
    "\n",
    "Then implement a unit test for your $\\text{integrator}$ function using this system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "692bfda705027a7fa4452a7e03c8fa7d",
     "grade": true,
     "grade_id": "cell-8b8cb282dfe95a03",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAHmCAYAAAAhuHt9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABu9ElEQVR4nO3deXxU5dn/8c81kw3CJmEnIFARRQiLYVFcQQWrgm2tu6K02s2ltbXFblrb/h5bfarVupTHDeuGu7jUXYpaF0BRdkFkCWsACQTIOtfvj5nEISQhQDJnknzfr9e8Zs59lvs7mRCu3LnPOebuiIiIiIhIcggFHUBERERERL6mAl1EREREJImoQBcRERERSSIq0EVEREREkogKdBERERGRJKICXUREREQkiahAFxERERFJIirQRUSkWTKz4Wb2vpnNNLPHzCw16EwiIqACXUREmq/VwGh3Pw5YAUwINo6ISJQKdBGRfWRmK8zspAPYf4GZnVB/iYJlZqebWWbQOfaVu69z912xxRIgEmQeEZEKKtBFpNkys2PM7L9mVmBmW8zsPTMbVs997FHMu/sR7j6jPvvZ1wz17DYgpQGP36DM7GDgFOCFejjWFWY228yKzezBOmzf3syeNbMdZrbSzM4/0Awi0vg12h+oIiIHwszaAC8CPwKeANKAY4HiIHM1NmZ2GPCOuxcEnWV/xL4P/gVc4u6le9n2BgB3v6GWzdYCfwLGAi3qEOFOoqP3nYHBwEtm9qm7L6jDviLSRGkEXUSaq0MB3P0xdy93913u/pq7fwZgZoeb2Qwz2xqbkjK+pgOZmZvZIXHLD5rZn8zsX0BP4AUzKzSzX8bWV45o762f2La/MLPPYiP908wsI7buV2a2xsy2m9kSMxtTTbaaMtTp/ZlZKzMrN7OucW0DzGydmbUGRgP/rMsXvJav3976OJBj/9XMnotbvtnM3jSzNDNLAR4H/uDuSw6knwru/oy7PwdsrkO2TOA7wO/cvdDd3wWmAxfVRxYRabxUoItIc/U5UG5mU83sVDM7qGJF7GoeLwCvAZ2AK4FHzKzfvnTg7hcBq4Az3L2Vu/81fv0+9HM2MA7oDeQAl8S2uQIY5u6tiY7YrqhLhn15f+5eCCwGhsY13wT8P3ffDrR29w+q7mdmL8aK/+oeL+5jHwfiL8CJZjbEzH5I9Ov4bXcvAc4DRgC/i/2ycs4B9rWvDgXK3P3zuLZPgSMSnENEkowKdBFpltx9G3AM4MD/AflmNt3MOgMjgVbATe5e4u5vEZ0Oc149x6hrP7e7+1p330K0sB4MlAPpQH8zS3X3Fe7+RT33W2EWseLZzI4D+hMbNXf3v1S3g7uf7u7tanicvi99HAh33wzcCkwFrgO+WTEdx93/5e5Z7n5C7DHtQPvbR62AbVXaCoAD+quBiDR+KtBFpNly90Xufom7ZwMDgG5ET3jsBqx29/ireqwEutdzhLr2sz7u9U6glbsvA34K3ABsNLPHzaxbPfdbobJ4Bv5KdEpGSR37qqs69REb6fYaHu/WcOxPgIHAde6+el9Cxf8lAJgMTK7pLwH7oRBoU6WtDXCgfzUQkUZOBbqICODui4EHiRbqa4EeZhb/M7InsKaG3XcCLeOWu8QfupZu97WfqpkfdfdjgINj/VQ7ml1Nhn3tdxYw1My+A2QAj+4tm5n9OzbnvbrHv/e3j9hIt9XwOKaaHAOBu4mOoE/aW+5q+qv8SwDRaTc37eUvAfvicyDFzPrGtQ0CdIKoSDOnAl1EmiUzO8zMfm5m2bHlHkSneHwAfEi06P6lmaVa9JrlZxA9obA6c4HzzSxsZuOA4+PWbQD61LDfvvYTn7+fmY02s3SgCNhFzdfxrpphX/v9lOgvHf9LdBS6tl86AHD3U2Nz3qt7nFoffeyNmXUnOiXoh8CPgYHWwNefN7OU2Em8YSBsZhmxk1H34O47gGeAG80s08xGEb1Z0r8aMqOIJD8V6CLSXG0neoLgh2a2g2hhPh/4eWxqxRnAqcAm4C7g4tgoe3Wujm2/FbgAeC5u3f8Av41NifhF/E770U+8dKIjupuIToHpRHSOdXV2y7Cv/bp7MTAPWOHu1Y1+H7D67sOil098Gfibu093953AzcCfD/TYe/Fbor8sTQYujL3+bVyuf5vZr+O2/zHRyzFuBB4DfqRLLIqI1cMghYiINGFmlgYsA86u7ootjaUPEZHGQiPoIiKyN9cD7zVw4ZyIPkREGgUV6CIiUi0zG2pmBcBxRK+V3ij7EBFpbDTFRUREREQkiWgEXUREREQkiahAFxERERFJIirQRURERESSSLU3T2jOOnTo4L169Qo6hoiIiIg0cXPmzNnk7h2rtqtAr6JXr17Mnj076BgiIiIi0sSZ2crq2jXFRUREREQkiahAFxERERFJIirQRURERESSSJMv0M1snJktMbNlZjY56DwiIiIiIrVp0gW6mYWBO4FTgf7AeWbWP9hUIiIiIiI1a+pXcRkOLHP35QBm9jgwAVgYaKoqIuXlRCLlQcdIKLMm/bvhbsws6AgJ05zeK4CFms/3sYiIJE5TL9C7A6vjlvOAEQFlqdGcF//JsE+uCzqGiOyHcjcihHDACRHB8NgjguFW8ToEFW27bRNrt/j2UOV+FctfbxP9pcAttMf2YJXtVPYdwi1MxMLRdZYSfQ6lVC4TW/76OYyHwmBhLJQSex2CUAoWCkMoBULRdfHPoXDFciqhcAhLSSeckoalpJKSkk4oNZ1QahopqemkpKSRkpYefZ2aTkpaGqlp6aSmphMKhwP7PEVEkkFTL9DrxMwuBy4H6NmzZ8L7z/rGkby/5YcJ71cSwD3oBAnUnN4r4BHMHcfBI7EHmJdHX8e2waMleeU2lftGAAf3yte7tVHDNsRv70A0h8Wev96mHNwJeST6a4CXRx9ECFFOyL9+DhNtr3hO8YrlCKmW+L/ulXqYMsKUkkKZpVAW91xuqZRbCmWWSnkojbJQOuWhdMrDGXg4jUhKBh5Oh5QWkJKBpWZAagtCqRmE0loQSmtBOK0F4bSWpKZnkJLektT0lqS1yKRFZltaZLYmLS0t4e9ZRCReUy/Q1wA94pazY227cfcpwBSA3NzchFcZfQaMoM+ApBvYFxEBwCPllJeVUVZWSiRSRllZOV5eSllZGV5eRll5GZGyMjxSRnl5tK28vIxIaQllpSVEyoopLyvGy0opL4sue1kJXlaKl5dEX5eXQlkJREqhvATKS7FICZSXEYpE20ORUkKRsthzCWEvoUVZAamRYtK8mFQvIY0S0ryUDEoI2f79OC/2VHZaBkVkUBRqQUnsURZuQVlKJpHUlkRSMiEtE0/LJJTeCkvPJCWjNWmZB5Heqh0t2mTRsk17Wrdtr4JfRPZZUy/QZwF9zaw30cL8XOD8YCOJiDQuFgqTkhYmJS096Ch15pEIpaXFFO3aSUnxDkqKdlFatJPS4h2UFRdRVryD8pIiykt2EikpIlKyk0jxDrxkB5TsIFQafYTLdpJSvovU8p1klhSQXrSLDC8iw4vItKI6ZdnhGRRaS3aGWlEUyqQ4pTWlqa0pT2tNJL0tpLfBWh5EaqsOpLftRMt2nWjdvjNt23ckvRF9zUWk/jTpAt3dy8zsCuBVIAzc7+4LAo4lIiINzEIhUtNbkJreAshqkD48Uk7xrh3sLCygaMd2incWUFz4FcWFX1G24yvKdxXguwqgeBuh4m2klGwjtayQVmVbyCheTabvoJXvqHUaUYFnsi3UhsJwW4pS21GSdhBlGe2hRXusVUfS23UhM6s7bTpmk9WpG+karRdpEpp0gQ7g7i8DLwedQ0REmhYLhcnIbENGZpv9P4g7RTu3U7h1E4VfbWTn1g2UbMuntHATkR2bsZ1bSCneQlrJVtqWbCSzaCnttm4j3Ur3OFS5G/nWlq2h9uxIzaIooyPlLTtC6y6ktu1CRlY2B3XpQ8euPclIVyEvksyafIEuIiKStMwqi/wO3fvUbR93indtZ9umdWzLX8OOr9ZQ8tU6ItvWE9qxkbSifFqVbKL7ti84aOtWUiyy2+6lHmatteerlE4UZnSltFVXrG02aVkH07pzL7K69SErqxOhUPO6bKpIMlGBLiIi0piYkd6yDR17tqFjz361buqRcrZ/tYGvNqymMH8VRZtWUb51NSnb19Bi1zp67viMDtvfJnX97tNsCjyT9eGuFLTIprh1L0JZvWnZpS9ZPQ6ja3YvUlN0KUyRhqQCXUREpImyUJjWWd1ondWNmm4D4uVlFGxay5a1y9m+cQUlm1diX60go3AV2bsW06lwJinrIxA7g2uXp7Eq1IWv0ruzs00frNPhtMk+gm6H5NAhK6vZ3bBMpCGoQBcREWnGLJxC2849adu5+vuAeFkJm9d8waa8xexYv4zIpi9I376KzrtW0WXjR6RuLIf50W3XkcX6tF4Utu4DHQ8jM7s/XQ4ZStfOnVW4i+wD82Z1I5W9y83N9dmzZwcdQ0REJOl5WQmb8z4nf/mn7FqzgNCWpbTd/gVdy1aTQUnldnl0Ym1GX3a2P5y07oPo1G84B/fqq6ky0uyZ2Rx3z63arhF0ERER2S+WkkaHXgPo0GvA7isiEQrWfcGG5Z+yY9WnhDfOp3vhYrqu+S+htQ6z4CtvxarUPmxteziWnUvn/sfQ5xv9VLSLoBH0PWgEXUREpGGU7Sxg/dKP2fzFHCLrPqPt1kV0L/2SdKKXjdzo7ViRcTiFHQfTsvcIDh54DF06dtD0GGmyahpBV4FehQp0ERGRxPGyYjYsncPGRe9B3mw6FHxGt/K1QPTa7sutJ2vaDoVeo8gefBJ9Du6lS0BKk6ECvY5UoIuIiASrZNsm8ua/w7ZlH5CxfjYH75xHC4oBWE53VrcZSqTn0XTNGUPfQw4lrIJdGikV6HWkAl1ERCS5eFkJG5d8yIZ5b5Ka9196Fn5GJrsA+ILurGx3FOG+Yzh0xFi6dsgKOK1I3alAryMV6CIiIkmuvIz8ZbNZ/9kbpK38DwcXziWDEoo9lfkpR7Cl67G0yzmVgYNHkJGm62FI8lKBXkcq0EVERBoXL9lJ3qdvsvWzf3PQunfJLlsJwBrvwOK2xxLufzoDjxpHVttWAScV2Z0K9DpSgS4iItK4FW1aycqPXsCX/JveBR+RTglbPZPPWoygpO+pHHr0mfTs2inomCIq0OtKBbqIiEjT4cWFrJr1EoWfPk+PTTNp49sp9lRmpQ2jsO8EBpzwXbI7ad66BEMFeh2pQBcREWmiysvYuGAG+R89Sbe1r3JQ5Cu2ewvmtDiK0sO+xeATvk3HdpoGI4mjAr2OVKCLiIg0A5FyNs57g03vP0KPDW/Q2newxVsxp9WJpORezFGjRusEU2lwKtDrSAW6iIhIM1NWzNrZL7Jt9uP03jSDdEpYwsEs7X4mfU6cxOHfOFh3M5UGoQK9jlSgi4iINF+RnV+xYsZDpH72CD2KllDsKXyQNpJdAy7gqJPOom1mWtARpQlRgV5HKtBFREQEYPvKT1jz1hS6rXqBNr6d5d6Ned3P4fBxP+DQnl2DjidNgAr0OlKBLiIiIrspLSLvvUeJfPhPeu5azHZvwbutTiHzmB9x9PARpIRDQSeURkoFeh2pQBcREZGabFv6Putf/zu9N75OKmX8N3QkW4b8hBNPHk9mRmrQ8aSRUYFeRyrQRUREZG/KCtbx5at30nnRQ7TxAubSj5WHX86xp11I+1YZQceTRkIFeh2pQBcREZE6K9nJqjen0HLOXXQo28BSz2Zer0sZMf5yume1CTqdJLmaCvSkmzRlZjeb2WIz+8zMnjWzdnHrrjOzZWa2xMzGxrWPi7UtM7PJce29zezDWPs0M9Op1yIiIlJ/0lrS89Sf0uG6Bawf83daZ6Ty7ZV/pOz2XJ66/39Z/9WOoBNKI5R0BTrwOjDA3XOAz4HrAMysP3AucAQwDrjLzMJmFgbuBE4F+gPnxbYF+Atwq7sfAnwFfC+h70RERESah3AqXY69hC6/+pjNZzxIWovWnLXqRgpvG8ZTD93Bpu27gk4ojUjSFeju/pq7l8UWPwCyY68nAI+7e7G7fwksA4bHHsvcfbm7lwCPAxMsekeB0cBTsf2nAmcm6G2IiIhIcxQKkXXkt+j6y1lsOnUKrTNSOGv5b8m/ZQRPPnovBTtKgk4ojUDSFehVTAL+HXvdHVgdty4v1lZTexawNa7Yr2gXERERaVihEB1GnEPnX33CxpNup0NaKd/9/OcsufkEnv/3vyktjwSdUJJYIAW6mb1hZvOreUyI2+Y3QBnwSALyXG5ms81sdn5+fkN3JyIiIs1FKEynYybScfJnrBv1Jw4PreaMD87jjf85i//M/hRdrEOqkxJEp+5+Um3rzewS4HRgjH/9nbsG6BG3WXasjRraNwPtzCwlNooev33VPFOAKRC9iss+vRkRERGRvQmn0vXkK/FjLmT1839izOIHKXvhZJ6ccTYDzv4d/Xt2DjqhJJGkm+JiZuOAXwLj3X1n3KrpwLlmlm5mvYG+wEfALKBv7IotaURPJJ0eK+zfBs6K7T8ReD5R70NERESkKmtxED3P/V/sills7HI8Zxf+i1b3jeKRf/0f24tKg44nSSLpCnTgH0Br4HUzm2tm9wC4+wLgCWAh8ArwE3cvj42OXwG8CiwCnohtC/Ar4BozW0Z0Tvp9iX0rIiIiIntK7dCHXj96ksJznyMjI5MLvvgFs/5yGq+/P0fTXkQ3KqpKNyoSERGRhCorYd0rN9N+9t8pc3iu3USOvuC39O7UNuhk0sAazY2KRERERJqVlDS6nv4bUq78kC2dRnBBwRR23HkCT738GuURDaQ2RyrQRURERJJAOKs3PX48nYLT76VneAvjPzyXx//3Kpat+yroaJJgKtBFREREkoUZbXO/S+trZpOffTIX7HiIontO5MmXXtFoejOiAl1EREQkyVirjnS/bBoFZ9xLz/BXTPjofB7+289Zvbkw6GiSACrQRURERJJU2yO/S+ufz2FztxOYWHgfq27/Jq9++GnQsaSBqUAXERERSWKW2YGulz/NlhP/Qq4t5siXT+fe++6hsLgs6GjSQFSgi4iIiCQ7M9of/0PCP5yBZ3bm+6t/xSs3T2T+yvygk0kDUIEuIiIi0kikdOlPx5+9y/rDL+Wsshcpve9Unv/PR0HHknqmAl1ERESkMUnNoMs5t7F9/H0cFs5j1Fvf4Z8PPkBRaXnQyaSeqEAXERERaYRaDz2LtB/9B2uZxfe//BmP/u0aVm3aEXQsqQcq0EVEREQaqXCnfmT99F029RzHpF0Psvgf3+a/i1YHHUsOkAp0ERERkcYsvRWdJz3GV6Ou5yQ+pPVjZ/D0DM1Lb8xUoIuIiIg0dmYcdPI1FJ31MIeE1zPq7bP557RnKCuPBJ1M9oMKdBEREZEmouWA00m7/HUy0tO5eOEPuPueW9lWVBp0LNlHKtBFREREmpBw14G0u+oddhx0GFfm38ijt01mfUFR0LFkH6hAFxEREWlqWnWiw09eZ1OPsfyw6F5ev/2HLNuwPehUUkcq0EVERESaotQMOlz6GJsPv5CLyp9l4d0X8MmKjUGnkjpQgS4iIiLSVIXCZJ39D7aO+AXj+Q8FD5zDfxasDDqV7IUKdBEREZGmzIx2p/6O7SfdzLE2l8xp3+Xfsz8POpXUQgW6iIiISDPQ+pjLKTnzXgaHltFl+nm8+NHCoCNJDVSgi4iIiDQTLQZ/h/KzHmRgaAUHv3gez703L+hIUg0V6CIiIiLNSPqA8UTOeYR+oTX0e/V8npz5SdCRpAoV6CIiIiLNTNrh4+C8x+kT3sCgNy5g2ttzgo4kcZK2QDezn5uZm1mH2LKZ2e1mtszMPjOzoXHbTjSzpbHHxLj2I81sXmyf283MgngvIiIiIskmrd9JhC54il7hTQx8+1Ke1nSXpJGUBbqZ9QBOAVbFNZ8K9I09Lgfujm3bHrgeGAEMB643s4Ni+9wNXBa337hE5BcRERFpDFIPOQ7OfYS+obX0efUSXpqlq7skg6Qs0IFbgV8CHtc2AXjIoz4A2plZV2As8Lq7b3H3r4DXgXGxdW3c/QN3d+Ah4MyEvgsRERGRJJfW72T8rAfICS0n64WLef3TL4OO1OwlXYFuZhOANe7+aZVV3YHVcct5sbba2vOqaRcRERGROGkDzqB0/D0MDy0m/emJzFyUt/edpMEEUqCb2RtmNr+axwTg18DvE5zncjObbWaz8/PzE9m1iIiISFLIGHoOReNu5bjQp2x/7Pt8vHJz0JGarUAKdHc/yd0HVH0Ay4HewKdmtgLIBj42sy7AGqBH3GGyY221tWdX015dninunuvuuR07dqyfNykiIiLSyLQceSk7jvsdp4XeZ/4DV7M8vzDoSM1SUk1xcfd57t7J3Xu5ey+i01KGuvt6YDpwcexqLiOBAndfB7wKnGJmB8VODj0FeDW2bpuZjYxdveVi4PlA3piIiIhII5F54s/ZljOJi3mBl6b8lo3bi4KO1OwkVYG+Fy8THWFfBvwf8GMAd98C/BGYFXvcGGsjts29sX2+AP6d4MwiIiIijYsZbc68ha29TuXK0gd44J//S2FxWdCpmhWLXuBEKuTm5vrs2bODjiEiIiISrNIitk45jRYb53JL55v45Q++R2q4MY3tJj8zm+PuuVXb9VUWERERkT2lZtBu0lMUte7JTzZcz9+ffBUN7CaGCnQRERERqV6Lg2g76RnSUlOYsPDnPP7O/KATNQsq0EVERESkZu17k3H+I/QOrafbGz/hnSXrg07U5KlAFxEREZFahfocS9mpt3B86FO+fOwaXX6xgalAFxEREZG9yhgxie2DvsfFvMTT9/4PBTtLg47UZKlAFxEREZE6aT3+rxR0O5ari+7h7w89RiSik0Ybggp0EREREambcAptL3yIkpad+N66G7jvNV2auiGoQBcRERGRumvZnswLH6VTaDv9/vszZi7WSaP1TQW6iIiIiOwT6z6EyKl/5bjQPBY//mvWbN0VdKQmRQW6iIiIiOyz9GGXsP3wc7icp7nv/nsoLisPOlKToQJdRERERPadGa2//Xe2tT2Mqwpu5h/Pzgg6UZNRa4FuZhlmdpaZ/d3MnjSzh8zsl2Z2RKICioiIiEiSSm1Bm4sfpUXYOWber3l9/tqgEzUJNRboZvYH4D3gKOBD4J/AE0AZcJOZvW5mOQlJKSIiIiLJKesbhE7/X0aEFrP06RtYX1AUdKJGL6WWdR+5+/U1rPubmXUCejZAJhERERFpRFKHns/2xa/zgyVPcONDI/j9TyYRDlnQsRqtGkfQ3f0liE5zqbrOzDq4+0Z318UvRURERITW3/47uzK7c9mm/8f9b84NOk6jVpeTRGeZ2ciKBTP7DvDfhoskIiIiIo1ORhsyz59KV/uKbjMn8/HKLUEnarTqUqCfD9xhZjeb2SPAZcDoho0lIiIiIo2NZedSevyvOS38AS89ejs7S8qCjtQo7bVAd/d5wJ+BHwInAle4e15DBxMRERGRxifj+J+xveORXFX0T+6a/k7QcRqlvRboZnYf8FMgB7gUeNHMftLAuURERESkMQqFaX3u/9EyVM7QT2/gvaX5QSdqdOoyxWUecKK7f+nurwIjgKENG0tEREREGq2sb8DJf2B0eC4zp93K9qLSoBM1KnWZ4nKbu3vccoG7f69hY4mIiIhIY5Y68gds73oUV5Tezz+eeTvoOI1KbTcqesHMzjCz1GrW9TGzG81sUsPGExEREZFGKRSi9dn/JC1sHLfoBt5etD7oRI1GbSPolwHHAovMbJaZvWxmb5nZcqJ3FZ3j7vcnJKWIiIiIND4HHUxo3J8ZFV7A+0/dpqkudVTbjYrWu/svgduBy4E/AtcAA9z9ZHd/vqFCmdmVZrbYzBaY2V/j2q8zs2VmtsTMxsa1j4u1LTOzyXHtvc3sw1j7NDNLa6jMIiIiIrKn1GGXsr3LSH5SNpV7XtStdOqiLieJdgKeBH4GdAF2NWQgMzsRmAAMcvcjgFti7f2Bc4EjgHHAXWYWNrMwcCdwKtAfOC+2LcBfgFvd/RDgK0Bz50VEREQSyYzWZ91JZqiMwz79f3y86qugEyW9upwk+lugL3AfcAmw1Mz+n5l9o4Ey/Qi4yd2LY/1vjLVPAB5392J3/xJYBgyPPZa5+3J3LwEeByaYmRG9odJTsf2nAmc2UGYRERERqUmHQyg/9uecEf6AZx+/n9LySNCJklpdRtCJXcVlfexRBhwEPBU//aQeHQocG5ua8h8zGxZr7w6sjtsuL9ZWU3sWsNXdy6q0i4iIiEiCpR93DYVtDuEHO+7igbfnBx0nqdXlRkVXm9kc4K/Ae8BAd/8RcCTwnf3p1MzeMLP51TwmAClAe2AkcC3wRGw0vMGY2eVmNtvMZufn62L6IiIiIvUuJY1WZ91Ftm0ideb/8OWmHUEnSlp1GUFvD3zb3ce6+5PuXgrg7hHg9P3p1N1PcvcB1TyeJzrS/YxHfQREgA7AGqBH3GGyY201tW8G2plZSpX26vJMcfdcd8/t2LHj/rwlEREREdmbniPYOegSLg79m3ufeIa4W+1InLrMQb/e3VfWsG5R/UfiOeBEADM7FEgDNgHTgXPNLN3MehOdF/8RMAvoG7tiSxrRE0mnx6blvA2cFTvuRKDBrjwjIiIiInvXctwfKEk7iG+v/zsvf7Y26DhJqU5z0BPsfqCPmc0nesLnxNho+gLgCWAh8ArwE3cvj80xvwJ4FVgEPBHbFuBXwDVmtozonPT7EvxeRERERCRei3akj/sjR4aWMueFu9lZUrb3fZoZ058Wdpebm+uzZ88OOoaIiIhI0xWJUHj3aHZtXM6jw5/h6tNyg04UCDOb4+57vPlkHEEXERERkaYsFKLVt24ly7bR5oP/1QmjVahAFxEREZHE6zaE4oEXclHoFf7v6Zd0wmgcFegiIiIiEogW4/5AWWorTs+7lTcXbgg6TtJQgS4iIiIiwcjMIvXk33F0eCEzXphKSZnuMAoq0EVEREQkQOHcSexo8w0u3fkAj76/LOg4SUEFuoiIiIgEJ5xCy9P+H98IrWP9W3dTsLM06ESBU4EuIiIiIoGyQ8eyo9soLo88wf+9/nHQcQKnAl1EREREgmVG5hk30c520Hb2HazavDPoRIFSgS4iIiIiweuaQ3H/s5kY+jf3vvh20GkCpQJdRERERJJCi3E3YOEww5bdzpyVXwUdJzAq0EVEREQkObTpho/8CWeEP2DaCy8225sXqUAXERERkaSRduzVFKe0YdyGe3ln6aag4wRCBbqIiIiIJI8W7Qgf9zNGh+fy0ovPNMtRdBXoIiIiIpJUUkb+gF3pHfj21vt5df66oOMknAp0EREREUkuaZmkj/4VI0KLmfHvaZRHmtcougp0EREREUk6oSMvYWfL7pxfOJXnPs4LOk5CqUAXERERkeSTkkbGSb8mJ/Qlc157iJKySNCJEkYFuoiIiIgkpdCgc9nRpg8Tix7liVkrg46TMCrQRURERCQ5hVNoedJ19AvlseitR5vNKLoKdBERERFJWjbgO+xs3ZsLiqfx7Merg46TECrQRURERCR5hcK0GH0t/UMr+eSNxygtb/qj6CrQRURERCSpWc7Z7MzswXlF03j+kzVBx2lwKtBFREREJLmFU2kx+loGhZYz640nmvx10ZOuQDezwWb2gZnNNbPZZjY81m5mdruZLTOzz8xsaNw+E81saewxMa79SDObF9vndjOzIN6TiIiIiBwYG3Qeu1p245ydj/Lip017FD3pCnTgr8Af3H0w8PvYMsCpQN/Y43LgbgAzaw9cD4wAhgPXm9lBsX3uBi6L229cYt6CiIiIiNSrlDTST/gFQ0PLeO/1p5v0KHoyFugOtIm9bgusjb2eADzkUR8A7cysKzAWeN3dt7j7V8DrwLjYujbu/oG7O/AQcGYi34iIiIiI1J/Q0AvZldGJCdun8dqC9UHHaTDJWKD/FLjZzFYDtwDXxdq7A/HX1smLtdXWnldNu4iIiIg0RinppB9zJaPCC3jjzVeIjsE2PYEU6Gb2hpnNr+YxAfgR8DN37wH8DLgvAXkuj813n52fn9/Q3YmIiIjIfgrlXkJJSitO3Pw4H325Jeg4DSKQAt3dT3L3AdU8ngcmAs/ENn2S6LxygDVAj7jDZMfaamvPrqa9ujxT3D3X3XM7dux4oG9PRERERBpKRhtCw77HqeGPeObNd4NO0yCScYrLWuD42OvRwNLY6+nAxbGruYwECtx9HfAqcIqZHRQ7OfQU4NXYum1mNjJ29ZaLgecT+k5EREREpN6lHPUjsDBHrPwXS9ZvDzpOvUvGAv0y4H/N7FPg/xG9YgvAy8ByYBnwf8CPAdx9C/BHYFbscWOsjdg298b2+QL4d4Leg4iIiIg0lDZdKRtwNt8N/4dH3vo46DT1LiXoAFW5+7vAkdW0O/CTGva5H7i/mvbZwID6zigiIiIiwUo/7qcw7xE6LJrKuoLhdG3bIuhI9SYZR9BFRERERGrX8VB29hnLhaFX+dd/Fgadpl6pQBcRERGRRqnlCT+nvRVSMvththWVBh2n3qhAFxEREZHGqecIdnYawnn8myc+Whl0mnqjAl1EREREGq2Wx/yEb4TWseS95yiPNI0bF6lAFxEREZHGq/8EijI6ctrO6by1eGPQaeqFCnQRERERabxS0kgdeTknhD/l1Rn/CTpNvVCBLiIiIiKNWjj3UsotlYFrn2Dx+m1BxzlgKtBFREREpHFr1ZGyI77DWeGZTHtnXtBpDpgKdBERERFp9NJH/ZhMKybts8f4akdJ0HEOiAp0EREREWn8ug5iZ9cRXGCv8PhHXwad5oCoQBcRERGRJqHlMT+mZyif5e89Q1l5JOg4+00FuoiIiIg0DYedRlFGJ04tfoUZS/KDTrPfUoIO0BiUlpaSl5dHUVFR0FGkkcrIyCA7O5vU1NSgo4iIiDRd4VRSh13CCe/czC/f+5CT+o8POtF+UYFeB3l5ebRu3ZpevXphZkHHkUbG3dm8eTN5eXn07t076DgiIiJNWjh3IpF3bqH3yqdYveUkerRvGXSkfaYpLnVQVFREVlaWinPZL2ZGVlaW/gIjIiKSCG2zKe5zMmeHZ/DEh18EnWa/qECvIxXnciD0/SMiIpI4LUZ+n45WwMZZz1JS1vhOFlWB3kiEw2EGDx5c+bjppptq3f7BBx/kiiuuSFC6Pc2YMYP//ve/lcuXXHIJTz311H4f70D3h7p9TVasWMGjjz5auTx79myuuuqqA+pXREREEuyQMezK7M740ld4feGGoNPsM81BbyRatGjB3LlzG+z4ZWVlpKTU37fDjBkzaNWqFUcffXS9HTMRKgr0888/H4Dc3Fxyc3MDTiUiIiL7JBQmffgkRr39R6559z1Oyzkr6ET7RCPojVyvXr3YtGkTEB3tPeGEE/bYJj8/n+985zsMGzaMYcOG8d577wFwww03cNFFFzFq1Cguuuii3faZMWMGxx9/PBMmTKBPnz5MnjyZRx55hOHDhzNw4EC++OKLGo+9YsUK7rnnHm699VYGDx7MO++8A8DMmTM5+uij6dOnT+VouLtz7bXXMmDAAAYOHMi0adMq26+44gr69evHSSedxMaNG6t9/7fffjv9+/cnJyeHc889F4AtW7Zw5plnkpOTw8iRI/nss8/22K/qiHyrVq0AmDx5Mu+88w6DBw/m1ltvZcaMGZx++um1HveGG25g0qRJnHDCCfTp04fbb799bx+biIiINLDQ0IsotzCHrX2GL/ILg46zTzSCvo/+8MICFq7dVq/H7N+tDdefcUSt2+zatYvBgwdXLl933XWcc845dTr+1Vdfzc9+9jOOOeYYVq1axdixY1m0aBEACxcu5N1336VFixZ77Pfpp5+yaNEi2rdvT58+ffj+97/PRx99xN///nfuuOMObrvtthqP/cMf/pBWrVrxi1/8AoD77ruPdevW8e6777J48WLGjx/PWWedxTPPPMPcuXP59NNP2bRpE8OGDeO4447j/fffZ8mSJSxcuJANGzbQv39/Jk2atEfGm266iS+//JL09HS2bt0KwPXXX8+QIUN47rnneOutt7j44ovr/NeHm266iVtuuYUXX3wRiP6iUqG24y5evJi3336b7du3069fP370ox/pkooiIiJBat2Z0r6n8d0lb3H3+0v59fghQSeqMxXojcSBTHF54403WLhwYeXytm3bKCyM/iY5fvz4aotzgGHDhtG1a1cAvvGNb3DKKacAMHDgQN5+++29HruqM888k1AoRP/+/dmwITof7N133+W8884jHA7TuXNnjj/+eGbNmsXMmTMr27t168bo0aOrPWZOTg4XXHABZ555JmeeeWblMZ9++mkARo8ezebNm9m27cB/qartuKeddhrp6emkp6fTqVMnNmzYQHZ29gH3KSIiIvsvY8QkMj6fzrZPnqPkm4NIS2kck0dUoO+jvY10J1pKSgqRSPTs5Jou4xeJRPjggw/IyMjYY11mZmaNx05PT698HQqFKpdDoRBlZWV7PXZtx3P3vW5fFy+99BIzZ87khRde4M9//jPz5s2r037xX7dIJEJJSckB5Yh/b+FwuPLrIyIiIgHqfTxFmd355rY3eXPRDzl1YNegE9VJ4/g1QmrUq1cv5syZA1A5ulvVKaecwh133FG5XJ8nm9Z07NatW7N9+/a97n/ssccybdo0ysvLyc/PZ+bMmQwfPpzjjjuusn3dunWVI/bxIpEIq1ev5sQTT+Qvf/kLBQUFFBYWcuyxx/LII48A0SkqHTp0oE2bNrvtG/91mz59OqWlpXvNXZfjioiISBIJhUg78kKOCc/njQ/mBJ2mzgIp0M3su2a2wMwiZpZbZd11ZrbMzJaY2di49nGxtmVmNjmuvbeZfRhrn2ZmabH29Njystj6Xgl7gw2gYg56xWPy5OiX4Prrr+fqq68mNzeXcDhc7b633347s2fPJicnh/79+3PPPffUW66ajn3GGWfw7LPP7naSaHW+9a1vkZOTw6BBgxg9ejR//etf6dKlC9/61rfo27cv/fv35+KLL+aoo47aY9/y8nIuvPBCBg4cyJAhQ7jqqqto164dN9xwA3PmzCEnJ4fJkyczderUPfa97LLL+M9//sOgQYN4//33K/+SkJOTQzgcZtCgQdx666277VOX44qIiEhyCQ05nxBOt5XPsWFb47hpoNXXVIN96tTscCAC/BP4hbvPjrX3Bx4DhgPdgDeAQ2O7fQ6cDOQBs4Dz3H2hmT0BPOPuj5vZPcCn7n63mf0YyHH3H5rZucC33H2vZ1Xm5ub67Nmzd2tbtGgRhx9++IG/cWnW9H0kIiISjF33nsbGVZ/z0okv8eMTD937DgliZnPcfY/rOQcygu7ui9x9STWrJgCPu3uxu38JLCNarA8Hlrn7cncvAR4HJlj09oyjgYrr5U0Fzow7VsUQ51PAGNPtHEVERESanRbDLubg0EY+/+i1ejsPriEl2xz07sDquOW8WFtN7VnAVncvq9K+27Fi6wti24uIiIhIc3L4GZSmtOLYwleYs/KroNPsVYMV6Gb2hpnNr+YxoaH63F9mdrmZzTaz2fn5+UHHEREREZH6lNYSBnyHU0Mf8fyHi4NOs1cNdplFdz9pP3ZbA/SIW86OtVFD+2agnZmlxEbJ47evOFaemaUAbWPbV5d1CjAFonPQ9yO3iIiIiCSx1CMvInXuVGzBs+woHk5mevJebTzZprhMB86NXYGlN9AX+IjoSaF9Y1dsSQPOBaZ7dBLR28BZsf0nAs/HHWti7PVZwFveGCYdiYiIiEj9y85lV9tDGM8MXp63Lug0tQrqMovfMrM84CjgJTN7FcDdFwBPAAuBV4CfuHt5bHT8CuBVYBHwRGxbgF8B15jZMqJzzO+Ltd8HZMXarwEqL80oIiIiIs2MGRnDLiY39DnvffRh0GlqFdRVXJ5192x3T3f3zu4+Nm7dn939G+7ez93/Hdf+srsfGlv357j25e4+3N0PcffvuntxrL0otnxIbP3yxL7LxuO5555j4cKFB3ycGTNmcPrpp9e6zdatW7nrrrsql9euXctZZ51Vyx4iIiIi9cMGnoVjHLz2ZdZu3RV0nBol2xQXCUB9Feh1UbVA79atG0899VQte4iIiIjUk7bdKc4+mgmh93j+kzV73z4gKtAbiYcffpjhw4czePBgfvCDH1BeXs6sWbPIycmhqKiIHTt2cMQRRzB//nwKCwsZM2YMQ4cOZeDAgTz//POVx3nooYcq79x50UUX8d///pfp06dz7bXXMnjwYL744ovd+n3yyScZMGAAgwYN4rjjjgOgqKiISy+9tPIOnm+//fYeeW+44QZuueWWyuUBAwawYsUKJk+ezBdffMHgwYO59tprWbFiBQMGDKj1uA8++CDf/va3GTduHH379uWXv/xlvX99RUREpHnIGHoefULrWTj77aS9Jnrynr6arP49GdbPq99jdhkIp95U4+pFixYxbdo03nvvPVJTU/nxj3/MI488wsUXX8z48eP57W9/y65du7jwwgsZMGAAZWVlPPvss7Rp04ZNmzYxcuRIxo8fz8KFC/nTn/7Ef//7Xzp06MCWLVto374948eP5/TTT692qsmNN97Iq6++Svfu3dm6dSsAd955J2bGvHnzWLx4Maeccgqff/55nd7qTTfdxPz585k7dy4AK1asqFxX23Hnzp3LJ598Qnp6Ov369ePKK6+kR48e1fQgIiIiUov+4yl/8RqGFrzOonXn0L9bm6AT7UEFeiPw5ptvMmfOHIYNGwbArl276NSpEwC///3vGTZsGBkZGdx+++0AuDu//vWvmTlzJqFQiDVr1rBhwwbeeustvvvd79KhQwcA2rdvv9e+R40axSWXXMLZZ5/Nt7/9bQDeffddrrzySgAOO+wwDj744DoX6LWp7bhjxoyhbdu2APTv35+VK1eqQBcREZF9l9GW8kPGcsaS//B/n6ykf7eBQSfagwr0fVXLSHdDcXcmTpzI//zP/+yxbvPmzRQWFlJaWkpRURGZmZk88sgj5OfnM2fOHFJTU+nVqxdFRUX71fc999zDhx9+yEsvvcSRRx7JnDlz6rRfSkoKkUikcnl/+6+Qnp5e+TocDlNWVlbL1iIiIiI1SxtyLh0+f4H1n7xC+akDCIcs6Ei70Rz0RmDMmDE89dRTbNy4EYAtW7awcuVKAH7wgx/wxz/+kQsuuIBf/epXABQUFNCpUydSU1N5++23K7cdPXo0Tz75JJs3b648DkDr1q3Zvn17tX1/8cUXjBgxghtvvJGOHTuyevVqjj32WB555BEAPv/8c1atWkW/fv12269Xr158/PHHAHz88cd8+eWXe+2rLscVEREROWB9T6YktQ3HF8/gg+XV3scyUCrQG4H+/fvzpz/9iVNOOYWcnBxOPvlk1q1bx0MPPURqairnn38+kydPZtasWbz11ltccMEFzJ49m4EDB/LQQw9x2GGHAXDEEUfwm9/8huOPP55BgwZxzTXXAHDuuedy8803M2TIkD1OEr322msZOHAgAwYM4Oijj2bQoEH8+Mc/JhKJMHDgQM455xwefPDB3Ua4Ab7zne+wZcsWjjjiCP7xj39w6KGHApCVlcWoUaMYMGAA11577W771OW4IiIiIgcsJZ3QgG8xLjSLl+YsCzrNHixZz14NSm5urs+ePXu3tkWLFnH44YcHlEiaCn0fiYiIJJEV78GD3+T6lJ/y+1/fEMg0FzOb4+65Vds1B11EREREmp+eR1HeOpvfd1qgOegiIiIiIoELhQgPOY9wOAUi5UGn2Y1G0EVERESkeTrxN2DJNXoOGkGvM83VlwOh7x8REZEklITFOahAr5OMjAw2b96sIkv2i7uzefNmMjIygo4iIiIijYCmuNRBdnY2eXl55OfnBx1FGqmMjAyys7ODjiEiIiKNgAr0OkhNTaV3795BxxARERGRZkBTXEREREREkogKdBERERGRJKICXUREREQkiZiuTLI7M8sHVgbQdQdgUwD9SmLpc24e9Dk3D/qcmz59xs1DkJ/zwe7esWqjCvQkYWaz3T036BzSsPQ5Nw/6nJsHfc5Nnz7j5iEZP2dNcRERERERSSIq0EVEREREkogK9OQxJegAkhD6nJsHfc7Ngz7npk+fcfOQdJ+z5qCLiIiIiCQRjaCLiIiIiCQRFehJwMzGmdkSM1tmZpODziP1z8x6mNnbZrbQzBaY2dVBZ5KGYWZhM/vEzF4MOos0DDNrZ2ZPmdliM1tkZkcFnUnqn5n9LPbzer6ZPWZmGUFnkgNnZveb2UYzmx/X1t7MXjezpbHng4LMCCrQA2dmYeBO4FSgP3CemfUPNpU0gDLg5+7eHxgJ/ESfc5N1NbAo6BDSoP4OvOLuhwGD0Ofd5JhZd+AqINfdBwBh4NxgU0k9eRAYV6VtMvCmu/cF3owtB0oFevCGA8vcfbm7lwCPAxMCziT1zN3XufvHsdfbif6H3j3YVFLfzCwbOA24N+gs0jDMrC1wHHAfgLuXuPvWQENJQ0kBWphZCtASWBtwHqkH7j4T2FKleQIwNfZ6KnBmIjNVRwV68LoDq+OW81Dh1qSZWS9gCPBhwFGk/t0G/BKIBJxDGk5vIB94IDaV6V4zyww6lNQvd18D3AKsAtYBBe7+WrCppAF1dvd1sdfrgc5BhgEV6CIJZWatgKeBn7r7tqDzSP0xs9OBje4+J+gs0qBSgKHA3e4+BNhBEvw5XOpXbA7yBKK/kHUDMs3swmBTSSJ49PKGgV/iUAV68NYAPeKWs2Nt0sSYWSrR4vwRd38m6DxS70YB481sBdGpaqPN7OFgI0kDyAPy3L3iL2BPES3YpWk5CfjS3fPdvRR4Bjg64EzScDaYWVeA2PPGgPOoQE8Cs4C+ZtbbzNKInoQyPeBMUs/MzIjOWV3k7n8LOo/UP3e/zt2z3b0X0X/Hb7m7RtyaGHdfD6w2s36xpjHAwgAjScNYBYw0s5axn99j0MnATdl0YGLs9UTg+QCzANE/1UmA3L3MzK4AXiV6lvj97r4g4FhS/0YBFwHzzGxurO3X7v5ycJFEZD9dCTwSG1RZDlwacB6pZ+7+oZk9BXxM9Cpcn5CEd5uUfWdmjwEnAB3MLA+4HrgJeMLMvgesBM4OLmGU7iQqIiIiIpJENMVFRERERCSJqEAXEREREUkiKtBFRERERJKICnQRERERkSSiAl1EREREJImoQBcRERERSSIq0EVEREREkogKdBER2S9mNszMPjOzDDPLNLMFZjYg6FwiIo2dblQkIiL7zcz+BGQALYA8d/+fgCOJiDR6KtCr6NChg/fq1SvoGCIiIiLSxM2ZM2eTu3es2p4SRJhk1qtXL2bPnh10DBERERFp4sxsZXXtmoMuIiIiIpJEVKCLiIiIiCQRFegiInXUa/JLXP/8/KBjiIhIE6c56CIi+2Dq+yv5wwRdSVCkMSotLSUvL4+ioqKgo0gzk5GRQXZ2NqmpqXXaXgW6iIiINAt5eXm0bt2aXr16YWZBx5Fmwt3ZvHkzeXl59O7du077aIqLiIiINAtFRUVkZWWpOJeEMjOysrL26S83KtBFRESk2VBxLkHY1+87FegiIiIiIklEBbqIiIhIkmnVqlWt67du3cpdd91Vubx27VrOOuusBsuzYsUKBgw48BPkb7jhBm655ZZat5k7dy4vv/zyAfe1L775zW+ydevWPdrrkrchqEAXERERSTB3JxKJ7Pf+VQv0bt268dRTT9VHtMAdaIFeVla2z/u8/PLLtGvXbr/7rG8q0EVEREQSYMWKFfTr14+LL76YAQMGsHr1am6++WaGDRtGTk4O119//R77FBYWMmbMGIYOHcrAgQN5/vnnAZg8eTJffPEFgwcP5tprr91thHvkyJEsWLCg8hgnnHACs2fPZseOHUyaNInhw4czZMiQymPVpT+IFr4XXHABhx9+OGeddRY7d+6szNK/f39ycnL4xS9+UfleR48eTU5ODmPGjGHVqlV79FWRC2DTpk306tWLkpISfv/73zNt2jQGDx7MtGnT6pR7xowZHHvssYwfP57+/ftTXl7OtddeW/m1/ec//wnAunXrOO644xg8eDADBgzgnXfeAaBXr15s2rQJgD//+c8ceuihHHPMMSxZsqTWvECNfR0IXWZRREREmp0/vLCAhWu31esx+3drw/VnHFHrNkuXLmXq1KmMHDmS1157jaVLl/LRRx/h7owfP56ZM2dy3HHHVW6fkZHBs88+S5s2bdi0aRMjR45k/Pjx3HTTTcyfP5+5c+cC0YK4wjnnnMMTTzzBH/7wB9atW8e6devIzc3l17/+NaNHj+b+++9n69atDB8+nJNOOonMzMy99gewZMkS7rvvPkaNGsWkSZO46667uPTSS3n22WdZvHgxZlY5TeTKK69k4sSJTJw4kfvvv5+rrrqK5557bq9fw7S0NG688UZmz57NP/7xD4A65Qb4+OOPmT9/Pr1792bKlCm0bduWWbNmUVxczKhRozjllFN45plnGDt2LL/5zW8oLy+v/CWjwpw5c3j88ceZO3cuZWVlDB06lCOPPLLWzPfdd1+1fdX1korV0Qi6iIiISIIcfPDBjBw5EoDXXnuN1157jSFDhjB06FAWL17M0qVLd9ve3fn1r39NTk4OJ510EmvWrGHDhg219nH22WdXTnd54oknKuemv/baa9x0000MHjyYE044gaKioj1Gtmvrr0ePHowaNQqACy+8kHfffZe2bduSkZHB9773PZ555hlatmwJwPvvv8/5558PwEUXXcS7776731+zuuQGGD58eGVR/Nprr/HQQw8xePBgRowYwebNm1m6dCnDhg3jgQce4IYbbmDevHm0bt16t2O88847fOtb36Jly5a0adOm8peTveWrrq8DkfARdDPrWcdNt7p7/f5qKyIiIgJ7HeluKPGjvu7Oddddxw9+8IMat3/kkUfIz89nzpw5pKam0qtXr71eT7t79+5kZWXx2WefMW3aNO65557K/p5++mn69eu3X/1VvVSgmZGSksJHH33Em2++yVNPPcU//vEP3nrrrb1+HQBSUlIq5+HX9p7qkhv2/NrecccdjB07do/tZs6cyUsvvcQll1zCNddcw8UXX3xAeWvra38FMYI+FXgw9lzT40HgzACyiYiIiCTE2LFjuf/++yksLARgzZo1bNy4cbdtCgoK6NSpE6mpqbz99tusXLkSgNatW7N9+/Yaj33OOefw17/+lYKCAnJycir7u+OOO3B3AD755JM99qupP4BVq1bx/vvvA/Doo49yzDHHUFhYSEFBAd/85je59dZb+fTTTwE4+uijefzxx4Fo0X/sscfu0VevXr2YM2cOwG4nuFZ9b3XJXdXYsWO5++67KS0tBeDzzz9nx44drFy5ks6dO3PZZZfx/e9/n48//ni3/Y477jiee+45du3axfbt23nhhRf2mremvg5EwkfQ3f3ERPcpIiIikmxOOeUUFi1axFFHHQVEL6348MMP06lTp8ptLrjgAs444wwGDhxIbm4uhx12GABZWVmMGjWKAQMGcOqpp/KTn/xkt2OfddZZXH311fzud7+rbPvd737HT3/6U3JycohEIvTu3ZsXX3xxt/1q6g+gX79+3HnnnUyaNIn+/fvzox/9iIKCAiZMmEBRURHuzt/+9jcA7rjjDi699FJuvvlmOnbsyAMPPLDH+//FL37B2WefzZQpUzjttNMq20888cTKKS3XXXddnXJX9f3vf58VK1YwdOhQ3J2OHTvy3HPPMWPGDG6++WZSU1Np1aoVDz300G77DR06lHPOOYdBgwbRqVMnhg0btte8NfV1IKzitxGJys3N9YozdEVE4vWa/BIAK246bS9bikgyWrRoEYcffnjQMaSZqu77z8zmuHtu1W2T6iRRM7s06AwiIiIiIkFKqgId+EPQAUREREREghTEVVw+q2kV0DmRWURERKR5cfc9rkYi0tD2dUp5EDcq6gyMBb6q0m7AfxMfR0RERJqDjIwMNm/eTFZWlop0SRh3Z/PmzWRkZNR5nyAK9BeBVu4+t+oKM5uR8DQiIiLSLGRnZ5OXl0d+fn7QUaSZycjIIDs7u87bB3GZxe/Vsu78RGYRERGR5iM1NfWAbr8ukijJdpKoiIiIiEiz1qgLdDMbZ2ZLzGyZmU2uZv01ZrbQzD4zszfN7OAgcoqIiIiI1FXgBbqZnbGf+4WBO4FTgf7AeWbWv8pmnwC57p4DPAX89UCyioiIiIg0tMALdODP+7nfcGCZuy939xLgcWBC/Abu/ra774wtfgDUfXa+iIiIiEgAkqFA39/rHHUHVsct58XaavI94N/72ZeIiIiISEIEcZnFqvbtyu37wcwuBHKB42tYfzlwOUDPnj0bOo6IiIiISI2SYQR9f60BesQtZ8fadmNmJwG/Aca7e3F1B3L3Ke6e6+65HTt2bJCwItK47etd4ERERPZXYy7QZwF9zay3maUB5wLT4zcwsyHAP4kW5xsDyCgiTYTqcxERSZRkKNA37M9O7l4GXAG8CiwCnnD3BWZ2o5mNj212M9AKeNLM5prZ9BoOJyIiIiKSFAKfg+7uJx/Avi8DL1dp+33c65MOIJqISCUNoIuISKIkwwi6iEjS0xx0ERFJFBXoIiIiIiJJJPAC3cwyY3cFFRFJWho/FxGRREl4gW5mITM738xeMrONwGJgnZktNLObzeyQRGcSEdkbzXAREZFECWIE/W3gG8B1QBd37+HunYBjgA+Av8RuLCQikjRcY+giIpIgQVzF5SR3L63a6O5bgKeBp80sNfGxRERERESCl/AR9Iri3Mz+bmZW2zYiIslCU1xERCRRgjxJdDsw3cwyAcxsrJm9F2AeEREREZHABXajInf/rZmdD8wwsxKgEJgcVB4RERERkWQQWIFuZmOAy4AdQFdgkrsvCSqPiEhtNMVFREQSJcgpLr8BfufuJwBnAdPMbHSAeUREaqSruIiISKIEOcVldNzreWZ2KtGruBwdVCYRERERkaAFcaOimq7csg4YU9s2IiJB0RQXERFJlEBuVGRmV5pZz/hGM0sDjjKzqcDEAHKJiNRI9bmIiCRKEFNcxgGTgMfMrDewFcgAwsBrwG3u/kkAuUREauQaQhcRkQRJeIHu7kXAXcBdsTuGdgB2ufvWRGcREREREUk2gV3FxcxeB/q7+zoV5yKS7DR+LiIiiRLkZRZ/BdxmZg+YWdcAc4iI7JVmuIiISKIEVqC7+8fufiLwIvCKmV1vZi2CyiMiIiIikgyCHEGvuJziEuBu4EpgqZldFGQmEZFqaQRdREQSJMg56O8Ba4Bbge7AJcAJwHAzmxJULhGR6uhOoiIikigJv4qLmR0FfABcDiz0Pa9ddqWZLUp0LhGR2mgOuoiIJEoQI+gXA3OA3wETzaxLNducVpcDmdk4M1tiZsvMbHI1648zs4/NrMzMzjqw2CIiIiIiDS+I66D/CMDMDgNOBR40s7bA28ArwHvuvnxvxzGzMHAncDKQB8wys+nuvjBus1VEp878ol7fhIg0OxpAFxGRRAnyKi6L3f1Wdx8HjAbeBb4LfFjHQwwHlrn7cncvAR4HJlTpY4W7fwZE6jG6iDRDupOoiIgkSsJH0Ksys0ygyN1fBl7eh127A6vjlvOAEfWZTUREREQk0RI+gm5mITM738xeMrONRC+zuN7MFprZzWZ2SACZLjez2WY2Oz8/P9Hdi0gjoPFzERFJlCCmuLwNfAO4Duji7tnu3hE4hujVXf5iZhfW4ThrgB5xy9mxtn3m7lPcPdfdczt27Lg/hxCRJk4zXEREJFGCmOJykruXVm109y3A08DTZpZah+PMAvqaWW+ihfm5wPn1mlREJEbXQRcRkURJ+Ah6RXFuZn+P3Um0xm32cpwy4ArgVWAR8IS7LzCzG81sfKyPYWaWR/Tk03+a2YL6eh8iIiIiIg0hyJNEtwPTzexcd99hZmOB37v7qLoeoLoTS93993GvZxGd+iIicmA0gC4iIgkSWIHu7r81s/OBGWZWAhQCe9xsSEQkGag+FxGRRAmsQDezMcBlwA6gKzDJ3ZcElUdEREREJBkEdqMi4DfA79z9BOAsYJqZjQ4wj4hIjXQVFxERSZQgp7iMjns9z8xOJXoVl6ODyiQiUhNdxUVERBIliBsV1XTllnXAmNq2EREJikbQRUQkUQK5UZGZXWlmPeMbzSwNOMrMpgITA8glIiIiIhK4IKa4jAMmAY+ZWR/gK6AF0V8WXgNuc/dPAsglIlIjDaCLiEiiJLxAd/ci4C7grtgdQzsAu9x9a6KziIjUlWuOi4iIJEiQl1m8BrgY2ALMM7NPgU+B+e5eHFQuEREREZEgBXkn0SuBk4FyIAcYDJwBDDCzYncfEGA2EZHdaABdREQSJcgCfQnwhUf/bvwl8HzFCjNrG1gqEREREZEABXmjoo3A/WbWu+oKdy8III+IiIiISOCCLNDnEx3Bf9bM8szsdTP73wDziIjUqLQ8EnQEERFpJoK8k+hfK16bWQrQDxgYVB4RkdoUlapAFxGRxEh4gV71BkVxtgP/jVu/1d23JSiWiEitisvKg44gIiLNRBAj6FOJ3vPDatnGgQeBhxIRSERkbzSCLiIiiRLEjYpOTHSfIiIHqkgj6CIikiBBniQqItJoFGsEXUREEkQFuohIHWwrKq18XaYruoiISANSgS4iUgebCosrX2/ZURJgEhERaepUoIuI1MHarbsqX+fHFesiIiL1TQW6iEgdzMsrIByKXnwqf7sKdBERaTiB3aioPpjZOODvQBi4191vqrI+neilGo8ENgPnuPuKROcUkcbtw+Wb+TSvgAtG9OSRD1exrqAo6Egi0oDcHXeIuBOJPX+9HG3zuHUV673K8u77x20fqfsx49fVeszK/as7/l62j8te2zaVx4zU7ZjxGao+V/+ea3sPNW9f+9e9lv0j0fXdD2rBKz89Luhvu9002gLdzMLAncDJQB4wy8ymu/vCuM2+B3zl7oeY2bnAX4BzEp9WRBqDSMTZUVLG9qIyNm4v5stNhXz05Vc883EevbJa8stxh/Hqgg28u3QT5w2v6Z5r0pCqFilONYUR4JHouvj/pD1un6r/sVc+E7ddNcVXdcd0vi5YaswTKzicfcyz2/qvC6KKPqPHq7JvNX3skTfuPdb8HuL23e097L7v131W/x4iDrBncQTVbFflPUTb9iy22O091L0oq/r1qKlwi3UhQMggZEbIDKt8zdfLIatss7h1Feut1v3jt69+/3DISA1ZjdsbseVQxf5V8u21P6N9ZmrQX+Y9NNoCHRgOLHP35QBm9jgwAYgv0CcAN8RePwX8w8zMvXH/04tEnHJ3yiN7/oCr+h9T/H8A8f9pVPwwjt+/6g/12M9Foq+o/IFV+cyeP8Ti1+2+XLHe415X7lXDttHtq7bFH5/q9qml7/hlan1v1e/jXnvemt571W+5Grf1r7ev6Krys4t7fxXLVPl67v65xeWv5ThfZ9p7P3u0x3/N47PX0E9N31cVbcR9n1aXN76Pr7+2X3/vRiJOWST6byP6HKF8t2WnrDz6XPFvqCwSobTM2VFcRmFJ2R7f0y3TwpyW05VfjTuMti1SmTC4Gw+89yVPzcnj20O6YwZf7Sxl684Stu4qZduuUorLIpTEHtHX5RSXRSiPZSyPQLlHC62KLBXFUXkkWjCUV4xQVayv8m8zvjBx9iycqhZE8Z9N1cKvalFU9XupusIo/mdJ9fvv+fOnukIt/vulLrnkwMUXL7sVOBXLFr2TYEXhZXxdOFUUPJVtITB23zcUd0yzavatXB+/zir7TNmtuIsrrojvI+74WCxr7UXZvheKsbbQ3rentv1366+aQrTWwvLr7Pt0zN2K571vb7VsI8FozAV6d2B13HIeMKKmbdy9zMwKgCxgU0IS1tEbCzdw86tLKItEKIsrICoKjLLy+KIjEvtPTiSxKv7TrviPtOI/xuhKdmurui3xy1XXxfZjt/32PE5lhrh1sUNXHjclZIRDodjz14/UcIiMVNutPSUUij1Hl1tlpNA6PYXWGam0ykghKzONPh0z6dk+k7SUr0/XufqkvsxdvZVfPPkpf3ppISVlEXaW7N9NjMwgbNHiIhxXDIRjy2ZGOETl6+qKofhCKFT5daq+mNp9ufaiqGKb6gqjaEEU+/x3K8jic8UVFrHtoeZcVmXd7setqZjc/X1VZPu6mKmmSGTPgif6Nfh6390Klt2+nl8fs2q+aJEW1xYrjCrec9XPqLIwqylPZe49C6j474HaC+O4Yi3uMxeR5NeYC/R6Y2aXA5cD9OyZ+D9bt8pIoXeHTMLhaLGQUlFghL8uHlJCRko4rvAw+3p0I+4/zd3/I4r95k/1P8Dj/zPZ7Yd/lTagSiH1dUtFgQBxRVSVbb8upL6uqPa6TzXHp5p1u+9bXaY981bdp7a+48Wvr+m97Xn8r/uvMe9ejl+1IK58f1WX2XNbKguk6o9TXb97bKv/0HfTJiOVaZeP5NUFG5j5eT6Z6SlkH9SC9plptG2ZSpuMVNJTQmSkhkgLh0lLCZGeEiI1JfrvN2TRf8ManRIRkZo05gJ9DdAjbjk71lbdNnlmlgK0JXqy6G7cfQowBSA3Nzfh49Mj+2Qxsk9WorsVkf2UEg5xWk5XTsvpGnQUERFpghrzZRZnAX3NrLeZpQHnAtOrbDMdmBh7fRbwVmOffy4iIiIiTVujHUGPzSm/AniV6GUW73f3BWZ2IzDb3acD9wH/MrNlwBaiRbyIiIiISNIyDSjvzszygZUBdN2BJDt5VRqEPufmQZ9z86DPuenTZ9w8BPk5H+zuHas2qkBPEmY2291zg84hDUufc/Ogz7l50Ofc9Okzbh6S8XNuzHPQRURERESaHBXoIiIiIiJJRAV68pgSdABJCH3OzYM+5+ZBn3PTp8+4eUi6z1lz0EVEREREkohG0EVEREREkogK9CRgZuPMbImZLTOzyUHnkfpnZj3M7G0zW2hmC8zs6qAzScMws7CZfWJmLwadRRqGmbUzs6fMbLGZLTKzo4LOJPXPzH4W+3k938weM7OMoDPJgTOz+81so5nNj2trb2avm9nS2PNBQWYEFeiBM7MwcCdwKtAfOM/M+gebShpAGfBzd+8PjAR+os+5yboaWBR0CGlQfwdecffDgEHo825yzKw7cBWQ6+4DiN4QUTc7bBoeBMZVaZsMvOnufYE3Y8uBUoEevOHAMndf7u4lwOPAhIAzST1z93Xu/nHs9Xai/6F3DzaV1DczywZOA+4NOos0DDNrCxxH9E7VuHuJu28NNJQ0lBSghZmlAC2BtQHnkXrg7jOJ3l0+3gRgauz1VODMRGaqjgr04HUHVsct56HCrUkzs17AEODDgKNI/bsN+CUQCTiHNJzeQD7wQGwq071mlhl0KKlf7r4GuAVYBawDCtz9tWBTSQPq7O7rYq/XA52DDAMq0EUSysxaAU8DP3X3bUHnkfpjZqcDG919TtBZpEGlAEOBu919CLCDJPhzuNSv2BzkCUR/IesGZJrZhcGmkkTw6OUNA7/EoQr04K0BesQtZ8fapIkxs1Sixfkj7v5M0Hmk3o0CxpvZCqJT1Uab2cPBRpIGkAfkuXvFX8CeIlqwS9NyEvClu+e7eynwDHB0wJmk4Wwws64AseeNAedRgZ4EZgF9zay3maURPQllesCZpJ6ZmRGds7rI3f8WdB6pf+5+nbtnu3svov+O33J3jbg1Me6+HlhtZv1iTWOAhQFGkoaxChhpZi1jP7/HoJOBm7LpwMTY64nA8wFmAaJ/qpMAuXuZmV0BvEr0LPH73X1BwLGk/o0CLgLmmdncWNuv3f3l4CKJyH66EngkNqiyHLg04DxSz9z9QzN7CviY6FW4PiEJ7zYp+87MHgNOADqYWR5wPXAT8ISZfQ9YCZwdXMIo3UlURERERCSJaIqLiIiIiEgSUYEuIiIiIpJEVKCLiIiIiCQRFegiIiIiIklEBbqIiIiISBJRgS4iIiIikkRUoIuIiIiIJBEV6CIiIiIiSUQFuoiIiIhIElGBLiIiIiKSRFKCDpBsOnTo4L169Qo6hoiIiIg0cXPmzNnk7h2rtqtAr6JXr17Mnj076BgiIiIi0sSZ2crq2jXFRUREREQkiahAFxEREZFm6c1FG/jX+ysoj3jQUXajAl1EREREmp3C4jJ++9x8Hv1oNe7JVaBrDnodlJaWkpeXR1FRUdBRpB5lZGSQnZ1Nampq0FFEREQkwW57/XPWbyvizguGkhJOrjFrFeh1kJeXR+vWrenVqxdmFnQcqQfuzubNm8nLy6N3795BxxEREZEEWrC2gAf+u4LzhvdkaM+Dgo6zh+T6dSFJFRUVkZWVpeK8CTEzsrKy9FcRERGRZiYScX7z7HzatUjlV2MPCzpOtRpVgW5m95vZRjObH9fW3sxeN7OlseeDYu1mZreb2TIz+8zMhh5g3wcaX5KMPlMREZHm57FZq5i7eiu/Pf1w2rZMzmmujapABx4ExlVpmwy86e59gTdjywCnAn1jj8uBuxOUMaFatWp1wMd48MEHueKKK2rdZsWKFTz66KMH3Ne++P73v8/ChQv3aK9LXhEREZGq1hcUcdO/F3NUnyzOHNw96Dg1alQFurvPBLZUaZ4ATI29ngqcGdf+kEd9ALQzs64JCdoEHWiB7u5EIpF92ufee++lf//++92niIiISAV35zfPzqO0PML/fHtgUv8lvVEV6DXo7O7rYq/XA51jr7sDq+O2y4u1NUpnnnkmRx55JEcccQRTpkzZbd3PfvYzjjjiCMaMGUN+fj4At99+O/379ycnJ4dzzz0XgC1btnDmmWeSk5PDyJEj+eyzz/bo55JLLuGpp56qXK4YoZ88eTLvvPMOgwcP5tZbb6W8vJxrr72WYcOGkZOTwz//+c89jrVixQr69evHxRdfzIABA1i9ejU333xz5T7XX389ADt27OC0005j0KBBDBgwgGnTpgFwwgknVN7V9YEHHuDQQw9l+PDhvPfee3vNC1Tbl4iIiDRPz89dy5uLN/KLU/rRq0Nm0HFq1aSu4uLubmb7fCFLM7uc6DQYevbsWeu2f3hhAQvXbtu/gDXo360N159xRK3b3H///bRv355du3YxbNgwvvOd75CVlcWOHTvIzc3l1ltv5cYbb+QPf/gD//jHP7jpppv48ssvSU9PZ+vWrQBcf/31DBkyhOeee4633nqLiy++mLlz59Yp40033cQtt9zCiy++CMCUKVNo27Yts2bNori4mFGjRnHKKafscUWUpUuXMnXqVEaOHMlrr73G0qVL+eijj3B3xo8fz8yZM8nPz6dbt2689NJLABQUFOx2jHXr1nH99dczZ84c2rZty4knnsiQIUNqzVtTX8cdd1yd3q+IiIg0HRu3F3HDCwsY2rMdl45K/qu3NYUR9A0VU1dizxtj7WuAHnHbZcfa9uDuU9w9191zO3bs2KBh99ftt9/OoEGDGDlyJKtXr2bp0qUAhEIhzjnnHAAuvPBC3n33XQBycnK44IILePjhh0lJif4e9u6773LRRRcBMHr0aDZv3sy2bfv3y8Zrr73GQw89xODBgxkxYgSbN2+uzBTv4IMPZuTIkZX7vPbaawwZMoShQ4eyePFili5dysCBA3n99df51a9+xTvvvEPbtm13O8aHH37ICSecQMeOHUlLS6t8v3vLV11fIiIi0vxc//wCdpaU89ezBhEOJe/UlgpNYQR9OjARuCn2/Hxc+xVm9jgwAiiImwqz3/Y20t0QZsyYwRtvvMH7779Py5YtOeGEE2q8PGDFfKqXXnqJmTNn8sILL/DnP/+ZefPm1amvlJSUyrnikUiEkpKSardzd+644w7Gjh1b6/EyM7/+E5K7c9111/GDH/xgj+0+/vhjXn75ZX77298yZswYfv/73x9Q3tr6EhERkebjxc/W8u/56/nVuMM4pNOBX1wjERrVCLqZPQa8D/Qzszwz+x7RwvxkM1sKnBRbBngZWA4sA/4P+HEAketFQUEBBx10EC1btmTx4sV88MEHlesikUjlHOxHH32UY445hkgkwurVqznxxBP5y1/+QkFBAYWFhRx77LE88sgjQLTo79ChA23atNmtr169ejFnzhwApk+fTmlpKQCtW7dm+/btlduNHTuWu+++u3L9559/zo4dO2p9H2PHjuX++++nsLAQgDVr1rBx40bWrl1Ly5YtufDCC7n22mv5+OOPd9tvxIgR/Oc//2Hz5s2Ulpby5JNP7jVvTX2JiIhI87GuYBe/eXY+g7LbctmxyT+1pUKjGkF39/NqWDWmmm0d+EnDJkqMcePGcc8993D44YfTr1+/yikjEB2h/uijj/jTn/5Ep06dmDZtGuXl5Vx44YUUFBTg7lx11VW0a9eOG264gUmTJpGTk0PLli2ZOnXqHn1ddtllTJgwgUGDBjFu3LjKEfCcnBzC4TCDBg3ikksu4eqrr2bFihUMHToUd6djx44899xztb6PU045hUWLFnHUUUcB0RM6H374YZYtW8a1115LKBQiNTWVu+/e/YqYXbt25YYbbuCoo46iXbt2DB48eK95a+qrU6dO+/z1FxERkcYnEnF+8eSnlJRFuPWcwaSEG8+4tEXrWKmQm5vrFVcOqbBo0SIOP/zwgBJJQ9JnKyIi0jT938zl/PnlRfzPtwdy3vDaLwISFDOb4+65Vdsbz68SIiIiIiJ1sHDtNm5+dQkn9+/MucN67H2HJKMCXURERESajKLScq5+/BPatkzlL9/JSeobEtUkIXPQzayuf1fY6u71e5FxEREREWk2/vjiQpZuLGTqpOG0z0wLOs5+SdRJolMBB2r7FcaBB4GHEhFoX7l7o/wNTGqm8y9ERESalufnruGRD1fxg+P7cPyhyXlvm7pISIHu7icmop+GkpGRwebNm8nKylKR3kS4O5s3byYjIyPoKCIiIlIPlm3cznXPzGNYr4P4xSn9go5zQAK5zKKZZQJF7l4eRP/7Kjs7m7y8PPLz84OOIvUoIyOD7OzsoGOIiIjIAdpZUsaPHv6YFqlh7jhvKKmN6JKK1UnUHPQQcC5wATAMKAHSzSwfeAn4p7svS0SW/ZGamkrv3o3n4vYiIiIizYW789tn57Msv5B/TRpBl7aN/6/jifr14m3gG8B1QBd3z3b3jsAxwAfAX8zswgRlEREREZEm4uEPVvLMJ2u4ekxfjunbIeg49SJRU1xOcvfSqo3uvgV4GnjazFITlEVEREREmoD3v9jMH15YyOjDOnHl6L5Bx6k3CRlBryjOzezvVsNZltUV8CIiIiIi1Vm9ZSc/fmQOB2e15LZzBxMONZ0LeSR6Bv12YHrsJFHMbKyZvZfgDCIiIiLSiO0oLuOyh2ZTHnHunTiMNhlNayJGQq/i4u6/NbPzgRlmVgIUApMTmUFEREREGq9IxPn5E5/y+YbtPHjpcHp3yAw6Ur1LaIFuZmOAy4AdQFdgkrsvSWQGEREREWm8bnplMa8sWM9vTzuc4xrxzYhqk+gpLr8BfufuJwBnAdPMbHSCM4iIiIhIIzT1vyuYMnM5Fx91MN87puleAjvRU1xGx72eZ2anEr2Ky9GJzCEiIiIijcurC9ZzwwsLOOnwzlx/xhFN+u7uCRlBr+XKLeuAMbVtIyIiIiLN28ervuKqxz4hJ7sdd5w3pEldsaU6CbtRkZldaWY94xvNLA04ysymAhMTlEVEREREGomlG7bzvQdn0blNBvdNzKVFWjjoSA0uUVNcxgGTgMfMrDewFcgAwsBrwG3u/kmCsoiIiIhII7Bq804uvO9DUsIhHpo0nA6t0oOOlBAJKdDdvQi4C7grdsfQDsAud9+aiP5FREREpHFZX1DE+fd+QHFZhGmXH0WvJng5xZok9CouZvY60N/d16k4FxEREZHqbCos5oJ7P2DrzlIemjScfl1aBx0poRJ9mcVfAbeZ2QNm1jXBfYuIiIhIkvtqRwkX3/cRa7bu4v5LhpGT3S7oSAmX0ALd3T929xOBF4FXzOx6M2uRyAwiIiIikpw2FxZz3v99wLL8Qv55US7De7cPOlIgEj2CXnE5xSXA3cCVwFIzuyjROUREREQkeeRvjxbnX27awX0Tczm+id4ltC4SPQf9PWANcCvQHbgEOAEYbmZTEplFRERERJLDhm1FnDvlfVZv2cUDlw7j2L7NtziHBN9JFLgcWOjuXqX9SjNblOAsIiIiIhKw1Vt2ctF9H5K/vZipk4Y322kt8RJaoLv7glpWn5awICIiIiISuEXrtjHx/o8oKi3noe+N4MiDDwo6UlJI9Ah6jdx9edAZRERERCQxPli+mcumzqZVRgpP/ehoDu3cvC6lWJuEnyQKYGZnBNGviIiIiATvlfnruPj+j+jcNoOnVZzvIZACHfhzQP2KiIiISEDcnfve/ZIfPfIxA7q14akfHkW3drridlVBTXGxgPoVERERkQCUlEW4fvp8HvtoNWOP6Mxt5wyhRVo46FhJKagCvepVXA6Yma0AtgPlQJm755pZe2Aa0AtYAZzt7l/Vd98iIiIiUrOvdpTww4fn8OGXW7jixEO45uRDCYU0XluToKa4NJQT3X2wu+fGlicDb7p7X+DN2LKIiIiIJMjSDduZcOd7fLJ6K7edM5hfjO2n4nwvmlqBXtUEYGrs9VTgzOCiiIiIiDQvz89dw4Q732NnSTmPXz6SM4d0DzpSoxDUFJcNDXBMB14zMwf+6e5TgM7uvi62fj3QuQH6FREREZE4xWXl/OnFRfzrg5XkHnwQ/zh/KF3aZgQdq9EIpEB395Mb4LDHuPsaM+sEvG5mi6v06bHifQ9mdjnRu5zSs2fPBogmIiIi0jzkfbWTnzzyMZ/mFXDZsb355bjDSA039Ukb9StpblR0oNx9Tex5o5k9CwwHNphZV3dfZ2ZdgY017DsFmAKQm5tb7yewioiIiDQH/563jsnPzCMSce65cCjjBnQNOlKj1CR+nTGzTDNrXfEaOAWYD0wHJsY2mwg8H0xCERERkaarsLiMXzz5KT965GMOzmrJ9CuPUXF+AAIZQY8V0UXuXl5Ph+wMPGtmEH1Pj7r7K2Y2C3jCzL4HrATOrqf+RERERASYs/IrfjZtLnlf7eSKEw/h6pP6akrLAUpIgW5mIeBc4AJgGFAMpJvZJuAloid1Ltvf47v7cmBQNe2bgTH7e1wRERERqV5RaTl3vLWUe/6znK5tM5j2g6MY1qt90LGahESNoL8NvAFcB8x39whA7EZCJwJ/MbNn3f3hBOURERERkf00a8UWfvX0ZyzP38F3j8zm92f0p3VGatCxmoxEFegnuXtp1UZ33wI8DTxtZvpURURERJLY9qJS/vrKEv71wUqyD2rBQ5OGc9yhHYOO1eQkpECvKM7N7O/AT919jyulVFfAi4iIiEjw3J1X5q/njy8uZN22IiaN6s3PTzmUzPQmc0HApJLor+p2YLqZnevuO8xsLPB7dx+V4BwiIiIiUgdLN2znhhcW8N6yzRzWpTX/uGAoQ3seFHSsJi2hBbq7/9bMzgdmmFkJUAhMTmQGEREREdm7bUWl3Pb6Uqa+v4JW6SncOOEIzh/ekxRdoaXBJbRAN7MxwGXADqArMMndlyQyg4iIiIjUrKQswuOzVnH7m0vZvKOEc4f15Nqx/WifmRZ0tGYj0VNcfgP8zt3fNbOBwDQzu8bd30pwDhERERGJE4k4L85bxy2vLmHVlp2M6N2eBy7pz8DstkFHa3YSPcVldNzreWZ2KtGruBydyBwiIiIiEuXuvLN0E399dTHz12zjsC6teeDSYZxwaEdiN4GUBEvUjYqshiu3rItNe6lxGxERERGpf+7OjCX53P7WUj5ZtZXu7Vpw6zmDmDCoO6GQCvMgJexGRWb2NPC8u6+qaDSzNOAoM5tI9GZGDyYoj4iIiEizFIk4ry/awD/eWsa8NQV0b9eCP505gO/mZpOeEg46npC4An0cMAl4zMz6AF8BLYAQ8Bpwm7t/kqAsIiIiIs1OcVk50+eu5b53v2Tx+u0cnNWSv34nh28N7U6qrsySVBJ1o6Ii4C7grtgdQzsAu9x9ayL6FxEREWmuNhcW8/AHq/jXByvZVFhMv86t+dvZgxg/qJsumZikEn2ZxWuAi4EtwDwz+xT4FJjv7sWJzCIiIiLSlM1fU8C/3l/Js3PXUFIW4cR+HfneMX0YdUiWTv5Mcom+zOKVwMlAOZADDAbOAAaYWbG7D0hwHhEREZEmo7C4jOlz1/LYR6uYt6aAjNQQZx2ZzaRRvTmkU6ug40kdJbpAXwJ8Ebtay5fA8xUrzEwX2RQRERHZR+7OZ3kFPD5rFc/PXcvOknIO69KaGyccwYTB3WnbIjXoiLKPEl2gbwTuN7Mb3f3L+BXuXpDgLCIiIiKN1srNO3juk7U8P3cNyzftoEVqmDMGdeW84T0Z3KOdprE0Yoku0OcDA4FnzawDsAj4zN1/nuAcIiIiIo3OpsJiXvpsHc9+soa5q7diBiN6t+fy4/rwzZyutMnQaHlTkOg7if614rWZpQD9iBbsIiIiIlKN1Vt28uqC9by2YAOzV24h4nBYl9Zcd+phnDGoG93atQg6otSzRN1JtGcNq7YD/41bv9XdtyUik4iIiEgycneWbNjOaws28Mr89SxcFy2NDuvSmitG9+W0gV3p16V1wCmlISVqBH0q4EBtk6Gc6J1EH0pEIBEREZFkUbCrlPeWbeI/S/KZuTSfdQVFmMHQngfx628extgjunBwVmbQMSVBEnWjohMT0Y+IiIhIY1BaHuGzvIJoUf55PnNXb6U84rTOSOGYQzpw1ZiOjDm8E51aZwQdVQKQ6JNERURERJqd4rJy5q7aykdfbuHDL7cwZ+VX7Cotxwxyurflxyd8g+MP7cjgHu10d09RgS4iIiJS3zZuL2Luqq18mreV2Su+4pPVWykpi2AG/Tq35pxhPRjRuz0j+mTRPjMt6LiSZFSgi4iIiByAnSVlzMsr4NO8rcxdvZVPVxewZusuAFJCxuFd23DxyIMZ0SeLYb0Ool1LFeRSOxXoIiIiInXg7qwrKGLRum2xx3YWrdvGl5t34B7dpkf7Fgzp2Y5LR/ViSM92HNGtLRmp4WCDS6OjAl1EREQkjruTv72YZfmFfJG/gy82FrJo3TYWr99Owa7Syu16tG/B4V3acMagbuRkt2VQj3Z0aJUeYHJpKlSgi4iISLO0s6SM1Vt28eWmrwvxL/ILWZ6/g+3FZZXbtUwLc2jn1nxzYFf6d23N4V3b0K9La1rrrp3SQFSgi4iISJNUWh5h7dZdrN6yi9Vf7WTVlp2s3rKT1V/tIm/LTjbvKNlt+y5tMvhGp0y+NbQ73+jYim90bEWfjpl0aZNBKFTbrVxE6pcKdBEREWl0dpaUsb6giPXbitiwrYh1BUVsiC2v31bM+oJd5G8vJuJf75MSMrq1a0GP9i04uX9nerRvSfZBLejToRW9O2bSKl1lkSQHfSeKiIhI4NydbbvK2LSjmM2FJWwuLGbTjhK2FJawOda2qbCYzTtK2LitiG1FZXsco3VGCl3aZNClbQZ9O3Wka9sMehzUkuz2LejZviVd2mToGuPSKKhAFxERkXpTXFZOwa5Stu0qpWBXKVt3Rp+rPrbFvd66s5QtO0ooix/ujtOuZSpZmWlkZabTt1Mrjv5GFl3aZkSL8TYZdI69ztQIuDQRTf472czGAX8HwsC97n5TwJFEREQCV1Yeoagswq6ScopKyykuK6eoNEJRaTk7S8rZUVxGYXEZO4rL2BFbjraVs7Mkbl1xOTtKvn5dUh6ptd9W6Sm0bZFKmxaptG2RQu8OmbRrkUZWqzSyWqXToVUa7WPFeIdWaRyUmUaqRr2lmWnSBbqZhYE7gZOBPGCWmU1394XBJhMRkabO3SmPOOUVz3GPsohTUhahtDxCablTWh6hpDxCaVmV5YpHme++XB6/f2y5PEJxaYSisnKKS78utneVRgvwotLIbkV4TaPVNQmHjMy0MK3SU8iMPVqlp9ChVfpuba0zUmLF956PNhkpmmIiUgdNukAHhgPL3H05gJk9DkwAkqpAX7t1F/PWFOx1O6/zz9K6bVjX49W127ofL6h8+/af0d6PV8ft6vv9Bvh1qedvwQC/F+p6vGDy1fWA9f5+G6DfSOyYEXciXqUtEmujYp1XrvfYPh5rj/iey1/vV9FWsX/sWLv1G9sXJxL5uo+KY1UtoMs9WkRHIvHPESJO9DkSfS6PQHkkstt+Fa/3sf7dL6lhIzUcqnykhY2MtDAZKWEyUkNkpIZp0yI1+jolvMe6yueUMOmpIVqkhslIDdMyLUzLtJRY4R0mMz2F9JQQZrqSiUgiNPUCvTuwOm45DxhRdSMzuxy4HKBnz56JSRbng+WbueaJTxPer4hIMjGDkBlG7Nm+bqtoN4NQqOo2RsjAiD3H2kO2+3LFPnvsZxAOhQgbpIRChEKQFgoTChnhinWhinVGSih6jJSQVS6Hqz5sz7bK/cLR5WhBXVFcG6kpIdLDIVJTvm6rXJ9SZTm2XgWzSNPU1Av0OnH3KcAUgNzc3ASMeexu9GGdeOmqY+q0bfS/pTpsV8ef2XXerr77rdtmdT5eXY9Y//nq2G+dj1fH7er5/e6Lun9vNY2vTT1/Cwb4bymgzyOuSA7FimbYvXiueBYRkaimXqCvAXrELWfH2pJKu5ZptGuZFnQMEREREUkCTf1MjVlAXzPrbWZpwLnA9IAziYiIiIjUqEmPoLt7mZldAbxK9DKL97v7goBjiYiIiIjUyOr7yhaNnZnlAysD6LoDsCmAfiWx9Dk3D/qcmwd9zk2fPuPmIcjP+WB371i1UQV6kjCz2e6eG3QOaVj6nJsHfc7Ngz7npk+fcfOQjJ9zU5+DLiIiIiLSqKhAFxERERFJIirQk8eUoANIQuhzbh70OTcP+pybPn3GzUPSfc6agy4iIiIikkQ0gi4iIiIikkRUoCcBMxtnZkvMbJmZTQ46j9Q/M+thZm+b2UIzW2BmVwedSRqGmYXN7BMzezHoLNIwzKydmT1lZovNbJGZHRV0Jql/Zvaz2M/r+Wb2mJllBJ1JDpyZ3W9mG81sflxbezN73cyWxp4PCjIjqEAPnJmFgTuBU4H+wHlm1j/YVNIAyoCfu3t/YCTwE33OTdbVwKKgQ0iD+jvwirsfBgxCn3eTY2bdgauAXHcfQPRmh+cGm0rqyYPAuCptk4E33b0v8GZsOVAq0IM3HFjm7svdvQR4HJgQcCapZ+6+zt0/jr3eTvQ/9O7BppL6ZmbZwGnAvUFnkYZhZm2B44D7ANy9xN23BhpKGkoK0MLMUoCWwNqA80g9cPeZwJYqzROAqbHXU4EzE5mpOirQg9cdWB23nIcKtybNzHoBQ4APA44i9e824JdAJOAc0nB6A/nAA7GpTPeaWWbQoaR+ufsa4BZgFbAOKHD314JNJQ2os7uvi71eD3QOMgyoQBdJKDNrBTwN/NTdtwWdR+qPmZ0ObHT3OUFnkQaVAgwF7nb3IcAOkuDP4VK/YnOQJxD9hawbkGlmFwabShLBo5c3DPwShyrQg7cG6BG3nB1rkybGzFKJFuePuPszQeeRejcKGG9mK4hOVRttZg8HG0kaQB6Q5+4VfwF7imjBLk3LScCX7p7v7qXAM8DRAWeShrPBzLoCxJ43BpxHBXoSmAX0NbPeZpZG9CSU6QFnknpmZkZ0zuoid/9b0Hmk/rn7de6e7e69iP47fsvdNeLWxLj7emC1mfWLNY0BFgYYSRrGKmCkmbWM/fweg04GbsqmAxNjrycCzweYBYj+qU4C5O5lZnYF8CrRs8Tvd/cFAceS+jcKuAiYZ2ZzY22/dveXg4skIvvpSuCR2KDKcuDSgPNIPXP3D83sKeBjolfh+oQkvNuk7Dszeww4AehgZnnA9cBNwBNm9j1gJXB2cAmjdCdREREREZEkoikuIiIiIiJJRAW6iIiIiEgSUYEuIiIiIpJEVKCLiIiIiCQRFegiIiIiIklEBbqIiIiISBJRgS4iIiIikkRUoIuIiIiIJJH/D62bgWYG68rYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    For the given initial value problem y' = (x + 1)^2 + 1.0, taking x from 0 to 10 in 3000 steps,\n",
      "    we find a large deviation between exact and estimate solutions at the point that y' = 0.\n",
      "    We could fix this by taking some other measure, like the absolute difference, but as\n",
      "    the solutions decrease so hard, that will suck later on.\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "# This here is for debugging. We thought the keep it here,\n",
    "# as debugging was part of the task.\n",
    "\n",
    "def ODEF(x, y):\n",
    "    return y - x**2 + 1.0\n",
    "\n",
    "x = np.linspace(0, 10, 3000)\n",
    "y0 = .5\n",
    "\n",
    "eta = integrator(x, y0, ODEF, phi_euler)\n",
    "y_exact = (x + 1)**2 - 0.5*np.exp(x)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "gs = gridspec.GridSpec(3, 1, height_ratios=[2, 1, 1])\n",
    "\n",
    "# Above plots\n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.plot(x, eta, label=\"Euler method solution\")\n",
    "ax0.plot(x, y_exact, label=\"exact solution\")\n",
    "ax0.legend()\n",
    "plt.title(\"Solutions to $y' = y - x^2 + 1.0$\")\n",
    "plt.xlabel(\"x\")\n",
    "ax0.set_ylabel(\"y(x)\")\n",
    "\n",
    "# Relative absolute residue plots\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.plot(x, np.abs(eta/y_exact - 1), label=\"relative absolute residue\")\n",
    "ax1.legend()\n",
    "ax1.set_ylabel(\"$|\\eta(x)/y(x) - 1|$\")\n",
    "\n",
    "# Absolute residue plots\n",
    "ax2 = plt.subplot(gs[2])\n",
    "ax2.plot(x, np.abs(eta - y_exact), label=\"absolute residue\")\n",
    "ax2.legend()\n",
    "ax2.set_ylabel(\"$|\\eta(x) - y(x)|$\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"\"\"\n",
    "    For the given initial value problem y' = (x + 1)^2 + 1.0, taking x from 0 to 10 in 3000 steps,\n",
    "    we find a large deviation between exact and estimate solutions at the point that y' = 0.\n",
    "    We could fix this by taking some other measure, like the absolute difference, but as\n",
    "    the solutions decrease so hard, that will suck later on.\n",
    "    \"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "aea0c6713e9bd3871c2e93e65e93d903",
     "grade": true,
     "grade_id": "cell-715c304f22415d45",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def test_integrator():\n",
    "    def ODEF(x, y):\n",
    "        return y - x**2 + 1.0\n",
    "    \n",
    "    x = np.linspace(0, 3, 100000)\n",
    "    y0 = .5\n",
    "\n",
    "    eta = integrator(x, y0, ODEF, phi_euler)\n",
    "    y_exact = (x + 1)**2 - 0.5*np.exp(x)\n",
    "    \n",
    "    # Assert the relative error of the approximation is lower than TOL.\n",
    "    TOL = 1e-4\n",
    "    assert np.all(np.abs(eta/y_exact - 1) < TOL)\n",
    "    \n",
    "test_integrator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "8048460a3a4d5e94b9e9f086d90ff0d0",
     "grade": false,
     "grade_id": "cell-99399435c164c96d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 3\n",
    "\n",
    "Extend the set of possible single-step integration schemes to the modified Euler (Collatz), Heun, and 4th-order Runge-Kutta approaches by implementing the functions $\\text{phi_euler_modified(x, y, f, i)}$, $\\text{phi_heun(x, y, f, i)}$, and $\\text{phi_rk4(x, y, f, i)}$. These can be used in your $\\text{integrator}$ function instead of $\\text{phi_euler}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "dc82ac823cbd262952a6990776e26d09",
     "grade": true,
     "grade_id": "cell-ea865d2efbc574d6",
     "locked": false,
     "points": 9,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_euler_modified(x, y, f, i):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    \n",
    "def phi_heun(x, y, f, i):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    \n",
    "def phi_rk4(x, y, f, i):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "996fa4ae58be6c15aed465d8d7e5618e",
     "grade": false,
     "grade_id": "cell-a783908a0db32a24",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 4\n",
    "\n",
    "Add the possibility to also handle the following multi-step integration schemes:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\Phi_{AB3}(x, y, f, i) = \n",
    "    \\frac{1}{12} \\left[ \n",
    "        23 f(x_i, y_i) \n",
    "        - 16 f(x_{i-1}, y_{i-1})\n",
    "        + 5 f(x_{i-2}, y_{i-2})\n",
    "    \\right]\n",
    "\\end{align*}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{align*}\n",
    "\\Phi_{AB4}(x, y, f, i) = \n",
    "    \\frac{1}{24} \\left[ \n",
    "        55 f(x_i, y_i) \n",
    "        - 59 f(x_{i-1}, y_{i-1})\n",
    "        + 37 f(x_{i-2}, y_{i-2})\n",
    "        - 9 f(x_{i-3}, y_{i-3})\n",
    "    \\right]\n",
    "\\end{align*}\n",
    "        \n",
    "In these cases the initial condition $\\text{y0}$ must be an array consisting of several initial values corresponding to $y_0$, $y_1$, $y_2$ (and $y_3$). Use the Runga-Kutta method to calculate all of the initial values before you start the AB3 and AB4 integrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "d301f6514b2bc50a77d929518f6928df",
     "grade": true,
     "grade_id": "cell-0bd05d93759bd652",
     "locked": false,
     "points": 6,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def phi_ab3(x, y, f, i):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "def phi_ab4(x, y, f, i):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "c621385b0e71eb6fd2e4d6041268bd15",
     "grade": false,
     "grade_id": "cell-3fac52fc46f0a497",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 5\n",
    "\n",
    "Plot the absolute errors $\\delta(x) = |\\tilde{y}(x) - y(x)|$ with a $y$-log scale for all approaches with $0 \\leq x \\leq 2$ and a step size of $h=0.02$ for the ODE from task 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "22f375f78eeb46befb192c08db0f48e9",
     "grade": true,
     "grade_id": "cell-e61d969c23914a5f",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "9fd7c159424b1d536c11d53764cce54c",
     "grade": false,
     "grade_id": "cell-3c99ac2c49457d60",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 6\n",
    "\n",
    "Study the accuracies of all approaches as a function of the step size $h$. To this end use your implementations to solve\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = y\n",
    "\\end{align*}\n",
    "        \n",
    "with $y(x=0) = 1.0$ for $0 \\leq x \\leq 2$. \n",
    "\n",
    "Plot $\\delta(2) = |\\tilde{y}(2) - y(2)|$ as a function of $h$ for each integration scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0060c86da7630a17a6da4769bb0c728f",
     "grade": true,
     "grade_id": "cell-552ef9b066bc9b1c",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "0e3fcb92129c0b3ef4e2b8f99a38ddef",
     "grade": false,
     "grade_id": "cell-3070782f40a2599f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 7\n",
    "\n",
    "Apply the Euler and the Runga-Kutta methods to solve the pendulum problem given by\n",
    "\n",
    "\\begin{align*}\n",
    "    \\vec{y}'(x) = \\vec{f}(x, \\vec{y})\n",
    "    \\quad \\leftrightarrow \\quad \n",
    "    \\begin{pmatrix}\n",
    "    y_0'(x) \\\\\n",
    "    y_1'(x)\n",
    "    \\end{pmatrix}\n",
    "    =\n",
    "    \\begin{pmatrix}\n",
    "        y_1(x) \\\\\n",
    "        -\\operatorname{sin}\\left[y_0(x)  \\right]\n",
    "    \\end{pmatrix}\n",
    "\\end{align*}\n",
    "        \n",
    "To this end the quantities $\\text{x}$ and $\\text{y0}$ in your $\\text{integrator}$ function must become vectors. Depending on your implementation in task 1 you might have to define a new $\\text{integrator}$ function that can handle this type of input.\n",
    "\n",
    "Plot $y_0(x)$ for several oscillation periods. Does the Euler method behave physically?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "9cb9ca1fb8e9c0722a7e6b5ec0ff0c2a",
     "grade": true,
     "grade_id": "cell-3a89944b55976666",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}