notebook/Solving ODE with SciPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a836bea1",
   "metadata": {},
   "source": [
    "# Solving ODEs with SciPy\n",
    "Gilbert François Duivesteijn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f0fb4980",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "MathJax.Hub.Config({\n",
       "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
       "});\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14241fc4",
   "metadata": {},
   "source": [
    "## Example: Radioactive decay\n",
    "\n",
    "Radioactive decay is described as:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dC}{dt} = -\\frac{C}{\\tau}\n",
    "\\end{align}\n",
    "\n",
    "The exact solution is:\n",
    "\\begin{align}\n",
    "C(t) = C_0 e^{\\frac{-t}{\\tau}},\n",
    "\\end{align}\n",
    "where $C_0$ is the initial condition. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "280b41db",
   "metadata": {},
   "source": [
    "Solve this problem numerically with $\\tau = 1$, $C_0 = 1$ and $t_{end} = 5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "89f0f96e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bcde35dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# times\n",
    "t = np.linspace(0, 5, 200)\n",
    "\n",
    "# initial condition\n",
    "C0 = 1.0\n",
    "\n",
    "# Constant\n",
    "tau = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "70688586",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Right hand side of function (1)\n",
    "\n",
    "def fn(C, t, tau):\n",
    "    return -C / tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2c1fca4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical solution\n",
    "\n",
    "C_num = odeint(fn, C0, t, args=(tau,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d52d4f5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytical solution\n",
    "\n",
    "C_exact = C0 * np.exp(-t/tau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c1056fef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAF1CAYAAADyXFc5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABLt0lEQVR4nO3dd3iUZd7F8fMjCSVUF5Aa6ShVpAuKSC+CBURwwS66i65rX110Vayoq6voImLvioUiSA9FQelKXZoKSO8QWsj9/jET3gABBlLuKd/PdXGRTCaTk/F993BmnnnGnHMCAAAAACC35PEdAAAAAAAQWxiiAAAAAIBcxRAFAAAAAOQqhigAAAAAIFcxRAEAAAAAuYohCgAAAADIVQxRIIeZ2cNmNjSLt1HRzJyZxWdjrpZmtjYL35/l3wsAgGhlZjeY2fQsfP9gM3skOzMB4STb/lELIHPOuad9Z8gqM2sp6UPnXPn0y6Lh9wIAwMzelbTWOdffY4YbJN3inLso/TLn3O2+8gC5gWdEgRyUnc9gAgAAANGCIYqoYGa/mtl9Zvazme00s8/MLH/wa8cdGhM8zLVq8ON3zex1MxtjZnvM7HszK21mL5vZdjNbamYXZPjesmb2pZltNrPVZva3DF97zMyGmdmHZrZL0g3Byz7McJ2LzOwHM9thZmuCj4LKzDqb2Twz2xW8/LHT+P0fNLN1ZrbbzJaZWevg5fmCv8cfwT8vm1m+E9zGkfskw/3ypJkVlDRGUtng/bMneB8c+3t1NbNFwd8r2cxqhPLfBwCAjE7Us2b2JzNba2Zdgp8XMrMVZnZd8POT9mhm/WtmfSX9WdIDwX4bmUkeM7OXzGxT8LZ/MbPawa8VNbP3g1l/M7P+Znbcv68tk5fYBLvylmBfDpZ0YTDDjuDX3zWzJzNc/9bg77vNzEaYWdkMX3NmdruZLQ/+fq+ZmZ35fwUg5zFEEU16SOogqZKkupJuOM3v7S+phKQDkmZImhv8fJikf0tSsFxGSlogqZyk1pL+bmbtM9zW5cHvKSbpo4w/xMwqKDDqXpVUUlI9SfODX94r6brg93WW9Bczu+JUwc3sXEl3SGrknCssqb2kX4Nf/qekpsGfc76kxsHfM2TOub2SOkr6wzlXKPjnj2MyVJf0iaS/B3+v0ZJGmlneDFfLyn8fAEAMOFnPOue2SbpJ0ptmdraklyTNd869H/z2E/boifrXOTdEga4eGOy3LpnEaiephaTqkooq0Gdbg197NXhZZUmXBH/+jafzOzvnlki6XdKMYIZimdwvrSQ9E/zZZST9JunTY652maRGCnRsDwX+PQCELYYooskrzrk/gkU1UoGSCdXXzrk5zrn9kr6WtN85975z7rCkzySlPyPaSFJJ59wTzrmDzrlVkt6U1DPDbc1wzn3jnEtzzu075udcK2mCc+4T59wh59xW59x8SXLOJTvnfgl+388KDLtLQsh+WFI+STXNLME596tzbmXwa3+W9IRzbpNzbrOkxyX1OY37JVTXSPrWOTfeOXdI0guSCkhqluE6WfnvAwCIDSftWefcOElfSJooqZOk29K/8RQ9esL+DcEhSYUlnSfJnHNLnHPrzSwumOsh59xu59yvkl5UzvTsnyW97Zyb65w7IOkhBZ5BrZjhOs8653Y4536XNFn0LMIcQxTRZEOGj1MkFTqN792Y4eN9mXyeflsVFDhEdUf6H0kPSyqV4fprTvJzkiStzOwLZtbEzCYHD+/ZqcCjoyVOFdw5t0KBZyIfk7TJzD7NcLhOWQUeNU33W/Cy7HbUz3HOpSlwP5TLcJ2s/PcBAMSGUHp2iKTakt51zqU/M3mqHj1h/56Kc26SpEGSXlOgZ4eYWZHgbSfo+J4td/ytZNmxPbtHgWdl6VlELIYoYsFeSYnpn5hZ6Szc1hpJq51zxTL8Keyc65ThOu4U31/lBF/7WNIISUnOuaIKvF4kpNd3OOc+Dp5pr0Lw5z8X/NIfwcvSnRO8LDMpynA/Scp4P53sdzru5wRfl5Ikad0pwwMA8P9O2rPBZyGHSHpf0l8znttAJ+/Rk/XvqTpOzrlXnHMNJNVU4BDd+yVtUeDZ0mN7NrPu2xv8O7t6tqCk4if4WUBEYIgiFiyQVMvM6gVPkPNYFm7rJ0m7LXByoAJmFmdmtc2sUYjf/5GkNmbWw8zizay4mdULfq2wpG3Ouf1m1liBw4hOyczONbNWFjgJ0X4FnsFNC375E0n9zaykmZWQ9KikD09wU/MlXRv8nTro6MOCN0oqbmZFT/C9n0vqbGatzSxB0r0KvNb2h1B+BwAAgk7Vsw8rMNpukvS8pPeD41Q6eY+erH83KvAaz0yZWaPgs60JCgzK/ZLSgi/f+VzSU2ZWOPg61HuUSc8GXx6zTlLv4O90k44exhsllT/m3AoZfSLpxuC/ZfJJelrSj8HDgYGIxBBF1HPO/U/SE5ImSFou6YzfXDpYOpcp8LqL1Qo8GjpUgRMVhPL9vyvwmpZ7JW1TYPydH/zyXyU9YWa7FRiMn4cYK5+kZ4NZNkg6W4HXjkjSk5JmS/pZ0i8KnIDpyUxuQ5LuktRF0g4FXovyTYbcSxUowVXBQ6WOOrzXObdMUm8FTtqwJXg7XZxzB0P8HQAAOGnPmlkDBYbedcHrPafAKP1H8NtP2KOn6N+3FDjPwg4z+yaTWEUUeJ3qdgUOj92qwAiWpDsVGKerFPj3xceS3j7Br3erAs+kbpVUS0c/WDtJ0iJJG8xsSyb3ywRJj0j6UtJ6BUZsz2OvB0QSc+6URyMAAAAAAJBteEYUAAAAAJCrGKIAAAAAgFzFEAUAAAAA5CqGKAAAAAAgVzFEAQAAAAC5Kt7XDy5RooSrWLGirx8PAIgyc+bM2eKcK+k7RySjmwEA2elk3extiFasWFGzZ8/29eMBAFHGzH7znSHS0c0AgOx0sm7m0FwAAAAAQK5iiAIAAAAAchVDFAAAAACQqxiiAAAAAIBcxRAFAAAAAOQqhigAAAAAIFcxRAEAAAAAuYohCgAAAADIVQxRAAAAAECuOuUQNbO3zWyTmS08wdfNzF4xsxVm9rOZ1c/+mAAAIB3dDACIdKE8I/qupA4n+XpHSdWCf/pK+m/WYwEAgJN4V3QzACCCnXKIOuemStp2kqtcLul9FzBTUjEzK5NdAU9m9+7dGjFihHbt2pUbPw4AgLAQzt28a9cuDR8+XHv27MmNHwcAiFDZ8RrRcpLWZPh8bfCy45hZXzObbWazN2/enOUfPHfuXF1++eWaMmVKlm8LAIAo4q2bZ82apSuuuELTpk3L8m0BAKJXrp6syDk3xDnX0DnXsGTJklm+vSZNmihfvnyaPHlyNqQDACD2ZHc3X3jhhUpISKCbAQAnlR1DdJ2kpAyflw9eluPy58+vZs2aUXYAABzNWzcnJiaqadOmdDMA4KSyY4iOkHRd8Ax9TSXtdM6tz4bbDcmll16qBQsWaNu2k71UBgCAmOK9m+fOnaudO3fm1o8EAESYUN6+5RNJMySda2ZrzexmM7vdzG4PXmW0pFWSVkh6U9JfcyxtJi699FI553idKAAgZkRCN6elpWnq1Km5+WMBABEk/lRXcM71OsXXnaR+2ZboNDVu3FiJiYmaPHmyrrzySl8xAADINeHezU2bNj1yDocuXbr4igEACGO5erKinJA3b141b96c16IAABAmOIcDAOBUIn6ISoFDgBYuXKjsOO08AADIupYtW3IOBwDACUXNEJXE60QBAAgT6edw4HWiAIDMRMUQbdCggQoVKsQhQAAAhInGjRurQIECdDMAIFNRMUQTEhJ08cUXU3YAAISJfPnycQ4HAMAJRcUQlQKvRVmyZIk2bNjgOwoAAFDg8NxffvmFczgAAI4TNUM0/XWiycnJfoMAAABJnMMBAHBiUTNEL7jgAhUpUoRDgAAACBMNGzZUwYIF6WYAwHGiZojGx8erRYsWlB0AAGEiISFBF110Ed0MADhO1AxRSWrVqpWWL1+uNWvW+I4CAAAU6OYlS5Zo/fr1vqMAAMJIVA3RNm3aSJImTpzoOQkAAJDoZgBA5qJqiNauXVulSpXS+PHjfUcBAACS6tWrp+LFi2vChAm+owAAwkhUDVEzU5s2bTRhwgQ553zHAQAg5uXJk0etW7fW+PHj6WYAwBFRNUSlwCFAmzZt0sKFC31HAQAACnTzH3/8oaVLl/qOAgAIE1E3RFu3bi1JHAIEAECYSH+dKN0MAEgXdUM0KSlJ5557Lq8TBQAgTFSqVEmVK1dmiAIAjoi6ISpJbdu21ZQpU3Tw4EHfUQAAgALdPHnyZKWmpvqOAgAIA1E5RNu0aaOUlBTNnDnTdxQAAKBAN+/evVuzZs3yHQUAEAaicoi2bNlSefLk4RAgAADCxKWXXioz46UzAABJUTpEixYtqsaNG1N2AACEieLFi6t+/fo8SAwAkBSlQ1QKvBblp59+0s6dO31HAQAACnTzjBkztGfPHt9RAACeRe0QbdOmjdLS0pScnOw7CgAAUKCbU1NTNXXqVN9RAACeRe0Qbdq0qRITEzk8FwCAMNG8eXPlz5+fbgYARO8QzZs3r1q2bKlx48b5jgIAACTlz59fF198Md0MAIjeISpJHTp00PLly7Vq1SrfUQAAgALdvHjxYq1Zs8Z3FACAR1E/RCVp7NixnpMAAACJbgYABET1EK1ataoqV66s7777zncUAAAgqUaNGkpKSqKbASDGRfUQNTO1b99eEydO1MGDB33HAQAg5qV38/jx43Xo0CHfcQAAnkT1EJUChwDt3btX33//ve8oAABAgW7etWuXfvzxR99RAACeRP0QvfTSS5WQkMAhQAAAhInWrVsrLi6ObgaAGBb1Q7Rw4cK66KKLOCkCAABholixYrrwwgvpZgCIYVE/RCWpffv2WrBggf744w/fUQAAgALdPHv2bG3atMl3FACABzExRNNPFc8baAMAEB7Su3n8+PGekwAAfIiJIVq3bl2VLl2a16IAABAm6tevrxIlStDNABCjYmKImpk6dOig8ePH6/Dhw77jAAAQ8/LkyaP27dtr3LhxSktL8x0HAJDLYmKISoHXomzbtk2zZ8/2HQUAACjQzZs2bdL8+fN9RwEA5LKYGaJt27ZVnjx5NHr0aN9RAACApHbt2kkS3QwAMShmhmjx4sXVtGlTffvtt76jAAAASaVKlVLjxo3pZgCIQTEzRCXpsssu05w5c7R+/XrfUQAAgKTOnTvrxx9/1ObNm31HAQDkopgaop07d5bEIUAAAISLzp07yzmnMWPG+I4CAMhFMTVE69Spo6SkJI0aNcp3FAAAIOmCCy5QmTJl6GYAiDExNUTNTJ07d9b48eN14MAB33EAAIh5efLkUadOnTR27FgdOnTIdxwAQC6JqSEqBQ4B2rt3r6ZOneo7CgAAUOAcDrt27dL333/vOwoAIJfE3BBt1aqV8ufPzyFAAACEiTZt2ihv3rx0MwDEkJgboomJiWrVqpVGjRol55zvOAAAxLxChQqpZcuWvI0LAMSQmBuiUuDw3FWrVmnZsmW+owAAAAW6eenSpVq5cqXvKACAXBCzQ1QSj7wCABAm6GYAiC0xOUQrVKig2rVrU3YAAISJKlWq6LzzzqObASBGxOQQlQJn6Js2bZp27NjhOwoAAFCgm5OTk7V7927fUQAAOSxmh2iXLl2UmpqqMWPG+I4CAAAU6OaDBw9q7NixvqMAAHJYzA7RJk2a6Oyzz9bw4cN9RwEAAJKaNWum4sWL080AEANCGqJm1sHMlpnZCjP7RyZfP8fMJpvZPDP72cw6ZX/U7BUXF6cuXbpo9OjROnDggO84AACclmjs5vj4eF122WUaNWqUDh065DsOACAHnXKImlmcpNckdZRUU1IvM6t5zNX6S/rcOXeBpJ6SXs/uoDnhiiuu0O7du5WcnOw7CgAAIYv2bt6xY4emTZvmOwoAIAeF8oxoY0krnHOrnHMHJX0q6fJjruMkFQl+XFTSH9kXMee0bt1aiYmJHAIEAIg0UdvNbdu2Vf78+elmAIhyoQzRcpLWZPh8bfCyjB6T1NvM1koaLenOzG7IzPqa2Wwzm7158+YziJu9ChQooPbt22vEiBFyzvmOAwBAqKK2mwsWLKi2bdvqm2++oZsBIIpl18mKekl61zlXXlInSR+Y2XG37Zwb4pxr6JxrWLJkyWz60VlzxRVXaN26dZozZ47vKAAAZKeI7ubff/9dCxYs8B0FAJBDQhmi6yQlZfi8fPCyjG6W9LkkOedmSMovqUR2BMxpnTt3VlxcnL755hvfUQAACFVUd/Nll10mM+PwXACIYqEM0VmSqplZJTPLq8AJD0Ycc53fJbWWJDOroUDZ+T++JwTFixfXRRddRNkBACJJVHfz2WefrWbNmvEgMQBEsVMOUedcqqQ7JI2VtESBM/AtMrMnzKxr8Gr3SrrVzBZI+kTSDS6CXthxxRVXaOHChVq5cqXvKAAAnFKsdPP8+fP122+/+Y4CAMgBIb1G1Dk32jlX3TlXxTn3VPCyR51zI4IfL3bONXfOne+cq+ecG5eTobPb5ZcHTjTIs6IAgEgRK908YsSxT/QCAKJBdp2sKKJVqlRJderU4RAgAADCRLVq1VSjRg19/fXXvqMAAHIAQzToqquu0vTp07Vx40bfUQAAgALdPHXqVG3ZssV3FABANmOIBnXr1k3OOZ4VBQAgTHTr1k2HDx/mpTMAEIUYokG1a9dW9erVNWzYMN9RAACApHr16qlSpUp0MwBEIYZokJmpW7dumjx5srZu3eo7DgAAMc/M1L17d02cOFE7duzwHQcAkI0Yohl0796dQ4AAAAgj3bp106FDhzRy5EjfUQAA2YghmsEFF1ygihUr6ssvv/QdBQAASGrcuLGSkpI4PBcAogxDNIP0Q4DGjx/PIUAAAIQBM9NVV12lsWPHavfu3b7jAACyCUP0GOmHAI0aNcp3FAAAoEA3HzhwQN9++63vKACAbMIQPUbjxo1Vvnx5DgECACBMNGvWTKVLl+alMwAQRRiix8iTJ4+uuuoqfffddxwCBABAGIiLi9OVV16p0aNHKyUlxXccAEA2YIhmonv37jpw4IBGjx7tOwoAAFCgm1NSUvTdd9/5jgIAyAYM0UykHwL0xRdf+I4CAAAktWjRQiVKlKCbASBKMEQzERcXp6uvvlrffvutdu3a5TsOAAAxLz4+Xt27d9eIESO0d+9e33EAAFnEED2Bnj17av/+/RoxYoTvKAAAQIFuTklJ4cz2ABAFGKIn0LRpU51zzjn67LPPfEcBAACSLrroIpUtW5ZuBoAowBA9gTx58qhHjx4aO3astm3b5jsOAAAxL/2lM6NHj9bOnTt9xwEAZAFD9CR69uypQ4cO6euvv/YdBQAAKNDNBw4c0PDhw31HAQBkAUP0JOrXr6+qVavq008/9R0FAABIatKkiSpUqEA3A0CEY4iehJnpmmuu0aRJk7Rp0ybfcQAAiHnp3Tx+/Hht3brVdxwAwBliiJ5Cz549lZaWpmHDhvmOAgAAFOjm1NRUffXVV76jAADOEEP0FGrXrq1atWpxCBAAAGGiXr16ql69Ot0MABGMIRqCa665RtOnT9fatWt9RwEAIOalH56bnJysDRs2+I4DADgDDNEQ9OrVS845HnkFACBM9OrVS2lpabynKABEKIZoCKpWraomTZrogw8+8B0FAABIqlGjhurXr083A0CEYoiGqE+fPvr555/1888/+44CAAAU6OY5c+ZoyZIlvqMAAE4TQzREPXr0UHx8vD766CPfUQAAgAJnz82TJw/dDAARiCEaopIlS6pDhw766KOPlJaW5jsOAAAxr3Tp0mrbti3dDAARiCF6Gnr37q1169YpOTnZdxQAAKBAN//666/6/vvvfUcBAJwGhuhp6Nq1qwoXLqwPP/zQdxQAACDpyiuvVMGCBelmAIgwDNHTUKBAAXXv3l3Dhg3Tvn37fMcBACDmFSxYUFdeeaU+//xzHThwwHccAECIGKKnqXfv3tq9e7dGjhzpOwoAAFCgm3fs2KHRo0f7jgIACBFD9DRdcsklKleuHIcAAQAQJlq3bq1SpUrRzQAQQRiipykuLk5//vOfNWbMGG3atMl3HAAAYl58fLyuvfZajRw5Ulu3bvUdBwAQAoboGbjhhhuUmprKI68AAISJG264QYcOHdLHH3/sOwoAIAQM0TNQo0YNNWnSRO+8846cc77jAAAQ8+rWrav69evrnXfe8R0FABAChugZuvHGG7Vw4ULNmTPHdxQAAKBAN8+bN08LFizwHQUAcAoM0TPUs2dP5c+fn0deAQAIE9dee63y5s1LNwNABGCInqGiRYvqqquu0scff6z9+/f7jgMAQMz705/+pMsvv1wffvihDh486DsOAOAkGKJZcMMNN2jHjh0aPny47ygAAECBw3O3bt3K+30DQJhjiGZBq1atlJSUxCFAAACEiXbt2qls2bJ0MwCEOYZoFsTFxen666/XuHHjtHbtWt9xAACIeXFxcbruuus0ZswYrV+/3nccAMAJMESz6IYbbpBzTu+//77vKAAAQIHDc9PS0vTBBx/4jgIAOAGGaBZVqVJFl1xyid566y2lpaX5jgMAQMyrXr26mjdvrqFDh/J+3wAQphii2aBv375atWqVJk2a5DsKAABQoJuXL1+uKVOm+I4CAMgEQzQbXHXVVSpevLiGDBniOwoAAJB09dVXq1ixYnQzAIQphmg2yJ8/v66//np9/fXX2rhxo+84AADEvAIFCqhPnz768ssvtWXLFt9xAADHYIhmk1tvvVWpqal67733fEcBAAAKdPPBgwc5aREAhCGGaDY577zz1KJFC7355puctAgAgDBQp04dXXjhhRoyZAgnLQKAMMMQzUZ9+/bVihUrlJyc7DsKAABQoJuXLl2q6dOn+44CAMggpCFqZh3MbJmZrTCzf5zgOj3MbLGZLTKzj7M3ZmTo1q2bzjrrLE6MAADIcXRzaHr06KGiRYvSzQAQZk45RM0sTtJrkjpKqimpl5nVPOY61SQ9JKm5c66WpL9nf9Twl37Soq+++kqbN2/2HQcAEKXo5tAlJiaqd+/e+uKLL7Rt2zbfcQAAQaE8I9pY0grn3Crn3EFJn0q6/Jjr3CrpNefcdklyzm3K3piR49Zbb9WhQ4f09ttv+44CAIhedPNp6Nu3rw4cOKB3333XdxQAQFAoQ7ScpDUZPl8bvCyj6pKqm9n3ZjbTzDpkV8BIU7NmTbVs2VL//e9/dfjwYd9xAADRiW4+DXXr1lXz5s31+uuvc0JBAAgT2XWyonhJ1SS1lNRL0ptmVuzYK5lZXzObbWazo/nQ1X79+um3337T6NGjfUcBAMQuujmDfv36aeXKlRo3bpzvKAAAhTZE10lKyvB5+eBlGa2VNMI5d8g5t1rS/xQov6M454Y45xo65xqWLFnyTDOHvcsvv1xly5bVa6+95jsKACA60c2nqVu3bipVqhTdDABhIpQhOktSNTOrZGZ5JfWUNOKY63yjwCOuMrMSChwOtCr7YkaWhIQE3XbbbRo7dqyWL1/uOw4AIPrQzacpb9686tu3r7799lutXr3adxwAiHmnHKLOuVRJd0gaK2mJpM+dc4vM7Akz6xq82lhJW81ssaTJku53zm3NqdCR4NZbb1V8fLxef/1131EAAFGGbj4zt912m/LkyaP//ve/vqMAQMwz55yXH9ywYUM3e/ZsLz87t/Tq1UtjxozRunXrVLBgQd9xACCqmdkc51xD3zkiWSx0c/fu3TV58mStXbtWBQoU8B0HAKLaybo5u05WhEz069dPO3fu1Mcfx+R7iAMAEHbuuOMObdu2TZ999pnvKAAQ0xiiOah58+aqW7euXnvtNfl65hkAAPy/Sy65RLVq1dKgQYPoZgDwiCGag8xMd9xxhxYsWKCpU6f6jgMAQMwzM/Xr109z5szRjBkzfMcBgJjFEM1hvXv3VvHixfXSSy/5jgIAACRdd911KlasGN0MAB4xRHNYgQIFdPvtt2vEiBFauXKl7zgAAMS8ggULqm/fvvrqq6/022+/+Y4DADGJIZoL/vrXvyo+Pl6vvPKK7ygAAECBkxaZmV599VXfUQAgJjFEc0HZsmV1zTXX6O2339bOnTt9xwEAIOYlJSXp6quv1ptvvqndu3f7jgMAMYchmkvuvvtu7dmzR0OHDvUdBQAAKNDNu3bt0jvvvOM7CgDEHIZoLqlfv75atGihV199Vampqb7jAAAQ8xo3bqwLL7xQr7zyig4fPuw7DgDEFIZoLrr77rv122+/6ZtvvvEdBQAAKNDNK1eu1KhRo3xHAYCYwhDNRV26dFHlypX14osv8ibaAACEgSuvvFLnnHOOXnzxRd9RACCmMERzUVxcnO655x7NnDlT33//ve84AADEvPj4eN19992aNm2aZs6c6TsOAMQMhmguu/HGG1WiRAk999xzvqMAAABJt9xyi8466ywNHDjQdxQAiBkM0VyWmJioO++8U6NGjdKiRYt8xwEAIOYVKlRI/fr10zfffKNly5b5jgMAMYEh6kG/fv2UmJio559/3ncUAAAg6c4771S+fPn0wgsv+I4CADGBIepB8eLFdcstt+ijjz7S2rVrfccBACDmnX322brxxhv1/vvva/369b7jAEDUY4h6cs8998g5p5deesl3FAAAIOnee+9VamqqXn75Zd9RACDqMUQ9qVChgnr27KkhQ4Zo+/btvuMAABDzqlSpou7du2vw4MHauXOn7zgAENUYoh498MAD2rNnj1577TXfUQAAgKQHH3xQu3bt0n//+1/fUQAgqjFEPapbt646d+6sl19+WXv27PEdBwCAmFe/fn21a9dO//73v5WSkuI7DgBELYaoZ4888oi2bt2qwYMH+44CAAAU6ObNmzdryJAhvqMAQNRiiHrWpEkTtW3bVi+88IL27dvnOw4AADHvoosuUsuWLTVw4EDt37/fdxwAiEoM0TDQv39/bdy4UW+++abvKAAAQIFnRdevX6933nnHdxQAiEoM0TDQokULtWjRQgMHDtSBAwd8xwEAIOZdeumluvDCC/Xss8/q4MGDvuMAQNRhiIaJ/v37a926dXr33Xd9RwEAIOaZmR555BH9/vvv+uCDD3zHAYCowxANE23atFGTJk30zDPP6NChQ77jAAAQ8zp06KAGDRro6aefVmpqqu84ABBVGKJhIv2R199++41nRQEACAPp3bxq1Sp9+OGHvuMAQFRhiIaRTp06qUmTJhowYACvFQUAIAx07dpVDRo00BNPPMERSwCQjRiiYcTM9MQTT2jNmjUaOnSo7zgAAMS89G5evXo1Z9AFgGzEEA0zbdu21UUXXaSnnnqK9xUFACAMdOzYUU2bNtWTTz7JEUsAkE0YomHGzDRgwACtX79egwcP9h0HAICYl97Na9as4T2/ASCbMETDUMuWLdWqVSs9++yz2rt3r+84AADEvNatW6tFixZ66qmnlJKS4jsOAEQ8hmiYGjBggDZt2qRBgwb5jgIAQMxLf1Z0w4YN+u9//+s7DgBEPIZomGrWrJk6duyo5557Tjt27PAdBwCAmNeiRQu1bdtWzz77rHbt2uU7DgBENIZoGHvqqae0fft2DRw40HcUAACgQDdv2bJFL7zwgu8oABDRGKJh7IILLtC1116rl19+WevWrfMdBwCAmNeoUSP16NFDL774ojZs2OA7DgBELIZomBswYIBSU1P1+OOP+44CAAAkPfnkkzp48KAGDBjgOwoARCyGaJirXLmy/vKXv+itt97S0qVLfccBACDmVatWTX379tWQIUO0fPly33EAICIxRCNA//79VbBgQT388MO+owAAAEmPPPKI8uXLp/79+/uOAgARiSEaAUqWLKn7779fX3/9tWbMmOE7DgAAMa906dK699579fnnn2vWrFm+4wBAxGGIRoi7775bpUuX1n333SfnnO84AADEvHvvvVclS5akmwHgDDBEI0ShQoU0YMAA/fDDD/riiy98xwEAIOYVKVJEjz/+uKZOnapvvvnGdxwAiCjm6xG8hg0butmzZ3v52ZHq8OHDatCggXbu3KklS5Yof/78viMBQNgwsznOuYa+c0Qyuvn0paamql69etq/f78WLVqkfPny+Y4EAGHjZN3MM6IRJC4uTv/+97/166+/6uWXX/YdBwCAmBcfH68XX3xRK1eu1KBBg3zHAYCIwRCNMK1atVLXrl319NNPa+PGjb7jAAAQ89q3b6+OHTtqwIAB2rx5s+84ABARGKIR6Pnnn9e+ffv0yCOP+I4CAAAkvfjii9qzZ48ee+wx31EAICIwRCNQ9erVdccdd+itt97SggULfMcBACDm1ahRQ7fffrveeOMNLVq0yHccAAh7DNEI9eijj+pPf/qT7rjjDk4ZDwBAGHjsscdUpEgR3XnnnXQzAJwCQzRCnXXWWXrmmWc0ffp0ffzxx77jAAAQ80qUKKGnnnpKkydP5q3WAOAUePuWCJaWlqamTZtqzZo1WrZsmYoUKeI7EgB4w9u3ZB3dnHWHDx9Wo0aNtGnTJi1dulSFChXyHQkAvOHtW6JUnjx59Nprr2njxo164oknfMcBACDmxcXFadCgQVq3bp2eeuop33EAIGyFNETNrIOZLTOzFWb2j5Ncr5uZOTPjEelc0qhRI9188836z3/+o8WLF/uOAwDIJXRz+GrWrJmuv/56vfjii1q2bJnvOAAQlk45RM0sTtJrkjpKqimpl5nVzOR6hSXdJenH7A6Jk3vmmWdUuHBhTlwEADGCbg5/zz33nAoUKKC//e1vdDMAZCKUZ0QbS1rhnFvlnDso6VNJl2dyvQGSnpO0PxvzIQQZT47AiYsAICbQzWGuVKlSGjBggMaNG8eJiwAgE6EM0XKS1mT4fG3wsiPMrL6kJOfct9mYDaehb9++atKkie6++25t27bNdxwAQM6imyNAv3791KBBA911113asWOH7zgAEFayfLIiM8sj6d+S7g3hun3NbLaZzd68eXNWfzQyiIuL05AhQ7Rt2zY98MADvuMAADyim8NDejdv2rRJDz30kO84ABBWQhmi6yQlZfi8fPCydIUl1ZaUbGa/SmoqaURmJ0Vwzg1xzjV0zjUsWbLkmadGpurWrat7771Xb731lqZOneo7DgAg59DNEaJ+/fq66667NHjwYM2YMcN3HAAIG6d8H1Ezi5f0P0mtFSi5WZKudc4tOsH1kyXd55w76RuR8V5lOWPv3r2qXbu28ufPr/nz5ytfvny+IwFAroil9xGlmyPLnj17VLNmTRUtWlRz585VQkKC70gAkCuy9D6izrlUSXdIGitpiaTPnXOLzOwJM+uavVGRVQULFtTrr7+upUuX6rnnnvMdBwCQA+jmyFKoUCENGjRICxcu1AsvvOA7DgCEhVM+I5pTeNQ1Z/Xq1UtfffWV5s6dq1q1avmOAwA5LpaeEc0pdHPO6t69u0aNGqX58+frvPPO8x0HAHJclp4RRWR65ZVXVKRIEd14441KTU31HQcAgJg3aNAgFSxYUDfddJMOHz7sOw4AeMUQjVIlS5bUoEGDNGvWLL300ku+4wAAEPNKly6t//znP5oxY4ZeffVV33EAwCuGaBTr0aOHrrjiCj3yyCNatmyZ7zgAAMS8P//5z+rcubMefvhhrVixwnccAPCGIRrFzEyvv/66EhMTdfPNN3MYEAAAnpmZ3njjDeXNm1e33HKL0tLSfEcCAC8YolGuTJkyevnll/X999/rlVde8R0HAICYV65cOf373//WlClT9Nprr/mOAwBeMERjQJ8+fdSlSxc99NBDWrx4se84AADEvBtvvFEdO3bUgw8+yMtnAMQkhmgMMDO9+eabKlKkiHr37q2DBw/6jgQAQEwzM7311lsqUKCA+vTpo0OHDvmOBAC5iiEaI0qVKqUhQ4Zo3rx5euKJJ3zHAQAg5pUpU0aDBw/WrFmz9PTTT/uOAwC5iiEaQ6644grdcMMNeuaZZzRjxgzfcQAAiHlXX321evfurQEDBuinn37yHQcAcg1DNMb85z//UVJSkvr06aM9e/b4jgMAQMx79dVXVaZMGfXp00d79+71HQcAcgVDNMYUKVJE77//vlatWqU777zTdxwAAGJesWLF9N5772n58uX6+9//7jsOAOQKhmgMatGihfr37693331XH3/8se84AADEvFatWunBBx/U0KFD9fnnn/uOAwA5jiEaox599FE1b95ct99+u1auXOk7DgAAMe+JJ55QkyZN1LdvX/3666++4wBAjmKIxqj4+Hh99NFHiouLU69evXhLFwAAPEtISNAnn3wi55x69erFW7oAiGoM0RhWoUIFDR06VLNmzVL//v19xwEAIOZVqlRJQ4YM0cyZM/Wvf/3LdxwAyDEM0RjXrVs33X777Xr++ec1fPhw33EAAIh511xzjW6++WY988wz+vbbb33HAYAcwRCFXnrpJTVo0EDXX3+9VqxY4TsOAAAx79VXX1W9evXUp08frV692nccAMh2DFEof/78GjZsmOLi4tStWzelpKT4jgQAQEwrUKCAvvzySznn1K1bN+3fv993JADIVgxRSJIqVqyojz76SL/88ov++te/yjnnOxIAADGtcuXK+uCDDzRv3jze+xtA1GGI4ogOHTro0Ucf1Xvvvac333zTdxwAAGLeZZddpn/+858aOnSo3n77bd9xACDbMERxlEcffVQdOnTQnXfeqVmzZvmOAwBAzHv88cfVpk0b9evXT/PmzfMdBwCyBUMUR8mTJ48+/PBDlSlTRt27d9eWLVt8RwIAIKbFxcXp448/VokSJdStWzdt377ddyQAyDKGKI5TvHhxDRs2TBs3btTVV1+tgwcP+o4EAEBMK1mypIYNG6a1a9eqR48eOnTokO9IAJAlDFFkqmHDhho6dKiSk5N1xx13cPIiAAA8a9KkiYYMGaIJEybo7rvv9h0HALIk3ncAhK/evXtryZIlevrpp1WrVi3dddddviMBABDTbrjhBi1ZskQDBw5UjRo11K9fP9+RAOCMMERxUgMGDNCSJUt0zz33qHr16urYsaPvSAAAxLSnn35aS5cu1V133aXq1aurbdu2viMBwGnj0FycVJ48efTBBx+obt26uuaaa7Ro0SLfkQAAiGlxcXH68MMPVatWLV199dVaunSp70gAcNoYojilggULasSIESpYsKC6dOnCmXQBAPCscOHCGjFihPLly6fLLrtMW7du9R0JAE4LQxQhSUpK0vDhw7V+/XpdddVV2r9/v+9IAADEtAoVKuibb77RmjVr1L17dx04cMB3JAAIGUMUIWvcuLHeffddTZs2Tb1799bhw4d9RwIAIKZdeOGFeuedd5ScnKzrrrtOaWlpviMBQEgYojgt11xzjV566SV9+eWX+tvf/sbbugAA4Nm1116rgQMH6vPPP9fdd99NNwOICJw1F6ft73//u/744w89//zzKlu2rP75z3/6jgQAQEy77777tH79er300ksqW7asHnzwQd+RAOCkGKI4I88++6w2bNig/v37q3Tp0rr55pt9RwIAIGaZmV544QVt2LBB//jHP1S6dGldf/31vmMBwAkxRHFG8uTJo7feekubN29W3759dfbZZ6tLly6+YwEAELPy5Mmjd999V5s3b9bNN9+skiVLqlOnTr5jAUCmeI0ozlhCQoK++OILNWjQQD169ND333/vOxIAADEtb968+uqrr3T++efr6quv1syZM31HAoBMMUSRJYUKFdK3336rc845R506ddKsWbN8RwIAIKYVLlxYo0ePVpkyZdShQwfNmTPHdyQAOA5DFFlWsmRJTZw4USVKlFC7du00b94835EAAIhppUqV0qRJk1SsWDG1a9dOP//8s+9IAHAUhiiyRfny5TVp0iQVKVJEbdu21cKFC31HAgAgpp1zzjmaNGmSEhMT1aZNGy1evNh3JAA4giGKbFOhQgVNmjRJ+fLlU+vWrbV06VLfkQAAiGmVK1fWxIkTFRcXp9atW2v58uW+IwGAJIYoslmVKlU0adIkmZlatWqlFStW+I4EAEBMq169uiZOnKjDhw+rVatWWr16te9IAMAQRfY799xzNXHiRB06dIgxCgBAGKhZs6YmTJiglJQUXXrppVq1apXvSABiHEMUOaJWrVqaMGGC9u3bp4svvliLFi3yHQkAgJhWt25djR8/Xrt371aLFi14CQ0ArxiiyDHnn3++pkyZIjPTJZdcorlz5/qOBABATKtfv76Sk5OVmpqqFi1aaMGCBb4jAYhRDFHkqJo1a2ratGkqVKiQLr30Uv3www++IwEAENPq1KmjqVOnKl++fGrZsqV+/PFH35EAxCCGKHJclSpVNG3aNJUqVUrt2rXTpEmTfEcCACCmVa9eXdOmTVPx4sXVpk0bTZkyxXckADGGIYpckZSUpKlTp6pSpUrq1KmTRo4c6TsSAAAxrWLFipo6daqSkpLUoUMHjRkzxnckADGEIYpcU7p0aSUnJ6tOnTq68sor9eabb/qOBABATCtbtqymTJmiGjVqqEuXLnrnnXd8RwIQIxiiyFXFixfX5MmT1a5dO/Xt21ePPvqonHO+YwEAELNKliyp5ORktWrVSjfddJOeeOIJuhlAjmOIItcVKlRIw4cP10033aQBAwbopptu0qFDh3zHAgAgZhUpUkTffvutrrvuOv3rX//SbbfdptTUVN+xAESxeN8BEJsSEhI0dOhQnXPOOXrssce0fv16ffHFFypcuLDvaAAAxKSEhAS9++67SkpK0lNPPaU//vhDn332mQoWLOg7GoAoxDOi8MbM9K9//UtDhw7VhAkTdMkll2jdunW+YwEAELPMTE8++aQGDx6sMWPGqGXLllq/fr3vWACiUEhD1Mw6mNkyM1thZv/I5Ov3mNliM/vZzCaaWYXsj4podfPNN2vkyJFavny5GjVqxPuZAUAI6GbkpNtuu03ffPONFi9erEaNGmn27Nm+IwGIMqccomYWJ+k1SR0l1ZTUy8xqHnO1eZIaOufqShomaWB2B0V069ixo2bMmKH8+fPrkksu0QcffOA7EgCELboZuaFLly764YcfFB8fr4svvliffPKJ70gAokgoz4g2lrTCObfKOXdQ0qeSLs94BefcZOdcSvDTmZLKZ29MxILatWtr1qxZatasma677jo98MADOnz4sO9YABCO6GbkivPPP1+zZs1S48aNde211+rhhx9WWlqa71gAokAoQ7ScpDUZPl8bvOxEbpaU6Tsim1lfM5ttZrM3b94cekrEjOLFi2vs2LHq16+fnn/+eXXt2lU7d+70HQsAwg3djFxTsmRJjR8/Xn379tUzzzyjK664Qrt27fIdC0CEy9aTFZlZb0kNJT2f2dedc0Occw2dcw1LliyZnT8aUSQhIUGDBg3S4MGDNW7cODVu3Fg///yz71gAEJHoZmSHvHnzavDgwXrttdc0evRoNWnSRIsWLfIdC0AEC2WIrpOUlOHz8sHLjmJmbST9U1JX59yB7ImHWHbbbbdp0qRJ2r17t5o0aaJ33nnHdyQACBd0M3Kdmemvf/2rJkyYoO3bt6tx48ac0wHAGQtliM6SVM3MKplZXkk9JY3IeAUzu0DSGwoU3absj4lYdfHFF2vevHlq1qyZbrrpJt10001KSUk59TcCQHSjm+FNy5YtNW/ePDVq1EjXXXedbr31Vu3bt893LAAR5pRD1DmXKukOSWMlLZH0uXNukZk9YWZdg1d7XlIhSV+Y2XwzG3GCmwNOW6lSpTRu3Dg98sgjevfdd9W0aVMtW7bMdywA8IZuhm9lypTRhAkT9PDDD2vo0KFq1qyZVqxY4TsWgAhizjkvP7hhw4aO96TC6fruu+/Uu3dvHThwQK+//rp69+4tM/MdC0AYMLM5zrmGvnNEMroZZ2L06NHq06ePDh06pDfeeEO9evXyHQlAmDhZN2fryYqAnNahQwfNmzdP9erV03XXXadevXpp+/btvmMBABCzOnXqpLlz56p27dq69tpr1bt3b+3YscN3LABhjiGKiJOUlKTk5GQ99dRT+vLLL1W3bl1NnjzZdywAAGJWhQoVNHXqVD3++OP69NNPdf7552vq1Km+YwEIYwxRRKS4uDg9/PDD+uGHH5SYmKjWrVvrgQce0IEDnBQSAAAf4uPj9eijj2r69OlKSEhQy5Yt9dBDD+ngwYO+owEIQwxRRLRGjRpp7ty5uvXWW/X888+rSZMmmjdvnu9YAADErKZNm2r+/Pm66aab9Oyzz6pp06ZasGCB71gAwgxDFBGvYMGCeuONNzR8+HBt3LhRjRo10j//+U/t37/fdzQAAGJSoUKFNHToUH311Vdat26dGjZsqEcffZQjlwAcwRBF1OjatasWL16sPn366Omnn1b9+vU1c+ZM37EAAIhZV155pRYvXqyePXtqwIABatCggX766SffsQCEAYYoospZZ52ld955R2PGjNGePXvUrFkz3XPPPdq7d6/vaAAAxKTixYvrgw8+0KhRo7Rjxw5deOGFuv/++5WSkuI7GgCPGKKISh06dNDChQt1++2366WXXlLNmjU1fPhw+XrfXAAAYl3nzp21aNEi3XLLLXrhhRdUq1YtjRw50ncsAJ4wRBG1ihQpotdff11Tp05VkSJFdMUVV6hLly5atWqV72gAAMSkokWL6o033lBycrISExPVtWtXXX755fr11199RwOQyxiiiHoXX3yx5s6dqxdeeEFTpkxRrVq1NGDAAE5mBACAJ5dcconmz5+vgQMHauLEiapZs6aeeuopTmYExBCGKGJCQkKC7r33Xi1dulRdu3bVo48+qjp16uibb77hcF0AADxISEjQ/fffryVLlqhTp07q37+/6tSpo5EjR9LNQAxgiCKmlCtXTp999pnGjRunvHnz6sorr1SrVq00d+5c39EAAIhJSUlJGjZsmL777jvFxcWpa9euatOmjebPn+87GoAcxBBFTGrbtq0WLFig119/XQsXLlTDhg114403at26db6jAQAQk9q3b6+ff/5ZgwYN0oIFC1S/fn3dcsstWr9+ve9oAHIAQxQxKz4+Xn/5y1+0YsUK3X///fr4449VvXp1PfLII9q5c6fveAAAxJyEhAT169dPy5cv1z333KP3339f1apV02OPPaZdu3b5jgcgGzFEEfOKFi2q5557TkuXLlWXLl305JNPqlKlSnruued4/1EAADw466yz9MILL2jJkiXq2LGjHn/8cVWuXFnPP/887z8KRAmGKBBUqVIlffrpp5o7d66aNWumf/zjH6pataoGDRrEWfwAAPCgSpUq+uKLLzR79mw1atRIDzzwgKpWrarXX39dBw8e9B0PQBYwRIFjXHDBBRo1apSmT5+uc889V3feeafOPfdcDR48mLd8AQDAgwYNGmjMmDGaOnWqqlSpon79+uncc8/VkCFDeLAYiFAMUeAEmjdvrsmTJ2vcuHEqU6aM/vKXv6hKlSp6+eWXOSwIAAAPLr74Yk2dOlVjxozR2Wefrdtuu01VqlTRK6+8QjcDEYYhCpyEmalt27b64YcfNHHiRFWvXl133323KlasqGeffZYTJwAAkMvMTB06dNDMmTM1btw4ValSRXfddZcqVaqkgQMHavfu3b4jAggBQxQIgZmpVatWmjx5sqZPn66GDRvqoYceUlJSku677z6tWbPGd0QAAGJK+oPFU6ZM0ZQpU1SvXj09+OCDSkpK0gMPPKC1a9f6jgjgJBiiwGlq3ry5Ro8erdmzZ6tTp056+eWXValSJV177bWaM2eO73gAAMScFi1aaOzYsfrpp5/Uvn17vfjii6pUqZJ69+6tefPm+Y4HIBMMUeAMNWjQQJ988olWrlypu+66S6NGjVLDhg3VsmVLDRs2TIcOHfIdEQCAmNKoUSN99tlnWrlype644w4NHz5c9evXV6tWrfTVV18pNTXVd0QAQQxRIIsqVKigF198UWvWrNELL7yg1atX6+qrr1aFChX0r3/9i0ODAADIZRUrVtRLL72kNWvWaODAgVqxYoW6deumihUr6vHHH9cff/zhOyIQ8xiiQDYpWrSo7r33Xq1atUojR45UvXr1NGDAAFWsWFFXXXWVJkyYoLS0NN8xAQCIGcWKFdP999+vVatWafjw4apdu7Yee+wxnXPOOerevbsmTZok55zvmEBMYogC2SwuLk6XXXaZRo8erRUrVui+++7TtGnT1LZtW5133nl6/vnntX79et8xAQCIGfHx8eratau+++47rVixQvfcc4+Sk5PVunVr1ahRQy+++KI2bNjgOyYQUxiiQA6qXLmynn32Wa1du1Yffvihzj77bD3wwAMqX768OnXqpM8++0z79+/3HRMAgJhRpUoVDRw4UGvXrtX777+vP/3pT7rvvvtUvnx5de7cWV988QXdDOQChiiQC/Lly6c///nPmj59upYtW6aHHnpICxcuVM+ePVW6dGnddtttmjFjBocHAQCQS/Lnz68+ffrohx9+0JIlS/TAAw/o559/Vo8ePVSmTBn95S9/0cyZM+lmIIeYr//natiwoZs9e7aXnw2Eg7S0NE2ePFnvvfeevvzyS6WkpKhatWq65ppr1KNHD9WuXVtm5jsmEDHMbI5zrqHvHJGMbkasO3z4sCZNmqT33ntPX331lfbt26fq1aurZ8+e6tGjh2rVquU7IhBRTtbNDFEgDOzevVvDhg3Thx9+qOTkZKWlpalGjRrq0aOHevTooZo1a/qOCIQ9hmjW0c3A/9u1a9eRbp4yZYrS0tJUs2bNI91co0YN3xGBsMcQBSLIpk2b9NVXX+nzzz9XcnKynHOqVauWevTooauuukq1atXimVIgEwzRrKObgcxt2LDhSDdPnTpVzjnVqVNHV199ta666irVrFmTbgYywRAFItSGDRv05Zdf6vPPP9e0adPknFOlSpXUtWtXde3aVRdffLESEhJ8xwTCAkM06+hm4NTWr19/pJunT58u55wqV658pJsvuugiuhkIYogCUWD9+vUaNWqURowYoQkTJmj//v0qWrSoOnbsqK5du6pDhw4666yzfMcEvGGIZh3dDJyeP/7446huPnDggIoVK6ZOnTod6eaiRYv6jgl4wxAFokxKSoomTJigESNGaOTIkdq0aZPy5MmjJk2aqH379mrXrp0aNWqk+Ph431GBXMMQzTq6GThze/fu1fjx4zVixAiNGjVKmzdvVlxcnJo2bap27dqpffv2atiwoeLi4nxHBXINQxSIYmlpafrxxx81ZswYjR07VrNmzZJzTsWKFVPr1q2PDNMKFSr4jgrkKIZo1tHNQPY4fPjwUd08e/ZsOed01llnqU2bNmrXrp3atWunc845x3dUIEcxRIEYsm3bNk2YMEHjxo3T2LFjtXbtWklS5cqV1bJlS11yySVq2bIl5YeowxDNOroZyBlbtmzRxIkTj3TzunXrJElVqlQ5qpuTkpI8JwWyF0MUiFHOOS1dulTjxo3TlClTNGXKFG3btk2SVKlSpSPFd/HFF6tSpUqc8Q8RjSGadXQzkPOcc1qyZInGjRun5ORkTZ06Vdu3b5cUeNA4YzdXrFiRbkZEY4gCkBQ4jHfhwoVKTk4+Mky3bt0qSTr77LN14YUXHvnTsGFDJSYmek4MhI4hmnV0M5D70tLS9Msvvyg5OfnIME1/0LhUqVLHdXOBAgU8JwZCxxAFkKn0YfrDDz9oxowZmjFjhpYvXy5JiouL0/nnn3+k/Jo2barKlSvzyCzCFkM06+hmwL/0Yfr9999r5syZmjFjhlasWCFJio+PP66bOaIJ4YwhCiBkW7ZsOVJ8M2bM0E8//aS9e/dKkooVK6YLLrhADRo0UP369VW/fn1Vq1ZNefLk8ZwaYIhmB7oZCE+bN28+rptTUlIkBbo5vZPT+7lq1ap0M8ICQxTAGTt8+LAWLlyoH3/8UfPmzdOcOXP0888/68CBA5KkQoUK6YILLjhSgnXq1FGNGjWUP39+z8kRaxiiWUc3A5EhNTU1024+ePCgpP/v5gYNGuiCCy5Q3bp1dd5559HNyHUMUQDZ6tChQ1qyZInmzp2ruXPnas6cOZo/f/6RR2fz5MmjatWqqXbt2kf+1KlTR1WqVOG9TZFjGKJZRzcDkevQoUNavHjxcd28b98+SYGX3BzbzbVr11bVqlV5b1PkGIYogBx3+PBh/e9//9PChQu1cOFC/fLLL1q4cKFWrFih9P+dyZcvn2rUqKEaNWqoevXqOvfcc3XuueeqWrVqKly4sOffAJGOIZp1dDMQXdK7Ob2T0/t55cqVR3VzzZo1j+rm6tWrq3r16ipUqJDn3wCRjiEKwJuUlBQtWbLkqAJcunSpfv/9d2X835+yZcseKb+Mf1esWJFnUREShmjW0c1AbEjv5vSB+ssvv2jZsmXHdXO5cuWO6+Vzzz1XFSpUoJsREoYogLCzb98+rVixQv/73/+0bNmyI38vW7bsyPupSYFDiZKSklSpUiVVqlRJlStXPvJxpUqVVKpUKc4WCEkM0exANwOxLb2bM/Zy+t/HdvM555xzVB9n7Oizzz6bboakk3czD2UA8KJAgQKqU6eO6tSpc9zXtmzZcqT4Vq5cqdWrV2v16tX69ttvtXHjxuNup2LFiqpUqZIqVqyo8uXLKykpSeXLl1f58uVVrlw53nMNAIAQnKibnXPaunXrkQeMV61apVWrVmn16tUaNWrUcd2cmJh4XDen/0lKSlK5cuU4cRIYogDCT4kSJVSiRAk1a9bsuK+lpKTo119/PTJO04tw9erVmjFjxlGP2Ga8vYwlmHGsli1bVqVLl1bRokV59BYAgEyY2ZFubt68+XFf37t371HdnLGjv//+e+3YseO470nv5owPHmd8ELlUqVJ0c5RjiAKIKImJiapZs6Zq1qyZ6df37NmjdevWae3atUf9WbNmjdasWaMZM2Zo69atx31fvnz5VKpUKZUuXfrI3xk/zngZJ28AAOD/FSxYULVq1VKtWrUy/XrGbl6zZs1R/fz777/rhx9+OGE3n6iLj72Mbo48DFEAUaVQoUJHTqZwIvv27dO6deu0Zs0arV+/Xhs2bNCGDRu0ceNGbdiwQb///rt++uknbdq0SZm9jj5//vxHHhkuXrz4Sf9O/7hgwYI8qgsAiEmhdHNKSspR3Zzeyen9/Ouvv2rmzJnavHlzpt1coECBE3bwifo5MTGRbvaIIQog5hQoUEBVq1ZV1apVT3q9w4cPa8uWLccN1c2bN2vLli3aunWrtmzZonnz5mnLli3avn17puUoBR7VLV68uP70pz+pWLFiKlq0qIoVK5bpn2O/VrRoUeXNmzcn7goAAMJCYmKiqlWrpmrVqp30eqmpqdqyZctRQzW9m9N7eevWrZo7d662bt2qbdu2nfC28uXLpxIlSuiss846aQ9ndjndnHUMUQA4gbi4OJUqVUqlSpXS+eeff8rrHz58WNu3bz9qpB779/bt27Vjxw6tX79eS5Ys0Y4dO7Rjxw6lpaWd9LYTExOPFGDhwoVVqFChU/59oq8VLFiQNy8HAESk+Pj4I4fnnm43Z9bLW7ZsOdLF69at06JFi7Rjxw7t3LkzpG5OH6WhdvOJ/o7FbmaIAkA2iYuLO3K4z+lwzmnPnj1HijC9ADN+nvHPnj17tHv3bm3ZskW7d+8+8vn+/ftD/pmJiYkqWLCgEhMTVaBAASUmJh71caiXHfv1/PnzK1++fMqfP/+Rj3mvOQCAL9nZzSfr6Nzo5jPp6fRuztjP4dLNISUwsw6S/iMpTtJQ59yzx3w9n6T3JTWQtFXSNc65X7M3KgBEJzNT4cKFVbhwYSUlJZ3x7aSmpmrv3r1HFeDJ/t6zZ4/27dunffv2KSUlRfv27dPu3bu1adMmpaSkHLksJSXltIr0WHFxcceN04x/v/3226pRo8YZ336sopsBIOdkZzend27G/j3dbt64ceNRl2VHN2fWyel/v//++6c8TDqrTjlEzSxO0muS2kpaK2mWmY1wzi3OcLWbJW13zlU1s56SnpN0TU4EBgBkLj4+XkWLFlXRokWz/bbT0tK0f//+40ow48f79u3T/v37deDAAe3fvz/kjxMSErI9b7SjmwEgMsTHxx95XWl2S+/mY7v42J5O79xj//bdzaE8I9pY0grn3CpJMrNPJV0uKWPZXS7pseDHwyQNMjNzJzprBwAgouTJk+fI4T7Fixf3HQd0MwDEvIzdHInyhHCdcpLWZPh8bfCyTK/jnEuVtFPScf9SMbO+ZjbbzGZv3rz5zBIDAAC6GQAQ0UIZotnGOTfEOdfQOdewZMmSufmjAQBAJuhmAIAPoQzRdZIyvkK3fPCyTK9jZvGSiipwYgQAAJD96GYAQEQLZYjOklTNzCqZWV5JPSWNOOY6IyRdH/y4u6RJvAYFAIAcQzcDACLaKU9W5JxLNbM7JI1V4BTxbzvnFpnZE5JmO+dGSHpL0gdmtkLSNgUKEQAA5AC6GQAQ6UJ6H1Hn3GhJo4+57NEMH++XdHX2RgMAACdCNwMAIlmunqwIAAAAAACGKAAAAAAgVzFEAQAAAAC5iiEKAAAAAMhVDFEAAAAAQK5iiAIAAAAAchVDFAAAAACQq8w55+cHm22W9Fs23VwJSVuy6baiFfdRaLifQsP9FBrup1PLzvuognOuZDbdVkyim3Md91FouJ9Cw/0UGu6nU8uVbvY2RLOTmc12zjX0nSOccR+FhvspNNxPoeF+OjXuo+jFf9tT4z4KDfdTaLifQsP9dGq5dR9xaC4AAAAAIFcxRAEAAAAAuSpahugQ3wEiAPdRaLifQsP9FBrup1PjPope/Lc9Ne6j0HA/hYb7KTTcT6eWK/dRVLxGFAAAAAAQOaLlGVEAAAAAQISI6CFqZh3MbJmZrTCzf/jOE47M7G0z22RmC31nCWdmlmRmk81ssZktMrO7fGcKN2aW38x+MrMFwfvocd+ZwpmZxZnZPDMb5TtLuDKzX83sFzObb2azfedB9qCbT41uDg3dfGp08+mhm08tN7s5Yg/NNbM4Sf+T1FbSWkmzJPVyzi32GizMmFkLSXskve+cq+07T7gyszKSyjjn5ppZYUlzJF3B/z39PzMzSQWdc3vMLEHSdEl3Oedmeo4WlszsHkkNJRVxzl3mO084MrNfJTV0zvF+blGCbg4N3RwauvnU6ObTQzefWm52cyQ/I9pY0grn3Crn3EFJn0q63HOmsOOcmyppm+8c4c45t945Nzf48W5JSySV85sqvLiAPcFPE4J/IvORrBxmZuUldZY01HcWIJfRzSGgm0NDN58a3Rw6ujn8RPIQLSdpTYbP14r/cUI2MLOKki6Q9KPnKGEneEjLfEmbJI13znEfZe5lSQ9ISvOcI9w5SePMbI6Z9fUdBtmCbkaOoJtPjG4O2cuim0ORa90cyUMUyHZmVkjSl5L+7pzb5TtPuHHOHXbO1ZNUXlJjM+OQsmOY2WWSNjnn5vjOEgEucs7Vl9RRUr/g4YoAcBS6+eTo5lOjm09LrnVzJA/RdZKSMnxePngZcEaCr634UtJHzrmvfOcJZ865HZImS+rgOUo4ai6pa/A1Fp9KamVmH/qNFJ6cc+uCf2+S9LUCh3UistHNyFZ0c+jo5pOim0OUm90cyUN0lqRqZlbJzPJK6ilphOdMiFDBF/u/JWmJc+7fvvOEIzMraWbFgh8XUOBkJEu9hgpDzrmHnHPlnXMVFfjfpUnOud6eY4UdMysYPPmIzKygpHaSOINo5KObkW3o5lOjm0NDN4cmt7s5Yoeocy5V0h2Sxirw4vXPnXOL/KYKP2b2iaQZks41s7VmdrPvTGGquaQ+CjxCNj/4p5PvUGGmjKTJZvazAv/YHO+c4/TnOFOlJE03swWSfpL0rXPuO8+ZkEV0c2jo5pDRzadGNyM75Wo3R+zbtwAAAAAAIlPEPiMKAAAAAIhMDFEAAAAAQK5iiAIAAAAAchVDFAAAAACQqxiiAAAAAIBcxRAFAAAAAOQqhigAAAAAIFcxRAEAAAAAuer/ABtzYQmvSPeLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(16, 6))\n",
    "axs[0].plot(t, C_num, label=\"numerical solution\", c=\"black\")\n",
    "axs[1].plot(t, C_exact, label=\"exact solution\", c=\"black\")\n",
    "axs[0].set_title(\"numerical solution\")\n",
    "axs[1].set_title(\"exact solution\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03944572",
   "metadata": {},
   "source": [
    "## Example 2: Rate equation\n",
    "\n",
    "Solve the rate equation for $v(t)$ and plot the result.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dv}{dt} = g - cv^2\n",
    "\\end{align}\n",
    "with gravitational acceleration $g = 9.81$, drag coefficient $c = 1$, initial velocity $v_0 = 0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "db542247",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time\n",
    "\n",
    "t = np.linspace(0, 2, 200)\n",
    "\n",
    "# gravitational acceleration\n",
    "g = 9.81\n",
    "\n",
    "# drag coefficient\n",
    "c = 1\n",
    "\n",
    "# initial velocity\n",
    "v0 = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bd86499f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fn(v, t, g, c):\n",
    "    return g - c*v**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5ac28054",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_num = odeint(fn, v0, t, args=(g, c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4376be54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAGFCAYAAADElR/1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5AklEQVR4nO3df5zM9f7///vDrkXktzcVEv0i50itzunn8aODE9EPOtWnoh+UKJXS7zh6v/VD8iMcKTlKoeRXfaNfOKUoy5FflUoiRX7kN8vuPr5/7HA2LXbZ2efM7O16ucxlZl6v18zeZ55mZ++er3mNubsAAAAAAIh1xUIHAAAAAAAgLyiwAAAAAIC4QIEFAAAAAMQFCiwAAAAAIC5QYAEAAAAAcYECCwAAAACICxRYAACizMweMrMXj/I+apmZm1lyAeZqbGY/HsXtj/pxAQCQHwX2JggAAHLn7n1DZzhaZtZY0hh3r75vWSI8LgBAfGEGFgCAKCrIGVMAAIo6CiwAICGY2Uozu9fMFpnZFjMbb2YlI+s6mtnsA7Z3Mzs5cvlfZjbMzKaZ2XYz+8TMqpnZQDP71cy+MrOGOW57vJm9aWbrzex7M7szx7reZjbBzMaY2VZJHSPLxuTY5gIz+9TMNpvZajPrGFneysz+Y2ZbI8t75+Px329ma8xsm5l9bWbNIstLRB7HT5HTQDMrcZD72P+c5Hhe/tfMSkuaJun4yPOzPfIcHPi42pjZ0sjjmmVmdfMyPgAA5BUFFgCQSK6S1FLSSZL+KKljPm/7iKTKktIlzZG0IHJ9gqRnJcnMikl6S9IXkk6Q1EzSXWbWIsd9tY3cprykV3P+EDM7Udll8DlJVSSdKWlhZPUOSTdEbtdKUhczu+xwwc3sNEndJDVy92MltZC0MrL6YUl/jvycBpLOiTzOPHP3HZL+Jukndy8TOf10QIZTJY2VdFfkcb0j6S0zS8mx2dGMDwAAFFgAQEIZ7O4/ufsmZZfMM/Nx20nuPt/dd0uaJGm3u7/s7pmSxkvaNwPbSFIVd+/j7nvcfYWkFyRdneO+5rj7ZHfPcvddB/ycayV94O5j3X2vu29094WS5O6z3H1x5HaLlF0I/5KH7JmSSkiqZ2bF3X2lu38XWff/JPVx91/cfb2kf0i6Ph/PS179XdL/5+7vu/teSc9IKiXpvBzbHM34AABAgQUAJJS1OS7vlFQmH7ddl+Pyrlyu77uvE5W9K+3mfSdJD0mqmmP71Yf4OTUkfZfbCjP7k5nNjOyavEXSbcqeAT4kd/9W2TOfvSX9YmbjzOz4yOrjJf2QY/MfIssK2m9+jrtnKft5OCHHNkczPgAAUGABAEXCDknH7LtiZtWO4r5WS/re3cvnOB3r7pfk2MYPc/s6B1n3mqSpkmq4ezlJwyVZXkK5+2vufoGyC7ZLeiqy6qfIsn1qRpblZqdyPE+Scj5Ph3pMv/s5ZmbKLutrDhseAIA8osACAIqCLySdYWZnRg4c1Pso7utzSdsiB00qZWZJZlbfzBrl8favSrrYzK4ys2Qzq2RmZ0bWHStpk7vvNrNzlL278WGZ2Wlm1jRycKbdyp4xzoqsHivpETOrYmaVJT0macxB7mqhpGsjj6mlfrv78jpJlcys3EFu+7qkVmbWzMyKS+qh7M8Sf5qXxwAAQF5QYAEACc/dl0vqI+kDSd9Imn3oWxzyvjIltVb25ze/l7RB0ouSDlbsDrz9KkmXKLvgbVJ2aWwQWX27pD5mtk3ZRfP1PMYqIenJSJa1kv5H0oORdf8rKU3SIkmLlX1gqv89yP10l3SppM3K/uzs5By5v1J2GV4R2XX6N7shu/vXkq5T9sGpNkTu51J335PHxwAAwGGZ++H2CAIAAAAAIDxmYAEAAAAAcYECCwAAAACICxRYAAAAAEBcoMACAAAAAOICBRYAAAAAEBeSQwfIr8qVK3utWrVCxwAAAAAARMH8+fM3uHuV3NbFXYGtVauW0tLSQscAAAAAAESBmf1wsHXsQgwAAAAAiAsUWAAAAABAXKDAAgAAAADiAgUWAAAAABAXKLAAAAAAgLhAgQUAAAAAxAUKLAAAAAAgLlBgAQAAAABxgQILAAAAAIgLFFgAAAAAQFygwAIAAAAA4kJy6AAAAACxwt2VlZX1m/MDl+3bbt95zssHOz/abQrjZxR0jsM9z/G8PhYy8Bjz9hiQuwsuuEDHHnts6BhHhAILAEBgmZmZysjI0N69ewv0PCMjQ1lZWcrMzFRmZuZvLh/uekFte+DlfafciuGhSmNhrOOPYQBFxZIlS3TGGWeEjnFEKLAAgCIjMzNTu3bt0s6dO393np6evv+0Z8+eg14/2OW8bHewkhkLxcnMlJSUpGLFiikpKel3l/OzLrdtS5QosX+ZmcnMVKxYsVzPY32dme1/zvad57asoLcpjJ9R0DkO928untfHQgYeY94eA37vpJNOCh3hiFFgAQAxIzMzUzt27NC2bdu0ffv2g55v37491xJ6qPNdu3Zpz549R50xOTlZJUqUUIkSJZSSkrL/8oHXS5cu/Zt1xYsX339KTk6O+vm+U16LJ38EAgDiAQUWAFAg0tPTtXnz5v2nX3/99TfXc562bNmSazHduXNnnn9eSkqKSpUqpWOOOeZ355UqVVKNGjVUqlSpg26T87xUqVIqWbLkYYtpSkqKihXj+IcAAIRCgQUA/EZ6ero2bNigjRs3asOGDbmeciunu3fvPuT9Fi9eXOXLl1f58uVVrlw5HXvssapZs6aOPfZYlSlTJl/npUuXVnIyb2EAABQ1vPsDQIJzd23dulVr16496ClnOd2+fftB76t8+fKqVKmSKlasqPLly6tGjRr7S+nhTqVKlWI3VQAAcFQosAAQx3bu3Kkff/xRq1ev/s15znK6bt26XGdHixcvrmrVqqlq1aqqUqWK6tatq8qVK//mVKlSpf2XK1asqOLFiwd4lAAAANkosAAQozIzM/Xjjz9q5cqVWr169W9O+8rqxo0bf3e7ypUr67jjjlO1atV06qmnqlq1armeKlSowIwoAACIKxRYAAhox44dWrFihb777rv95/sur1y5Unv37v3N9hUrVlT16tVVo0YN/fnPf1aNGjX2n6pXr67q1aurZMmSgR4NAABAdFFgASDKdu3apW+++UZfffWVvvrqK33zzTf7y+q6det+s2358uVVp04dnXnmmbryyitVu3ZtnXTSSfsLaunSpQM9CgAAgPCiVmDNrKSkjySViPycCe7e64BtSkh6WdLZkjZK+ru7r4xWJgCIpi1btmjJkiVatmzZ/rL61Vdf6fvvv5e7S8r+wvUaNWqodu3aat26tWrXrq06deqoTp06ql27tipWrBj4UQAAAMSuaM7Apktq6u7bzay4pNlmNs3d5+bY5mZJv7r7yWZ2taSnJP09ipkA4Kjt3btXX3/9tRYvXqzFixdr0aJFWrx4sVatWrV/m1KlSum0007TOeecoxtuuEGnn366Tj/9dJ166qkqVapUwPQAAADxK2oF1rOnG/Z9F0PxyMkP2KytpN6RyxMkDTEz831TFQAQ2K5du7Ro0SKlpaVp/vz5WrBggZYtW7b/s6nJyck6/fTTdcEFF+iPf/yj/vCHP6hevXqqWbOmihUrFjg9AABAYonqZ2DNLEnSfEknSxrq7p8dsMkJklZLkrtnmNkWSZUkbTjgfjpL6ixJNWvWjGZkAEVYenq6Fi5cqPnz5+8vrEuXLlVmZqak7KP7nn322WrRosX+snr66acrJSUlcHIAAICiIaoF1t0zJZ1pZuUlTTKz+u6+5AjuZ4SkEZKUmprK7CyAArF27VrNmTNHn376qT799FPNnz9f6enpkv5bVlu3bq2zzz5bqampqlGjBl87AwAAEFChHIXY3Teb2UxJLSXlLLBrJNWQ9KOZJUsqp+yDOQFAgXJ3LV26VB999NH+0rpixQpJUkpKilJTU9WtWzede+65atSoEWUVAAAgBkXzKMRVJO2NlNdSkv6q7IM05TRVUgdJcyS1kzSDz78CKAjuruXLl2vmzJmaMWOGZs2apfXr10uSqlWrpvPPP1+33367zjvvPJ111lkqUaJE4MQAAAA4nGjOwB4naXTkc7DFJL3u7m+bWR9Jae4+VdJISa+Y2beSNkm6Oop5ACS4H374Qe+//75mzpypmTNn6ueff5YkVa9eXS1btlSTJk3UuHFj1apVi9lVAACAOBTNoxAvktQwl+WP5bi8W1L7aGUAkNjS09P18ccfa9q0aZo2bZq+/PJLSVLVqlXVpEkTNW3aVE2aNFGdOnUorAAAAAmgUD4DCwAF5YcffthfWD/88EPt2LFDKSkp+stf/qJOnTqpRYsWqlu3LoUVAAAgAVFgAcQ0d9fixYv15ptvauLEiVqyJPs4cLVq1dINN9ygv/3tb2ratKlKly4dOCkAAACijQILIOZkZWVp3rx5+0vrd999JzPThRdeqGeeeUaXXHKJTj/9dGZZAQAAihgKLICYkJWVpdmzZ2vChAmaOHGi1qxZo+TkZDVr1kw9e/ZU27ZtVbVq1dAxAQAAEBAFFkBQy5Yt05gxY/Tqq69q1apVKlmypFq0aKEnnnhCrVu3VoUKFUJHBAAAQIygwAIodGvXrtXYsWM1ZswYLViwQElJSWrevLmeeOIJtWnTRmXKlAkdEQAAADGIAgugUOzevVsTJ07UK6+8ovfee09ZWVk6++yzNWDAAF1zzTXsHgwAAIDDosACiKrly5drxIgRGjVqlDZt2qSaNWvq/vvv13XXXad69eqFjgcAAIA4QoEFUOD27NmjKVOm6Pnnn9eHH36o5ORkXXbZZbr11lvVtGlTFStWLHREAAAAxCEKLIACs3LlSr3wwgsaOXKk1q1bpxNPPFH/93//p5tuuknVqlULHQ8AAABxjgIL4KjNnTtX/fr106RJk2RmatWqlW677Ta1aNFCSUlJoeMBAAAgQVBgARyRzMxMvfXWW3rmmWf0ySefqEKFCnrggQd02223qWbNmqHjAQAAIAFRYAHky65du/Tyyy+rf//++uabb1SrVi0NHjxYN954I19/AwAAgKiiwALIk02bNum5557TkCFDtGHDBqWmpmr8+PG64oorlJzMrxIAAABEH391AjikzZs369lnn9XAgQO1bds2tW7dWvfee68uuugimVnoeAAAAChCKLAAcrV161YNGjRI/fv315YtW9SuXTv16tVL9evXDx0NAAAARRQFFsBvbNu2TUOGDFG/fv3066+/6rLLLlPv3r3VoEGD0NEAAABQxFFgAUiSdu7cqaFDh+rpp5/Whg0b1Lp1a/Xu3Vtnn3126GgAAACAJKlY6AAAwsrKytLLL7+sU089VT179lRqaqo+++wzvfXWW5RXAAAAxBQKLFCEzZ49W3/605/UoUMHHX/88fr44481bdo0nXPOOaGjAQAAAL9DgQWKoO+//15XXXWVLrzwQv3888965ZVXNHfuXF1wwQWhowEAAAAHxWdggSJk69at6tu3rwYMGKDk5GT17t1b9957r0qXLh06GgAAAHBYFFigCMjKytKoUaP00EMP6ZdfftH111+vvn37qnr16qGjAQAAAHlGgQUS3OLFi3Xbbbfp008/1Xnnnae3335bjRo1Ch0LAAAAyDc+AwskqO3bt+u+++5Tw4YNtXz5co0aNUqzZ8+mvAIAACBuMQMLJKB3331XnTt31qpVq3TLLbfoySefVKVKlULHAgAAAI4KM7BAAtm4caM6dOigli1bqnTp0po9e7ZeeOEFyisAAAASAgUWSBATJkxQvXr19Nprr+nRRx/Vf/7zH51//vmhYwEAAAAFhl2IgTi3adMmde3aVePGjdPZZ5+t9957Tw0aNAgdCwAAAChwzMACcWzatGmqX7++JkyYoMcff1xz586lvAIAACBhUWCBOLRjxw7ddtttuuSSS1SpUiV9/vnneuSRR5SczE4VAAAASFwUWCDOLFy4UKmpqRoxYoTuvfdezZs3Tw0bNgwdCwAAAIg6CiwQJ7KysjRgwAD96U9/0tatW/XBBx+oX79+KlmyZOhoAAAAQKFgf0MgDmzYsEE33HCDpk2bpjZt2mjkyJGqXLly6FgAAABAoaLAAjHu008/1VVXXaUNGzZo6NCh6tKli8wsdCwAAACg0LELMRCj3F39+/fXX/7yF5UsWVJz5szR7bffTnkFAABAkcUMLBCDtmzZoo4dO2ry5Mm64oor9NJLL6lcuXKhYwEAAABBUWCBGPPll1/qsssu04oVKzRgwAB1796dWVcAAABAUdyF2MxqmNlMM1tmZkvNrHsu2zQ2sy1mtjByeixaeYB4MHnyZJ1zzjnavHmzPvzwQ911112UVwAAACAimp+BzZDUw93rSfqzpK5mVi+X7T529zMjpz5RzAPErKysLD322GO6/PLLVbduXaWlpemiiy4KHQsAAACIKVHbhdjdf5b0c+TyNjP7UtIJkpZF62cC8WjHjh26/vrrNWnSJN14440aNmwY3+0KAAAA5KJQjkJsZrUkNZT0WS6rzzWzL8xsmpmdcZDbdzazNDNLW79+fTSjAoVqzZo1uvDCCzVlyhQNGDBAI0eOpLwCAAAABxH1gziZWRlJb0q6y923HrB6gaQT3X27mV0iabKkUw68D3cfIWmEJKWmpnp0EwOFIy0tTW3atNH27ds1depUtWrVKnQkAAAAIKZFdQbWzIoru7y+6u4TD1zv7lvdfXvk8juSiptZ5WhmAmLBlClTdNFFFyklJUWffvop5RUAAADIg2gehdgkjZT0pbs/e5BtqkW2k5mdE8mzMVqZgFgwdOhQXXHFFfrjH/+ozz//XPXr1w8dCQAAAIgL0dyF+HxJ10tabGYLI8seklRTktx9uKR2krqYWYakXZKudnd2EUZCysrK0kMPPaSnnnpKbdq00dixY3XMMceEjgUAAADEjWgehXi2pEN+gaW7D5E0JFoZgFiRnp6um266Sa+99pq6dOmiwYMHKzk56h9BBwAAABIKf0EDUbZ9+3ZdccUVev/99/XEE0/o/vvvV2TPeQAAAAD5QIEFomjTpk1q1aqVPv/8c40aNUodO3YMHQkAAACIWxRYIEp++ukntWjRQsuXL9ebb76pyy67LHQkAAAAIK5RYIEo+O677/TXv/5V69ev17Rp09S0adPQkQAAAIC4R4EFCtjXX3+tpk2bavfu3ZoxY4YaNWoUOhIAAACQECiwQAFaunSpmjVrJnfXrFmz9Ic//CF0JAAAACBhFAsdAEgUX3zxhRo3bqxixYpRXgEAAIAooMACBWDBggVq0qSJSpYsqX//+9+qW7du6EgAAABAwmEXYuAoLVy4UBdffLHKli2rmTNn6qSTTgodCQAAAEhIzMACR2HJkiW6+OKLVaZMGcorAAAAEGUUWOAIffnll2rWrJlKlCihGTNmUF4BAACAKKPAAkdg+fLlatq0qYoVK6YZM2bo5JNPDh0JAAAASHh8BhbIp1WrVuniiy9WZmamZs2apdNOOy10JAAAAKBIoMAC+bBu3TpdfPHF2rp1q/7973+rXr16oSMBAAAARQYFFsijzZs3q0WLFlqzZo3ee+89NWjQIHQkAAAAoEihwAJ5sHPnTrVu3VrLli3TW2+9pfPPPz90JAAAAKDIocACh7F37161a9dOc+bM0fjx49WiRYvQkQAAAIAiiQILHIK7q1OnTpo2bZpGjBihdu3ahY4EAAAAFFl8jQ5wCI888ohGjx6t3r17q1OnTqHjAAAAAEUaBRY4iGHDhqlv377q1KmTHnvssdBxAAAAgCKPAgvkYtKkSerWrZsuvfRSDRs2TGYWOhIAAABQ5FFggQN8/vnnuvbaa3XOOedo3LhxSk7mo+IAAABALKDAAjn88MMPatOmjY477ji99dZbOuaYY0JHAgAAABDB1BIQsWXLFrVq1Uq7d+/WzJkzVaVKldCRAAAAAORAgQWU/V2vV111lb7++mtNnz5ddevWDR0JAAAAwAEosCjy3F133nmn3nvvPb344otq1qxZ6EgAAAAAcsFnYFHk/fOf/9Tw4cPVs2dP3XzzzaHjAAAAADgICiyKtJkzZ6p79+5q3bq1+vbtGzoOAAAAgEOgwKLIWrFihdq3b69TTjlFr776qpKSkkJHAgAAAHAIFFgUSdu2bVPbtm2VlZWlqVOnqmzZsqEjAQAAADgMDuKEIicrK0sdOnTQsmXLNH36dJ188smhIwEAAADIAwosipwnn3xSkyZN0rPPPqu//vWvoeMAAAAAyCN2IUaR8t577+mRRx7Rtddeq7vuuit0HAAAAAD5QIFFkbFy5Updc801OuOMMzRixAiZWehIAAAAAPKBAosiYffu3WrXrp0yMjI0ceJElS5dOnQkAAAAAPnEZ2BRJHTr1k3z58/X1KlTdcopp4SOAwAAAOAIMAOLhPevf/1LI0eO1MMPP6xLL700dBwAAAAAR4gCi4S2ZMkS3X777WrSpIn+8Y9/hI4DAAAA4ChErcCaWQ0zm2lmy8xsqZl1z2UbM7PBZvatmS0ys7OilQdFz/bt29W+fXuVLVtWr732mpKSkkJHAgAAAHAUovkZ2AxJPdx9gZkdK2m+mb3v7stybPM3SadETn+S9M/IOXBU3F1dunTR8uXL9f7776tatWqhIwEAAAA4SlGbgXX3n919QeTyNklfSjrhgM3aSnrZs82VVN7MjotWJhQdL730ksaMGaNevXqpadOmoeMAAAAAKACF8hlYM6slqaGkzw5YdYKk1Tmu/6jfl1wgX5YuXapu3brp4osv1sMPPxw6DgAAAIACEvUCa2ZlJL0p6S5333qE99HZzNLMLG39+vUFGxAJZdeuXbr66qtVtmxZjRkzhs+9AgAAAAkkqgXWzIoru7y+6u4Tc9lkjaQaOa5Xjyz7DXcf4e6p7p5apUqV6IRFQrjvvvu0ZMkSjR49WlWrVg0dBwAAAEABiuZRiE3SSElfuvuzB9lsqqQbIkcj/rOkLe7+c7QyIbFNnTpVQ4cO1T333KOWLVuGjgMAAACggEXzKMTnS7pe0mIzWxhZ9pCkmpLk7sMlvSPpEknfStop6cYo5kECW7NmjW666SY1bNhQffv2DR0HAAAAQBRErcC6+2xJdphtXFLXaGVA0ZCVlaUbbrhBu3bt0tixY1WiRInQkQAAAABEQTRnYIFC8eyzz2rGjBkaOXKkTjvttNBxAAAAAERJoXyNDhAtixYt0sMPP6zLL79cN97IHugAAABAIqPAIm6lp6fruuuuU4UKFfT8888r+7hhAAAAABIVuxAjbj366KNavHix3n77bfH1SgAAAEDiYwYWcenf//63nnnmGd16661q1apV6DgAAAAACgEFFnFn69at6tChg+rUqaNnnnkmdBwAAAAAhYRdiBF3evToodWrV2v27NkqU6ZM6DgAAAAACgkzsIgr06dP14svvqj77rtP5557bug4AAAAAAqRuXvoDPmSmprqaWlpoWMggM2bN6t+/foqV66c5s+fr5IlS4aOBAAAAKCAmdl8d0/NbR27ECNu3HPPPVq7dq0mTZpEeQUAAACKIHYhRlx45513NGrUKPXs2VONGjUKHQcAAABAAOxCjJi3efNmnXHGGapQoYLmz5+vEiVKhI4EAAAAIErYhRhx7d5779W6des0ZcoUyisAAABQhLELMWLahx9+qJEjR6pHjx5KTc31P2EAAAAAFBHsQoyYtXPnTv3hD39QsWLFtGjRIpUqVSp0JAAAAABRxi7EiEuPPfaYVqxYoZkzZ1JeAQAAALALMWLTvHnzNGDAAHXu3FmNGzcOHQcAAABADKDAIubs2bNHN998s6pVq6ann346dBwAAAAAMYJdiBFznnnmGS1evFiTJ09WuXLlQscBAAAAECOYgUVM+e677/T444/ryiuvVNu2bUPHAQAAABBDKLCIGe6u22+/XcWLF9egQYNCxwEAAAAQY9iFGDFj/Pjxeu+99zR48GCdcMIJoeMAAAAAiDF8DyxiwubNm3X66aerRo0amjt3rpKSkkJHAgAAABAA3wOLmPfggw9q/fr1mjZtGuUVAAAAQK74DCyCmzNnjoYPH67u3burYcOGoeMAAAAAiFEUWASVkZGh22+/XdWrV1efPn1CxwEAAAAQw9iFGEENHz5cCxcu1BtvvKEyZcqEjgMAAAAghjEDi2DWrVunRx55RM2bN9eVV14ZOg4AAACAGEeBRTA9e/bUzp079dxzz8nMQscBAAAAEOMosAji448/1ssvv6z77rtPp556aug4AAAAAOIABRaFLiMjQ127dlXNmjX10EMPhY4DAAAAIE5wECcUuqFDh2rx4sWaOHGiSpcuHToOAAAAgDjBDCwK1fr169WrVy81b95cl112Weg4AAAAAOIIBRaF6rHHHtP27ds1YMAADtwEAAAAIF8OW2DN7BMza1IYYZDYFi1apBEjRuj2229XvXr1QscBAAAAEGfyMgN7q6RuZvahmZ0b7UBITO6uu+++W+XLl1fv3r1DxwEAAAAQhw57ECd3XyLpSjM7S1KfyG6fD7v7F9EOh8QxZcoUzZgxQ88995wqVqwYOg4AAACAOGTunrcNzcpKqivpJkk3u3uQIxinpqZ6WlpaiB+NI5Senq4zzjhDJUqU0BdffKHkZA5+DQAAACB3Zjbf3VNzW3fYJmFmMyWdImmXpGWRU8eCDIjENmjQIH333Xd69913Ka8AAAAAjlhe2kQPSV+6+6783LGZvSSptaRf3L1+LusbS5oi6fvIoonu3ic/PwOxb+3atXr88cfVunVrNW/ePHQcAAAAAHEsL5+BXXCE9/0vSUMkvXyIbT5299ZHeP+IAw8//LDS09PVv3//0FEAAAAAxLmofQ+su38kaVO07h+xb8GCBRo1apTuvPNOnXrqqaHjAAAAAIhzUSuweXSumX1hZtPM7IyDbWRmnc0szczS1q9fX5j5cITcXd27d1flypX16KOPho4DAAAAIAGEPKLOAkknuvt2M7tE0mRlHyzqd9x9hKQRUvZRiAstIY7YG2+8odmzZ+v5559XuXLlQscBAAAAkACCzcC6+1Z33x65/I6k4mZWOVQeFJxdu3bpvvvuU4MGDXTzzTeHjgMAAAAgQQSbgTWzapLWubub2TnKLtMbQ+VBwenfv79WrVql0aNHKykpKXQcAAAAAAkiagXWzMZKaiypspn9KKmXpOKS5O7DJbWT1MXMMpT9HbNXuzu7B8e5X375RU899ZQuv/xyNW7cOHQcAAAAAAkkagXW3a85zPohyv6aHSSQPn36aNeuXXriiSdCRwEAAACQYEIfhRgJ5Ntvv9Xzzz+vW265RaeddlroOAAAAAASDAUWBebhhx9WSkqKevXqFToKAAAAgAREgUWBmDdvnl5//XX16NFDxx13XOg4AAAAABIQBRZHzd3Vs2dPValSRffee2/oOAAAAAASVLCv0UHimD59umbNmqXBgwerbNmyoeMAAAAASFAWb99ck5qa6mlpaaFjICIzM1MNGzbUzp07tWzZMqWkpISOBAAAACCOmdl8d0/NbR0zsDgqr776qhYvXqxx48ZRXgEAAABEFTOwOGK7d+/WaaedpipVqujzzz9XsWJ8pBoAAADA0WEGFlExZMgQrVq1SqNGjaK8AgAAAIg6WgeOyK+//qq+ffuqZcuWatq0aeg4AAAAAIoACiyOyJNPPqnNmzfrySefDB0FAAAAQBFBgUW+rV69WoMGDdJ1112nBg0ahI4DAAAAoIigwCLfevXqJXfX448/HjoKAAAAgCKEAot8+eqrrzR69Gh17dpVJ554Yug4AAAAAIoQCizypXfv3ipVqpQefPDB0FEAAAAAFDEUWOTZokWLNH78eHXv3l1VqlQJHQcAAABAEUOBRZ716tVL5cqV07333hs6CgAAAIAiiAKLPElLS9PkyZPVo0cPVahQIXQcAAAAAEUQBRZ58uijj6pixYrq3r176CgAAAAAiigKLA7rk08+0fTp03X//ferbNmyoeMAAAAAKKIosDisRx99VFWrVlXXrl1DRwEAAABQhCWHDoDYNmPGDM2cOVODBg1S6dKlQ8cBAAAAUIQxA4uDcnc98sgjql69ujp37hw6DgAAAIAijhlYHNT06dM1Z84cDR8+XCVLlgwdBwAAAEARxwwscrVv9vWkk07SjTfeGDoOAAAAADADi9xNnjxZCxYs0L/+9S+lpKSEjgMAAAAAMncPnSFfUlNTPS0tLXSMhJaVlaUGDRpoz549Wrp0qZKT+X8OAAAAAIXDzOa7e2pu62gm+J3XX39dS5Ys0dixYymvAAAAAGIGM7D4jYyMDJ1xxhlKSUnRF198oWLF+Jg0AAAAgMLDDCzybPz48Vq+fLnefPNNyisAAACAmMIMLPbLzMxU/fr1lZyczOwrAAAAgCCYgUWeTJgwQV999ZXGjx9PeQUAAAAQc5iBhaT/Hnk4MzNTixcvVlJSUuhIAAAAAIogZmBxWJMnT9aSJUv06quvUl4BAAAAxCRmYCF3V8OGDbVr1y4tW7aMAgsAAAAgGGZgcUhvvfWWvvjiC40ePZryCgAAACBmcaSeIs7d1adPH9WuXVvXXntt6DgAAAAAcFDMwBZx06ZN0/z58/Xiiy8qOZl/DgAAAABiV9RmYM3sJTP7xcyWHGS9mdlgM/vWzBaZ2VnRyoLc7Zt9PfHEE3X99deHjgMAAAAAhxTNXYj/JanlIdb/TdIpkVNnSf+MYhbk4oMPPtBnn32mBx98UCkpKaHjAAAAAMAhRa3AuvtHkjYdYpO2kl72bHMllTez46KVB7/l7vrHP/6h6tWrq2PHjqHjAAAAAMBhhTyI0wmSVue4/mNk2e+YWWczSzOztPXr1xdKuEQ3a9YsffLJJ3rggQdUokSJ0HEAAAAA4LDi4ijE7j7C3VPdPbVKlSqh4ySEPn366LjjjtPNN98cOgoAAAAA5EnIw86ukVQjx/XqkWWIso8++kizZs3SwIEDVbJkydBxAAAAACBPQs7ATpV0Q+RoxH+WtMXdfw6Yp8h4/PHH9T//8z/q1KlT6CgAAAAAkGdRm4E1s7GSGkuqbGY/SuolqbgkuftwSe9IukTSt5J2SroxWlnwX3PnztUHH3ygp59+Wsccc0zoOAAAAACQZ1ErsO5+zWHWu6Su0fr5yN0TTzyhihUrqkuXLqGjAAAAAEC+xMVBnFAwlixZoqlTp+qOO+5QmTJlQscBAAAAgHyhwBYhTz31lEqXLq077rgjdBQAAAAAyDcKbBHx/fffa+zYsbr11ltVqVKl0HEAAAAAIN8osEXEM888o2LFiumee+4JHQUAAAAAjggFtghYu3atRo4cqQ4dOuiEE04IHQcAAAAAjggFtggYOHCg9u7dq549e4aOAgAAAABHjAKb4DZv3qxhw4apXbt2OuWUU0LHAQAAAIAjRoFNcMOGDdO2bdv0wAMPhI4CAAAAAEeFApvAdu7cqYEDB6ply5Zq2LBh6DgAAAAAcFQosAnspZde0vr16/Xggw+GjgIAAAAAR40Cm6D27t2rfv366bzzztOFF14YOg4AAAAAHLXk0AEQHWPHjtWqVas0bNgwmVnoOAAAAABw1MzdQ2fIl9TUVE9LSwsdI6ZlZWWpfv36Kl68uBYuXEiBBQAAABA3zGy+u6fmto4Z2AQ0depUffnll3rttdcorwAAAAASBp+BTTDurr59+6p27dpq37596DgAAAAAUGCYgU0wM2fO1Lx58zR8+HAlJzO8AAAAABIHM7AJ5umnn1bVqlXVoUOH0FEAAAAAoEBRYBPIokWL9O677+rOO+9UyZIlQ8cBAAAAgAJFgU0g/fr1U+nSpdWlS5fQUQAAAACgwFFgE8SqVas0btw4derUSRUqVAgdBwAAAAAKHAU2QQwcOFDurrvvvjt0FAAAAACICgpsAti8ebNeeOEFXX311apZs2boOAAAAAAQFRTYBDB8+HBt375d9913X+goAAAAABA1FNg4l56erkGDBql58+Zq0KBB6DgAAAAAEDXJoQPg6IwZM0Zr167VK6+8EjoKAAAAAEQVM7BxLCsrS/369VPDhg3VrFmz0HEAAAAAIKqYgY1jb7/9tr7++mu99tprMrPQcQAAAAAgqpiBjWP9+vXTiSeeqPbt24eOAgAAAABRxwxsnJozZ45mz56tQYMGKTmZYQQAAACQ+JiBjVP9+vVThQoVdNNNN4WOAgAAAACFggIbh5YvX67Jkyera9euKlOmTOg4AAAAAFAoKLBxqH///kpJSVG3bt1CRwEAAACAQkOBjTPr1q3T6NGj1aFDB1WtWjV0HAAAAAAoNBTYODNkyBDt2bNHPXr0CB0FAAAAAAoVBTaO7Ny5U//85z/Vtm1bnXrqqaHjAAAAAEChosDGkdGjR2vjxo265557QkcBAAAAgEJHgY0TWVlZGjBggBo1aqQLLrggdBwAAAAAKHTJoQMgb95++2198803GjdunMwsdBwAAAAAKHRRnYE1s5Zm9rWZfWtmD+SyvqOZrTezhZHTLdHME8/69++vmjVr6sorrwwdBQAAAACCiNoMrJklSRoq6a+SfpQ0z8ymuvuyAzYd7+58oekhpKWl6aOPPlL//v2VnMykOQAAAICiKZozsOdI+tbdV7j7HknjJLWN4s9LWM8++6zKli2rW25hghoAAABA0RXNAnuCpNU5rv8YWXagK81skZlNMLMaUcwTl1atWqXXX39dnTp1UtmyZUPHAQAAAIBgQh+F+C1Jtdz9j5LelzQ6t43MrLOZpZlZ2vr16ws1YGiDBw+WJN15552BkwAAAABAWNEssGsk5ZxRrR5Ztp+7b3T39MjVFyWdndsdufsId09199QqVapEJWws2rp1q1544QW1b99eNWvWDB0HAAAAAIKKZoGdJ+kUMzvJzFIkXS1pas4NzOy4HFfbSPoyinnizsiRI7V161b16NEjdBQAAAAACC5qh7R19wwz6ybpXUlJkl5y96Vm1kdSmrtPlXSnmbWRlCFpk6SO0coTbzIyMjRo0CBdeOGFSk1NDR0HAAAAAIKL6neyuPs7kt45YNljOS4/KOnBaGaIV2+++aZ++OEHDRo0KHQUAAAAAIgJoQ/ihFy4u/r3769TTjlFl156aeg4AAAAABATojoDiyPzySefaN68eRo2bJiKFeP/GAAAAABAYgY2JvXv318VK1ZUhw4dQkcBAAAAgJhBgY0x33zzjaZMmaIuXbromGOOCR0HAAAAAGIGBTbGDBo0SMWLF1e3bt1CRwEAAACAmEKBjSG//vqrRo0apWuuuUbVqlULHQcAAAAAYgoFNoa8+OKL2rlzp+6+++7QUQAAAAAg5lBgY0RGRoaee+45NWnSRA0aNAgdBwAAAABiDl+jEyMmTpyo1atXa8iQIaGjAAAAAEBMYgY2RgwcOFB16tRRq1atQkcBAAAAgJjEDGwM+OyzzzRnzhwNHjxYSUlJoeMAAAAAQExiBjYGDBw4UGXLllXHjh1DRwEAAACAmEWBDWz16tV644031KlTJx177LGh4wAAAABAzKLABjZ06FC5u7p16xY6CgAAAADENApsQDt27NCIESN0xRVXqFatWqHjAAAAAEBMo8AG9PLLL+vXX3/VXXfdFToKAAAAAMQ8CmwgWVlZGjhwoBo1aqTzzjsvdBwAAAAAiHl8jU4g06dP1/Lly/Xqq6/KzELHAQAAAICYxwxsIAMGDNDxxx+v9u3bh44CAAAAAHGBAhvAkiVL9MEHH6hbt24qXrx46DgAAAAAEBcosAEMHDhQpUqVUufOnUNHAQAAAIC4QYEtZOvXr9eYMWN0ww03qFKlSqHjAAAAAEDcoMAWsueff17p6enq3r176CgAAAAAEFcosIUoPT1dQ4cOVcuWLVW3bt3QcQAAAAAgrlBgC9Hrr7+utWvX6q677godBQAAAADiDgW2kLi7BgwYoLp166p58+ah4wAAAABA3EkOHaCo+Pjjj/Wf//xHzz//vMwsdBwAAAAAiDvMwBaSAQMGqFKlSrr++utDRwEAAACAuESBLQQrVqzQlClTdOutt6pUqVKh4wAAAABAXKLAFoLBgwcrKSlJXbt2DR0FAAAAAOIWBTbKtm7dqpdeekl///vfdfzxx4eOAwAAAABxiwIbZSNHjtS2bdv46hwAAAAAOEoU2CjKzMzU4MGDdcEFFyg1NTV0HAAAAACIaxTYKJo6dapWrlzJ7CsAAAAAFAAKbBQNGDBAtWrV0mWXXRY6CgAAAADEPQpslMyfP18ff/yx7rjjDiUlJYWOAwAAAABxjwIbJQMHDlSZMmV08803h44CAAAAAAmBAhsFP/30k8aPH6+bbrpJ5cqVCx0HAAAAABICBTYKhg0bpoyMDN15552howAAAABAwohqgTWzlmb2tZl9a2YP5LK+hJmNj6z/zMxqRTNPYdi1a5eGDx+uNm3aqE6dOqHjAAAAAEDCiFqBNbMkSUMl/U1SPUnXmFm9Aza7WdKv7n6ypAGSnopWnsIyZswYbdy4UXfffXfoKAAAAACQUKI5A3uOpG/dfYW775E0TlLbA7ZpK2l05PIESc3MzKKYKarcXQMHDtSZZ56piy66KHQcAAAAAEgoyVG87xMkrc5x/UdJfzrYNu6eYWZbJFWStCHnRmbWWVJnSapZs2a08h61rVu3ql69emrbtq3iuIcDAAAAQEyKZoEtMO4+QtIISUpNTfXAcQ6qXLlyeuONN0LHAAAAAICEFM1diNdIqpHjevXIsly3MbNkSeUkbYxiJgAAAABAnIpmgZ0n6RQzO8nMUiRdLWnqAdtMldQhcrmdpBnuHrMzrAAAAACAcKK2C3HkM63dJL0rKUnSS+6+1Mz6SEpz96mSRkp6xcy+lbRJ2SUXAAAAAIDfiepnYN39HUnvHLDssRyXd0tqH80MAAAAAIDEEM1diAEAAAAAKDAUWAAAAABAXKDAAgAAAADiAgUWAAAAABAXKLAAAAAAgLhAgQUAAAAAxAUKLAAAAAAgLlBgAQAAAABxgQILAAAAAIgL5u6hM+SLma2X9EPoHIdRWdKG0CHwG4xJbGJcYg9jEpsYl9jDmMQmxiX2MCaxJx7G5ER3r5LbirgrsPHAzNLcPTV0DvwXYxKbGJfYw5jEJsYl9jAmsYlxiT2MSeyJ9zFhF2IAAAAAQFygwAIAAAAA4gIFNjpGhA6A32FMYhPjEnsYk9jEuMQexiQ2MS6xhzGJPXE9JnwGFgAAAAAQF5iBBQAAAADEBQpsPphZSzP72sy+NbMHcllfwszGR9Z/Zma1cqx7MLL8azNrUajBE1wexuUeM1tmZovM7EMzOzHHukwzWxg5TS3c5IkrD2PS0czW53jub8mxroOZfRM5dSjc5IktD+MyIMeYLDezzTnW8VqJAjN7ycx+MbMlB1lvZjY4MmaLzOysHOt4rURBHsbk/0XGYrGZfWpmDXKsWxlZvtDM0govdeLLw7g0NrMtOX5PPZZj3SF/9+HI5GFM7ssxHksi7yMVI+t4rUSBmdUws5mRv3uXmln3XLaJ//cVd+eUh5OkJEnfSaotKUXSF5LqHbDN7ZKGRy5fLWl85HK9yPYlJJ0UuZ+k0I8pEU55HJcmko6JXO6yb1wi17eHfgyJdsrjmHSUNCSX21aUtCJyXiFyuULox5QIp7yMywHb3yHppRzXea1EZ1wuknSWpCUHWX+JpGmSTNKfJX0WWc5rJdyYnLfvuZb0t31jErm+UlLl0I8hEU95GJfGkt7OZXm+fvdxKrgxOWDbSyXNyHGd10p0xuQ4SWdFLh8raXkuf4PF/fsKM7B5d46kb919hbvvkTROUtsDtmkraXTk8gRJzczMIsvHuXu6u38v6dvI/eHoHXZc3H2mu++MXJ0rqXohZyxq8vJaOZgWkt53903u/quk9yW1jFLOoia/43KNpLGFkqwIc/ePJG06xCZtJb3s2eZKKm9mx4nXStQcbkzc/dPIcy7xnlJo8vBaOZijeU/CIeRzTHhPKQTu/rO7L4hc3ibpS0knHLBZ3L+vUGDz7gRJq3Nc/1G//wexfxt3z5C0RVKlPN4WRya/z+3Nyv5fp31Kmlmamc01s8uikK8oyuuYXBnZdWWCmdXI522Rf3l+biO72Z8kaUaOxbxWwjjYuPFaiQ0Hvqe4pPfMbL6ZdQ6UqSg718y+MLNpZnZGZBmvlcDM7BhlF6E3cyzmtRJllv1RxoaSPjtgVdy/rySHDgAUFjO7TlKqpL/kWHyiu68xs9qSZpjZYnf/LkzCIuUtSWPdPd3MblX2ngtNA2fCf10taYK7Z+ZYxmsFyMHMmii7wF6QY/EFkdfJ/0h638y+isxSIfoWKPv31HYzu0TSZEmnhI2EiEslfeLuOWdrea1EkZmVUfZ/GNzl7ltD5ylozMDm3RpJNXJcrx5Zlus2ZpYsqZykjXm8LY5Mnp5bM7tY0sOS2rh7+r7l7r4mcr5C0ixl/08Vjs5hx8TdN+YYhxclnZ3X2+KI5ee5vVoH7OrFayWYg40br5WAzOyPyv7d1dbdN+5bnuN18oukSeLjQoXG3be6+/bI5XckFTezyuK1EgsO9Z7Ca6WAmVlxZZfXV919Yi6bxP37CgU27+ZJOsXMTjKzFGW/GA88EudUSfuO2NVO2R9W98jyqy37KMUnKft/BD8vpNyJ7rDjYmYNJT2v7PL6S47lFcysRORyZUnnS1pWaMkTV17G5LgcV9so+zMakvSupOaRsakgqXlkGY5eXn6HycxOV/bBG+bkWMZrJZypkm6IHDXyz5K2uPvP4rUSjJnVlDRR0vXuvjzH8tJmduy+y8oek1yPzoqCZ2bVIscdkZmdo+y/cTcqj7/7EB1mVk7Ze75NybGM10qURF4DIyV96e7PHmSzuH9fYRfiPHL3DDPrpuyBTFL20TmXmlkfSWnuPlXZ/2BeMbNvlf2h9qsjt11qZq8r+w++DEldD9g1D0coj+PST1IZSW9E3ttWuXsbSXUlPW9mWcp+o3vS3fmj/CjlcUzuNLM2yn49bFL2UYnl7pvM7HFl/8EhSX0O2OUIRyiP4yJl/94aF/nPt314rUSJmY1V9tFTK5vZj5J6SSouSe4+XNI7yj5i5LeSdkq6MbKO10qU5GFMHlP28S2GRd5TMtw9VVJVSZMiy5Ilvebu0wv9ASSoPIxLO0ldzCxD0i5JV0d+j+X6uy/AQ0g4eRgTSbpc0nvuviPHTXmtRM/5kq6XtNjMFkaWPSSpppQ47yv2279RAAAAAACITexCDAAAAACICxRYAAAAAEBcoMACAAAAAOICBRYAAAAAEBcosAAAAACAuECBBQAAAADEBQosAAAAACAuUGABAIhRZlbdzP4eOgcAALGCAgsAQOxqJums0CEAAIgV5u6hMwAAgAOY2QWSpkjaLGmbpCvcfUXQUAAABEaBBQAgRpnZdEn3uvuS0FkAAIgF7EIMAEDsOk3SV6FDAAAQKyiwAADEIDOrLGmLu2eEzgIAQKygwAIAEJtqSfopdAgAAGIJBRYAgNj0laTKZrbEzM4LHQYAgFjAQZwAAAAAAHGBGVgAAAAAQFygwAIAAAAA4gIFFgAAAAAQFyiwAAAAAIC4QIEFAAAAAMQFCiwAAAAAIC5QYAEAAAAAcYECCwAAAACIC/8/dUWYNOi/KEIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 1, figsize=(16, 6))\n",
    "axs.plot(t, v_num, label=\"numerical solution\", c=\"black\")\n",
    "axs.set_title(\"numerical solution\")\n",
    "axs.set_ylabel(\"$v$\")\n",
    "axs.set_xlabel(\"$t$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39b51865",
   "metadata": {},
   "source": [
    "## Example 3: Multiple equations\n",
    "\n",
    "We have multiple rate equations in multiple unknowns, we solve the problem in the same way, using vector notation:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dx}{dt} &= v,\\\\\n",
    "\\frac{dv}{dt} &= g - cv^2.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9d3fbffe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time\n",
    "\n",
    "t = np.linspace(0, 2, 200)\n",
    "\n",
    "# gravitational acceleration\n",
    "g = 9.81\n",
    "\n",
    "# drag coefficient\n",
    "c = 1\n",
    "\n",
    "# initial velocity and position\n",
    "X0 = np.array([0.0, 0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fdeae5f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fn(X, t, g, c):\n",
    "    x = X[0]\n",
    "    v = X[1]\n",
    "    dxdt = v\n",
    "    dvdt = g - c*v**2\n",
    "    return np.array([dxdt, dvdt])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5a1f51f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_num = odeint(fn, X0, t, args=(g, c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "86873fcb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6YAAAGFCAYAAADuJe5aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPj0lEQVR4nO3deZzNdf//8efLGKRUltFCttYrbaSubJFkjH2P0qJQqVyK6yqtV/Vt1VXpQmXpKiSEbIksw5iFZkaWRIuikspeGBMz798fc/SbhAZzzvucM4/77XZuzvKZz+fpOOZ9nufz+byPOecEAAAAAIAvxXwHAAAAAAAUbRRTAAAAAIBXFFMAAAAAgFcUUwAAAACAVxRTAAAAAIBXFFMAAAAAgFcUUyBCmNnrZvboER5/yMxGhjITAADRzsxuNbPk41xHQzP7vLAyAdHI+B5TIPKYWWNJY51zlT1HAQAgqpnZrZJ6OucaFOI61wfWOa+w1glEOvaYAgAAAAC8opgCQWJm681soJl9Zmbbzex/ZlYq8FgvM/vKzLaZ2XQzOzNwv5nZy2b2s5n9YmarzOyiwGNvmdn/mdmJkj6UdKaZ7QpczjSzf5vZ2Hzbb2Nmq81sh5ktNLO/HZRtgJmtNLOdZjbhQDYAAKKRmT1gZpMOum+wmb1qZqeY2Sgz22RmGwPjbcxh1lPPzNID42e6mdXL91i5wHj/Q2Dsnxq4v7GZfR+4PkZSFUkzAmP4v8zsAzO796DtrDSz9oX8NABhi2IKBNeNkuIlnS3pPEmPmFkTSc9K6iLpDEkbJI0PLN9M0tWBZU8JLLM1/wqdc7slJUj6wTl3UuDyQ/5lzOw8Se9K6icpTtIs5Q2AJfIt1kVSc0nVJV0i6dZC+RsDABCexktqYWZlJClQPLtIGifpLUn7JZ0jqZbyxuOeB6/AzMpJ+kDSq5LKS3pJ0gdmVj6wyBhJpSXVlFRR0ssHr8M5d5OkbyW1DozhL0h6W1L3fNu5VFKlwLaAIoFiCgTXEOfcd865bZKeltRNeWX1TefcMudctqSBkuqaWTVJ+ySVkXSB8s4BX+Oc23QM271e0gfOubnOuX2SXpR0gqR6+ZZ51Tn3QyDbDEmXHdtfEQCA8Oec2yBpmaQDeyGbSNoj6RtJLST1c87tds79rLxC2fUQq2kp6Uvn3Bjn3H7n3LuS1kpqbWZnKO+D4zudc9udc/ucc4sKGG+6pPPM7NzA7ZskTXDO/XYMf1UgIlFMgeD6Lt/1DZLODFw2HLjTObdLeXtFKznnFkgaImmopJ/NbLiZnXwM2z14G7mBLJXyLfNjvut7JJ10DNsBACCSjFPeh8SSdEPgdlVJsZI2BU5/2SHpDeXt8TzYH8bXgA3KG1/PkrTNObf9aEM55/ZKmiCpu5kVC2Qcc7TrASIZxRQIrrPyXa8i6YfApeqBOwPnjJaXtFGSnHOvOucul3Sh8g7p/ech1vtX02kfvA0LZNl49H8FAACixnuSGptZZeXtOR2nvA9usyVVcM6dGric7JyreYif/8P4GlBFeePrd5LKmdmpBchxqHH8beUdVXWtpD3OubSC/IWAaEExBYLrbjOrHDgn5WHlfRr6rqQeZnaZmZWU9Iykpc659WZ2hZn93cxiJe2WtFdS7iHW+5Ok8mZ2ymG2O1FSSzO7NrCu/sobdFML968HAEDkcM5tlrRQ0v8kfZPvlJmPJP3HzE42s2JmdraZNTrEKmYp75DbG8ysuJldr7wPkmcG1vOhpGFmVtbMYs3s6sNE+UlSjYOypSlvzP+P2FuKIohiCgTXOOUNdl9LWifp/wLfWfaopMmSNilvYqQD57GcLGmEpO3KOzRoq6RBB6/UObdWeQX368BhR2ce9PjnyptE4b+StkhqrbxJFjhXBQBQ1I2T1DTw5wE3Syoh6TPljcGTlDdB4R8457ZKaqW8D3y3SvqXpFbOuS2BRW5S3nwRayX9rLxJCA/lWeVNiLjDzAbku3+0pIsljT30jwHRy5z7qyMCARwLvjwbAAAcDTO7WVJv51wD31mAUGOPKQAAAOCZmZWW1EfScN9ZAB8opgAAAIBHZhYvabPyzj0d9xeLA1GJQ3kBAAAAAF6xxxQAAAAA4BXFFAAAAADgVXHfAfKrUKGCq1atmu8YAIAokJmZucU5F+c7R6RjbAYAFJYjjc1hVUyrVaumjIwM3zEAAFHAzDb4zhANGJsBAIXlSGMzh/ICAAAAALyimAIAAAAAvKKYAgAAAAC8opgCAAAAALyimAIAAAAAvKKYAgAAAAC8opgCAAAAALyimAIAAAAAvKKYAgAAAAC8opgCAAAAALyimAIAAAAAvCruOwAAAAdkZmaqdOnS+tvf/uY7CgBEJeeccnNzf7/k5OT84fZf3eec+/1yYH0HXz/SY8ey3PGuu6DPS2EuF23bPu2001S7du0Cb/9YUEwBAGFh7ty56tChg2rVqqVFixbJzHxHAoBDcs4pOztbu3bt0q5du7R3794jXrKzsw97/759+7R//37t27fvD9eP9s/8ZfJQxfLA/cCxaN++vaZMmRLUbVBMAQDejR8/XjfffLMuuOACjR8/nlIKIGhyc3P1yy+/aNu2bdq+ffsfLgfu27Fjh3799dffi+euXbv+cPvXX389rpJXokQJlSxZUiVLllRsbKxiY2NVvHjxQ/4ZGxurkiVL6qSTTjrsMsWLF1fx4sVVrFixP1xiYmIK9b4DF0kys99/Vx/q+pEeO5bljnfdBVHYy0XTtsuWLVvg7R4riikAwKtXXnlF9913n66++mpNmzZNp556qu9IACKMc07btm3TDz/8oB9//FE//fSTfvzxxz9cP/Dn1q1blZube9h1lShRQmXLllWZMmV00kknqUyZMqpQoYKqVav2++2TTjrpD5cTTjhBpUqVUqlSpVSyZMnfrx/qUrJkyd/LHYD/j2IKAPDCOaeBAwfq+eefV4cOHfTOO++oVKlSvmMBCFNbt27VF198ofXr12vDhg2//7lhwwZ9++232r17959+plSpUjr99NN1+umnq0aNGqpXr57i4uJUrlw5lS1b9vdL/tsnnHACR20AHlBMAQAht2/fPvXq1Utvv/227rzzTg0ZMkQxMTG+YwHwLCcnR+vXr9fatWv/dNmyZcsfli1XrpyqVaumCy64QPHx8apataoqVaqk008/XaeddppOP/10lSlThpIJRAiKKQAgpHbv3q0uXbpo1qxZeuKJJ/Too4/yxhEogrKzs7V69Wp98sknWrZsmT755BOtWLFCe/bs+X2ZuLg4XXDBBWrfvr0uuOACnXfeeapevbqqVq2qk046yWN6AIWNYgoACJktW7aoVatWSk9P1+uvv6477rjDdyQAIbJlyxalpKQoJSVFycnJysjI0L59+yRJZcqU0WWXXaZevXrp4osv1t/+9jedf/75Kl++vOfUAEKFYgoACIkNGzYoPj5e69ev16RJk9S+fXvfkQAE0c6dOzV//nzNmTNHSUlJWrt2raS8yYXq1Kmjfv36qU6dOqpVq5bOPvtsJgQCijiKKQAg6FatWqXmzZtr9+7d+uijj3T11Vf7jgSgkDnntGzZMs2ePVuzZ89WWlqacnJyVKZMGV199dW65ZZb1KBBA9WpU4eJzgD8CcUUABBUSUlJatOmjU488UQtXrxYF198se9IAAqJc07p6emaMGGC3nvvPX333XeSpNq1a+uBBx5QfHy86tatq9jYWM9JAYQ7iikAIGjef/99devWTdWqVdOcOXNUtWpV35GilpmVkpQkqaTyxvdJzrnHD1qmpKTRki6XtFXS9c659SGOigjnnNMnn3yi8ePHa+LEidqwYYNiY2MVHx+vp556SgkJCapYsaLvmAAiDMUUABAUb7zxhvr06aMrrrhCM2fOVIUKFXxHinbZkpo453aZWaykZDP70Dm3JN8yt0va7pw7x8y6Snpe0vU+wiLybN++XePGjdPIkSO1fPlyFS9eXM2aNdMTTzyhtm3b6tRTT/UdEUAEo5gCAAqVc05PPfWUHn/8cbVo0UITJ07UiSee6DtW1HPOOUm7AjdjAxd30GJtJf07cH2SpCFmZoGfBQ5pxYoVeuWVVzR+/Hjt3btXtWvX1rBhw3T99derXLlyvuMBiBIUUwBAocnJydE999yj119/XbfccotGjBjBuWUhZGYxkjIlnSNpqHNu6UGLVJL0nSQ55/ab2U5J5SVtOWg9vSX1lqQqVaoEOzbCUG5urmbOnKlXXnlFiYmJOvHEE9WjRw/16tVLtWrV8h0PQBRiXm4AQKHYu3evunTpotdff10PPvig/ve//1FKQ8w5l+Ocu0xSZUlXmtlFx7ie4c65Os65OnFxcYWaEeEtNzdXEydO1KWXXqq2bdvqq6++0gsvvKDvvvtOw4YNo5QCCBr2mAIAjtuOHTvUtm1bJSUl6eWXX1a/fv18RyrSnHM7zCxRUnNJn+Z7aKOksyR9b2bFJZ2ivEmQUMTl5uZq0qRJeuKJJ/TZZ5/pggsu0NixY9WlSxc+YAIQEuwxBQAclx9++EFXX3210tLSNG7cOEqpJ2YWZ2anBq6fIOk6SWsPWmy6pFsC1ztJWsD5pUhOTtZVV12l66+/Xs45vfvuu/r000914403UkoBhAx7TAEAx+zzzz9XfHy8tm7dqg8++EDXXXed70hF2RmS3g6cZ1pM0kTn3Ewze1JShnNuuqRRksaY2VeStknq6i8ufPv666/1r3/9S5MnT1alSpX01ltvqXv37oqJifEdDUARRDEFAByTpUuXqmXLloqJidHChQt1+eWX+45UpDnnVkr60wmAzrnH8l3fK6lzKHMh/Pz222968cUX9dRTTykmJkZPPvmk+vfvr9KlS/uOBqAIo5gCAI7ahx9+qE6dOun000/XnDlzdM455/iOBKAAUlNT1bt3b61evVqdOnXS4MGDdeaZZ/qOBQCcYwoAODqjR49WmzZtdP755yslJYVSCkSAvXv36r777lODBg30yy+/aPr06XrvvfcopQDCBsUUAFAgzjkNGjRIt9xyixo1aqSFCxfq9NNP9x0LwF9YtWqVrrjiCr3yyivq06ePPvvsM7Vu3dp3LAD4A4opAOAv5ebmqn///vrXv/6lLl266IMPPtDJJ5/sOxaAI3DOafDgwbriiiu0efNmzZo1S0OGDNFJJ53kOxoA/AnnmAIAjui3335Tjx49NG7cOPXt21cvv/yyihXjc00gnP3666/q0aOHJk+erNatW2vkyJGqWLGi71gAcFgUUwDAYf3666/q2LGj5s6dq2effVYPPPCAzMx3LABH8MUXX6h9+/Zau3atXnzxRd1///38vwUQ9iimAIBD+vnnn9WiRQstX75cb775pnr06OE7EoC/MHPmTN14442KjY3VRx99pGuvvdZ3JAAoEI7FAgD8yddff6369evrs88+09SpUymlQAQYNmyY2rRpo7PPPluZmZmUUgARJah7TM1svaRfJeVI2u+cqxPM7QEAjt8nn3yihIQE7du3T/Pnz1fdunV9RwJwBM45PfLII3rmmWfUqlUrTZgwQaVLl/YdCwCOSigO5b3GObclBNsBABynBQsWqF27djr11FOVmJiov/3tb74jATiCffv2qXfv3nrrrbfUq1cvDRs2TMWLc6YWgMjDobwAAEnSxIkTlZCQoCpVqig1NZVSCoS57OxsderUSW+99ZaeeOIJvfHGG5RSABEr2MXUSfrIzDLNrHeQtwUAOEZDhgxR165ddeWVV2rx4sWqXLmy70gAjiA7O1udO3fW9OnTNXToUD322GPMvAsgogW7mDZwztWWlCDpbjO7+uAFzKy3mWWYWcbmzZuDHAcAkJ9zTg8//LDuvfdetWnTRh999JHKli3rOxaAIzhQSmfMmKFhw4apT58+viMBwHELajF1zm0M/PmzpPclXXmIZYY75+o45+rExcUFMw4AIJ/9+/erZ8+eeuaZZ9SrVy9NmjRJJ5xwgu9YAI7g4FJ61113+Y4EAIUiaMXUzE40szIHrktqJunTYG0PAFBwe/bsUYcOHfTmm2/q0Ucf5dw0IALk5ubqlltuoZQCiErBfBdymqT3A+c7FJc0zjk3O4jbAwAUwLZt29SqVSstWbKEN7dABBkwYIAmTJig559/nv+3AKJO0Iqpc+5rSZcGa/0AgKP33XffKT4+XuvWrdPEiRPVqVMn35EAFMB//vMfvfzyy+rbt6/++c9/+o4DAIWO47YAoIhYvXq1mjdvrl9++UVz5sxR48aNfUcCUADjxo3TgAED1LlzZ7388svMvgsgKvE9pgBQBKSkpKhhw4bav3+/kpKSKKVAhEhLS9Ott96qRo0aafTo0SpWjLduAKITv90AIMpNnz5dTZs2VYUKFZSamqpLL+UsCyAS/PDDD+rYsaOqVKmiKVOmqFSpUr4jAUDQUEwBIIqNGjVK7du318UXX6yUlBRVr17ddyQABZCdna2OHTvql19+0dSpU1WuXDnfkQAgqCimABCFnHN6+umn1bNnT1133XVasGCB+K5oIHLce++9WrJkid566y1ddNFFvuMAQNAx+REARJmcnBz94x//0NChQ9W9e3e9+eabio2N9R0LQAENHz5cI0aM0MCBA5k5G0CRwR5TAIgi2dnZ6tq1q4YOHaoBAwbo7bffppQCEWTVqlXq27ev4uPj9dRTT/mOAwAhwx5TAIgSO3fuVPv27ZWYmKgXX3xR/fv39x0JwFHIyspSt27ddOqpp2r06NGKiYnxHQkAQoZiCgBRYNOmTUpISNDq1as1duxY3Xjjjb4jAThK//znP7V69WrNnj1bFStW9B0HAEKKYgoAEe7LL79Us2bNtHnzZs2cOVPx8fG+IwE4SjNmzNDQoUN1//33838YQJFEMQWACJaRkaEWLVrIOafExERdccUVviMBOEo//PCDevToocsuu0zPPPOM7zgA4AWTHwFAhProo4/UuHFjnXjiiUpJSaGUAhHIOac77rhDe/bs0bvvvquSJUv6jgQAXlBMASACvfPOO2rZsqXOOeccpaam6rzzzvMdCcAxmDBhgmbOnKmnnnpKF1xwge84AOANxRQAIsxLL72k7t27q0GDBlq0aJHOOOMM35EAHIOtW7eqb9++qlOnjv7xj3/4jgMAXnGOKQBEiNzcXD344IMaNGiQOnXqpDFjxqhUqVK+YwE4Rvfdd5+2b9+uefPmqXhx3pIBKNrYYwoAEWDfvn269dZbNWjQIN19990aP348pRSIYLNnz9aYMWP04IMP6pJLLvEdBwC84+M5AAhzu3btUufOnTV79mw99dRTevjhh2VmvmMBOEa7du3SHXfcoQsuuECPPPKI7zgAEBYopgAQxrZs2aKWLVsqIyNDI0aMUM+ePX1HAnCcnnnmGX377bdKTk5mFl4ACKCYAkCYWr9+veLj4/Xtt9/q/fffV5s2bXxHAnCc1q1bp//85z/q3r276tev7zsOAIQNiikAhKGVK1eqefPmysrK0rx583gDC0SJ/v37KzY2Vs8//7zvKAAQVpj8CADCzKJFi9SwYUPFxMQoOTmZUgpEiY8++kjTpk3TI488ojPPPNN3HAAIKxRTAAgjU6ZMUXx8vCpVqqTU1FTVrFnTdyREADM7y8wSzewzM1ttZn/6Ukwza2xmO81seeDymI+sRdW+ffvUr18/nX322brvvvt8xwGAsMOhvAAQJl577TXdfffduuqqqzRz5kyVK1fOdyREjv2S+jvnlplZGUmZZjbXOffZQcstds618pCvyBs6dKjWrFmjadOmMeERABwCe0wBwDPnnB5//HH16dNHLVu21Lx58yilOCrOuU3OuWWB679KWiOpkt9UOGD79u164okn1KxZM7Vu3dp3HAAISxRTAPBo//79uvPOO/Xkk0/qtttu0/vvv6/SpUv7joUIZmbVJNWStPQQD9c1sxVm9qGZHfY4cTPrbWYZZpaxefPmYEUtMgYNGqQdO3bohRde4DuIAeAwKKYA4ElWVpY6d+6s4cOH66GHHtLIkSNVvDhnWODYmdlJkiZL6uec++Wgh5dJquqcu1TSfyVNPdx6nHPDnXN1nHN14uLigpa3KPjxxx81ePBgdevWTZdeeqnvOAAQtiimAODB9u3b1axZM02bNk2vvvqqnn76afak4LiYWazySuk7zrkpBz/unPvFObcrcH2WpFgzqxDimEXO//3f/+m3337Tk08+6TsKAIQ1PpoHgBDbuHGjmjdvri+++ELjx49Xly5dfEdChLO8TzVGSVrjnHvpMMucLukn55wzsyuV9+H01hDGLHK++eYbDR8+XLfffrvOOecc33EAIKxRTAEghNasWaP4+Hjt2LFDH374oZo0aeI7EqJDfUk3SVplZssD9z0kqYokOedel9RJ0l1mtl9SlqSuzjnnIWuR8fjjjysmJkaPPvqo7ygAEPYopgAQImlpaWrVqpViY2O1aNEi1apVy3ckRAnnXLKkIx4L7pwbImlIaBLh008/1dixYzVgwABVqsQEyQDwVzjHFABC4IMPPtC1116rcuXKKTU1lVIKRLknnnhCZcqU0QMPPOA7CgBEBIopAATZW2+9pbZt2+rCCy9USkqKatSo4TsSgCBas2aNJk+erHvvvVfly5f3HQcAIgLFFACCxDmn5557Tj169FCTJk2UmJioihUr+o4FIMief/55nXDCCfrHP/7hOwoARAyKKQAEQW5urvr166eBAweqW7dumjlzpsqUKeM7FoAgW79+vcaOHavevXuL74AFgIKjmAJAIcvOztaNN96oV199Vf369dPYsWNVokQJ37EAhMCgQYNUrFgx9e/f33cUAIgozMoLAIXol19+UYcOHTR//ny98MILGjBggPK+YhJAtPvxxx81atQo3XLLLapcubLvOAAQUSimAFBIfvrpJyUkJGjlypV6++23dfPNN/uOBCCEXnrpJe3bt4+ZeAHgGFBMAaAQrFu3TvHx8dq0aZNmzJihhIQE35EAhNC2bdv02muv6frrr9c555zjOw4ARByKKQAcp2XLlikhIUE5OTlasGCB/v73v/uOBCDEhg8frl27dunBBx/0HQUAIhKTHwHAcZg3b54aNWqkUqVKKSUlhVIKFEH79u3T0KFDde211+qSSy7xHQcAIhLFFACO0fjx49WiRQtVr15dqampOv/8831HAuDBlClT9P3336tfv36+owBAxKKYAsAxGDx4sLp166a6desqKSlJlSpV8h0JgCeDBw/WOeecoxYtWviOAgARi2IKAEfBOaeBAweqX79+6tChg+bMmaNTTz3VdywAnnz88cdKS0vTvffeq2LFeFsFAMeKyY8AoID27dunXr166e2339add96pIUOGKCYmxncsAB4NHjxYJ598snr06OE7CgBEND7aA4AC2L17t9q1a6e3335bTzzxhIYNG0YpBYq4H374QRMnTtRtt92mMmXK+I4DABEt6HtMzSxGUoakjc65VsHeHgAUti1btqhVq1ZKT0/XG2+8od69e/uOBCAMDBs2TDk5Obr33nt9RwGAiBeKQ3n/IWmNpJNDsC0AKFQbNmxQfHy81q9fr8mTJ6tdu3a+IwEIA9nZ2XrjjTfUpk0b1ahRw3ccAIh4QT2U18wqS2opaWQwtwMAwbBq1SrVq1dPP/30k+bOnUspBfC7999/X1u2bFGfPn18RwGAqBDsc0xfkfQvSbmHW8DMeptZhpllbN68OchxAKBgkpKS1LBhQ5mZFi9erIYNG/qOBCCMDB8+XNWrV1fTpk19RwGAqBC0YmpmrST97JzLPNJyzrnhzrk6zrk6cXFxwYoDAAX2/vvvq1mzZjrjjDOUmpqqiy66yHckAGHkyy+/VGJionr27MlXxABAIQnmb9P6ktqY2XpJ4yU1MbOxQdweABy3119/XZ06dVKtWrWUnJysKlWq+I4EIMyMHDlSMTExfEUMABSioBVT59xA51xl51w1SV0lLXDOdQ/W9gDgeDjn9O9//1t33XWXEhISNH/+fJUvX953LABh5rffftP//vc/tWnTRmeccYbvOAAQNUIxKy8AhLWcnBz16dNHw4cPV48ePfTGG28oNjbWdywAYWjatGnavHmzevXq5TsKAESVkBRT59xCSQtDsS0AOBpZWVm64YYbNHXqVA0cOFBPP/20zMx3LABhasSIEapSpYqaNWvmOwoARBXO2AdQZO3YsUPx8fGaNm2aXn31VT3zzDOUUgCH9fXXX2vu3Lnq2bOnYmJifMcBgKjCobwAiqSNGzeqefPm+vzzzzV+/Hh16dLFdyQAYW7kyJEqVqyYbrvtNt9RACDqUEwBFDlr165VfHy8tm/frtmzZ6tJkya+IwEIczk5ORo9erQSEhJUqVIl33EAIOpwKC+AImXJkiWqX7++srOztWjRIkopgAJZuHChNm7cqJtvvtl3FACIShRTAEXGBx98oCZNmqhcuXJKTU1VrVq1fEcCECHGjBmjk08+Wa1bt/YdBQCiEsUUQJHwv//9T23bttWFF16olJQU1ahRw3ckABFi9+7dmjRpkjp37qwTTjjBdxwAiEoUUwBRzTmn5557TrfddpuaNGmixMREVaxY0XcsABFk6tSp2r17t2666SbfUQAgalFMAUSt3Nxc9evXTwMHDtQNN9ygmTNnqkyZMr5jAYgwY8aMUdWqVdWwYUPfUQAgalFMAUSl7Oxs3XDDDXr11Vd13333acyYMSpRooTvWAAizKZNmzR37lx1795dxYrxtgkAgoWviwEQdX755Rd16NBB8+fP16BBgzRgwADfkQBEqHHjxik3N5fDeAEgyCimAKLKjz/+qBYtWmjVqlUaPXo0byYBHJcxY8boiiuu0Pnnn+87CgBENYopgKjx1VdfKT4+Xj/++KNmzJih5s2b+44EIIKtWrVKK1as0Kuvvuo7CgBEPYopgKiQmZmphIQEOeeUmJioK6+80nckABHu3XffVUxMjLp27eo7CgBEPc7iBxDx5s6dq8aNG6t06dJKSUmhlAI4bs45TZw4Uddee63i4uJ8xwGAqEcxBRDR3n33XbVs2VI1atRQamqqzjvvPN+RgJAzs7PMLNHMPjOz1Wb2j0MsY2b2qpl9ZWYrzay2j6yR4pNPPtG6devUuXNn31EAoEigmAKIWC+//LJuuOEG1atXT0lJSTrzzDN9RwJ82S+pv3PuQklXSbrbzC48aJkESecGLr0lvRbaiJFl4sSJiomJUfv27X1HAYAigWIKIOI45/TAAw/o/vvvV8eOHTV79mydcsopvmMB3jjnNjnnlgWu/yppjaRKBy3WVtJol2eJpFPN7IwQR40IBw7jbdq0qcqXL+87DgAUCRRTABFl3759uvXWW/XCCy+oT58+mjBhgkqVKuU7FhA2zKyapFqSlh70UCVJ3+W7/b3+XF4hadmyZfrmm2/UpUsX31EAoMigmAKIGLt371bbtm01evRoPfXUUxoyZIhiYmJ8xwLChpmdJGmypH7OuV+OYz29zSzDzDI2b95ceAEjxMSJE1W8eHG1a9fOdxQAKDL4uhgAEWHLli1q2bKlMjIyNGLECPXs2dN3JCCsmFms8krpO865KYdYZKOks/Ldrhy470+cc8MlDZekOnXquEKOGtYOHMZ73XXXqVy5cr7jAECRwR5TAGFv/fr1ql+/vlauXKkpU6ZQSoGDmJlJGiVpjXPupcMsNl3SzYHZea+StNM5tylkISNERkaG1q9fz2G8ABBi7DEFENZWrlyp5s2bKysrS/PmzVP9+vV9RwLCUX1JN0laZWbLA/c9JKmKJDnnXpc0S1ILSV9J2iOpR+hjhr/33ntPsbGxatu2re8oAFCkUEwBhK1FixapTZs2Ovnkk5WcnKyaNWv6jgSEJedcsiT7i2WcpLtDkygy5T+Mt2zZsr7jAECRwqG8AMLS5MmTFR8fr0qVKik1NZVSCiDoli1bpg0bNqhz586+owBAkUMxBRB2XnvtNXXu3Fm1a9dWcnKyzjrrrL/+IQA4TlOnTlWxYsXUunVr31EAoMihmAIIG845PfbYY+rTp49atmypefPmMSsmgJCZNm2aGjZsqPLly/uOAgBFDsUUQFjYv3+/7rjjDj311FO67bbb9P7776t06dK+YwEoIr7++mutWrWKSY8AwBOKKQDvsrKy1KlTJ40YMUIPP/ywRo4cqeLFmZsNQOhMmzZNkiimAOAJ7/wAeLV9+3a1adNGKSkp+u9//6t77rnHdyQARdDUqVN18cUXq0aNGr6jAECRxB5TAN58//33atiwoT7++GNNmDCBUgrAiy1btig5OZm9pQDgEXtMAXixZs0axcfHa8eOHZo9e7auueYa35EAFFEffPCBcnNz1a5dO99RAKDIYo8pgJBLTU1V/fr1tW/fPiUlJVFKAXg1depUVa5cWbVr1/YdBQCKLIopgJCaMWOGmjZtqgoVKig1NVWXXXaZ70gAirA9e/Zozpw5atu2rczMdxwAKLIopgBC5s0331T79u1Vs2ZNpaSkqHr16r4jASji5s2bp6ysLM4vBQDPKKYAgs45p2eeeUa33367mjZtqsTERMXFxfmOBQCaNm2aTj75ZDVq1Mh3FAAo0iimAIIqJydHffv21cMPP6zu3btr+vTpOumkk3zHAgDl5uZq5syZatGihUqUKOE7DgAUaRRTAEGTnZ2tbt26aciQIerfv7/efvtt3vwBCBvLli3Tzz//rJYtW/qOAgBFHl8XAyAodu7cqfbt2ysxMVEvvvii+vfv7zsSAPzBhx9+KDNTfHy87ygAUORRTAEUuh9//FEJCQn69NNPNWbMGHXv3t13JAD4k1mzZumKK67gnHcACAMcygugUH355ZeqV6+evvzyS82cOZNSCiAsbdmyRUuXLlWLFi18RwEAiD2mAApRenr672/yEhMTdcUVV3hOBACHNmfOHDnnKKYAECbYYwqgUMyZM0fXXHONTjrpJKWkpFBKAYS1WbNmKS4uTpdffrnvKAAAUUwBFIJ33nlHrVq10jnnnKPU1FSdd955viMBwGHl5ORozpw5at68uYoV460QAIQDfhsDOC7/+c9/1L17dzVo0ECLFi3SGWec4TsSABxRenq6tm7dymG8ABBGglZMzayUmX1sZivMbLWZPRGsbQEIvdzcXA0YMEADBgxQ586dNXv2bJ1yyim+YwHAX5o1a5aKFSumZs2a+Y4CAAgI5uRH2ZKaOOd2mVmspGQz+9A5tySI2wQQAvv27dNtt92msWPH6u6779bgwYMVExPjOxYAFMisWbNUt25dlStXzncUAEBA0PaYujy7AjdjAxcXrO0BCI1du3apdevWGjt2rJ5++mn997//pZQCiBg//vijMjMzOYwXAMJMUL8uxsxiJGVKOkfSUOfc0kMs01tSb0mqUqVKMOMAOE6bN29Wy5YtlZmZqZEjR+r222/3HQkAjsqcOXMkSQkJCZ6TAADyC+rkR865HOfcZZIqS7rSzC46xDLDnXN1nHN14uLighkHwHH45ptvVL9+fa1atUpTp06llAKISLNnz9Zpp52myy67zHcUAEA+IZmV1zm3Q1KipOah2B6AwrV8+XLVq1dPW7Zs0fz589W6dWvfkQDgqOXm5mrevHm67rrrZGa+4wAA8gnmrLxxZnZq4PoJkq6TtDZY2wMQHImJiWrUqJGKFy+u5ORk1atXz3ckADgmK1as0JYtW3Tdddf5jgIAOEgw95ieISnRzFZKSpc01zk3M4jbA1DI3nvvPTVv3lyVK1dWWlqaLrzwQt+RAOCYzZ07V5LUtGlTz0kAAAcL2uRHzrmVkmoFa/0Agmvo0KG69957Va9ePU2fPp2vVQAQ8ebOnauaNWvqzDPP9B0FAHCQkJxjCiByOOf0yCOP6J577lHr1q01d+5cSimAiJeVlaXFixdzGC8AhKmgfl0MgMiyf/9+3XnnnRo1apR69uyp1157TcWL82sCQORLTk5WdnY2xRQAwhR7TAFIkvbs2aMOHTpo1KhRevTRRzV8+HBKKYCoMXfuXMXGxqpRo0a+owAADoF3nQC0bds2tW7dWmlpaRo6dKj69OnjOxIAFKq5c+eqXr16OvHEE31HAQAcwl/uMTWzwcaXfQFR67vvvlODBg2UkZGhiRMnUkoBz8wsxcyu8Z0jmvz8889avnw5h/ECQBgryKG8v0qabmYnSpKZxZtZSnBjAQiF1atXq169etq4caPmzJmjTp06+Y4EQLpD0j1mNt/M6voOEw3mz58vSRRTAAhjf3kor3PuETO7QdJCM/tN0i5JDwY9GYCgSklJUatWrVSqVCklJSXp0ksv9R0JgCTn3KeSOppZbUlPBg5aetg5t8Jvssg1d+5clS1bVpdffrnvKACAwyjIobzXSuolabekCpL6OucWBzsYgOCZPn26mjZtqooVKyotLY1SCoSnryQ9JWmjpEzPWSKWc05z585VkyZNFBMT4zsOAOAwCnIo78OSHnXONZbUSdIEM2sS1FQAgmbkyJFq3769LrnkEiUnJ6tatWq+IwHIx8wSzex75ZXRhyRtk3Sr11AR7PPPP9f333/PYbwAEOYKcihvk3zXV5lZgqTJkuoFMxiAwuWc09NPP61HH31UzZs313vvvaeTTjrJdywAf9Zf0hrnXNbR/JCZvSmplaSfnXMXHeLxxpKmSfomcNcU59yTxxc1/C1YsECSdO2113pOAgA4kqP+uhjn3KbA4b0AIkROTo769u2rYcOG6aabbtKoUaMUGxvrOxaAQ3DOLTvGH31L0hBJo4+wzGLnXKtjXH9ESkxMVOXKlXX22Wf7jgIAOIKCHMr7J0f7KS4Af/bu3auuXbtq2LBh+uc//6m33nqLUgpEIedckvIO+0WAc04LFy7UNddcI775DgDC2zEVUwCRYefOnWrevLkmTZqkl156SS+88IKKFeO/PVCE1TWzFWb2oZnVPNxCZtbbzDLMLGPz5s2hzFeoVq9erS1btuiaa/haWAAId0d9KC+AyPDDDz8oISFBa9as0TvvvKMbbrjBdyQAfi2TVNU5t8vMWkiaKuncQy3onBsuabgk1alTx4UsYSFLTEyUJIopAEQAdp0AUeiLL75Q/fr1tW7dOs2cOZNSCkDOuV+cc7sC12dJijWzCp5jBVViYqKqVavG7OMAEAEopkCUWbp0qerVq6fdu3dr4cKFatasme9IAMKAmZ1ugRMtzexK5b0H2Oo3VfDk5uZq0aJF7C0FgAjBobxAFPnggw/UpUsXnXHGGZo9e7bOOecc35EAhIiZvSupsaQKge9BfVxSrCQ5515X3neR32Vm+yVlSerqnIvYw3T/ysqVK7Vt2zaKKQBECIopECX+97//qVevXrr00ks1a9YsnXbaab4jAQgh51y3v3h8iPK+TqZIOHB+aePGjf0GAQAUCIfyAhHOOadnnnlGt912m5o0aaKFCxdSSgEUeYmJiTr77LN11lln+Y4CACgAiikQwXJyctS3b189/PDDuvHGGzVz5kyVKVPGdywA8ConJ0dJSUkcxgsAEYRiCkSovXv3qmvXrhoyZIgGDBig0aNHq0SJEr5jAYB3y5cv186dOymmABBBOMcUiEA7duxQu3bttGjRIv3nP//R/fff7zsSAIQNvr8UACIPxRSIMBs3blRCQoLWrl2rcePGqVu3I853AgBFTmJios4//3ydccYZvqMAAAqIQ3mBCLJmzRrVq1dP33zzjWbNmkUpBYCD7N+/X4sXL2ZvKQBEGPaYAhEiLS1NrVq1UmxsrJKSklSrVi3fkQAg7KxYsUK//vqrrr76at9RAABHgT2mQASYMWOGrr32WpUrV06pqamUUgA4jMWLF0uSGjZs6DkJAOBoUEyBMDdy5Ei1a9dOF110kVJSUlSjRg3fkQAgbC1evFjVq1dX5cqVfUcBABwFiikQppxz+r//+z/16tVLzZo104IFC1SxYkXfsQAgbDnntHjxYvaWAkAE4hxTIAzl5OTonnvu0euvv66bb75ZI0eOVGxsrO9YABDWPv/8c23evJliCgARiD2mQJjJyspS586d9frrr+vBBx/UW2+9RSkFgAJISkqSJCY+AoAIxB5TIIxs375dbdq0UUpKigYPHqy+ffv6jgQAEWPx4sWqWLGizj33XN9RAABHiWIKhInvvvtOCQkJ+vLLLzV+/Hh16dLFdyQAiCgHzi81M99RAABHiUN5gTCwevVq1atXT99++61mz55NKQWAo/Ttt99qw4YNHMYLABGKYgp4lpycrAYNGmj//v1KSkrSNddc4zsSAEQcvr8UACIbxRTwaOrUqbruuutUsWJFpaWl6bLLLvMdCQAi0uLFi3XyySfrkksu8R0FAHAMKKaAJ2+88YY6duyoSy+9VCkpKapWrZrvSAAQsZKSklS/fn3FxMT4jgIAOAYUUyDEnHP697//rTvvvFPNmzfX/PnzVaFCBd+xACBibdmyRWvWrOEwXgCIYMzKC4TQ/v371adPH40YMUI9evTQG2+8wXeUAsBxSk5OlsT3lwJAJGOPKRAie/bsUceOHTVixAg9/PDDGjVqFKUUAApBUlKSSpYsqTp16viOAgA4RuwxBUJg27Ztat26tdLS0jR06FD16dPHdyQAiBopKSm68sorVbJkSd9RAADHiD2mQJB9++23atCggTIyMvTee+9RSgGgEGVlZWnZsmWqX7++7ygAgOPAHlMgiFatWqXmzZtr9+7d+uijj9SoUSPfkQAgqqSnp2v//v2qV6+e7ygAgOPAHlMgSJKSkn6fIXLx4sWUUgAIgtTUVElS3bp1PScBABwPiikQBJMnT1azZs10xhlnKC0tTRdffLHvSAAQlVJSUnT++efztVsAEOGCVkzN7CwzSzSzz8xstZn9I1jbAsLJsGHD1LlzZ9WuXVvJycmqUqWK70gAEJWcc0pNTeX8UgCIAsHcY7pfUn/n3IWSrpJ0t5ldGMTtAV455/TII4/o7rvvVuvWrTVv3jyVL1/edywAiFqff/65tm3bxvmlABAFglZMnXObnHPLAtd/lbRGUqVgbQ/waf/+/erZs6eefvpp9erVS5MnT1bp0qV9xwKAqHbg/FL2mAJA5AvJOaZmVk1SLUlLQ7E9IJR2796tdu3a6c0339Rjjz2mN954Q8WLM+E1AARbSkqKypUrp/POO893FADAcQr6u2czO0nSZEn9nHO/HOLx3pJ6S+JcPEScLVu2qFWrVkpPT9drr72mO++803ckACgyUlNTVa9ePRUrxlyOABDpgvqb3MxilVdK33HOTTnUMs654c65Os65OnFxccGMAxSq9evXq0GDBlq+fLkmTZpEKQWAENq6davWrl3L+aUAECWCtsfUzEzSKElrnHMvBWs7gA+ffPKJWrRooezsbM2dO/f37ysFAIRGWlqaJM4vBYBoEcw9pvUl3SSpiZktD1xaBHF7QEjMnTtXV199tWJjY5WcnEwpBQAPUlJSVLx4cdWpU8d3FABAIQjaHlPnXLIkC9b6AR/GjBmj2267TRdeeKFmzZqlSpWYaBoAfEhNTVXt2rWZAR0AogSzBQAF4JzTs88+q5tvvllXX321kpKSKKUA4Mlvv/2mjz/+mPNLASCKUEyBv5CTk6N77rlHDz30kG644QZ9+OGHOuWUU3zHAoAia/ny5dq7dy/nlwJAFKGYAkeQlZWlzp07a9iwYfrnP/+pMWPGqESJEr5jAUCRdmDio7p163pOAgAoLBRT4DC2bt2qpk2baurUqRo8eLBeeOEFvisPQNgyszfN7Gcz+/Qwj5uZvWpmX5nZSjOrHeqMhWXJkiWqXLkyp1QAQBThXTZwCOvXr1f9+vWVmZmpiRMnqm/fvr4jAcBfeUtS8yM8niDp3MClt6TXQpApKJYuXaqrrrrKdwwAQCGimAIHWb58uerWrauffvpJH330kTp16uQ7EgD8JedckqRtR1ikraTRLs8SSaea2RmhSVd4fv75Z33zzTcUUwCIMhRTIJ958+b94TtKr776at+RAKCwVJL0Xb7b3wfu+xMz621mGWaWsXnz5pCEK6ilS5dKkv7+9797TgIAKEwUUyBg7NixSkhIULVq1ZSWlqaaNWv6jgQAXjjnhjvn6jjn6sTFxfmO8wdLlixR8eLFVbt2xJ4iCwA4BIopijznnJ5//nnddNNNatiwoRYvXsyEGgCi0UZJZ+W7XTlwX0RZsmSJLrnkEpUuXdp3FABAIaKYokjLyclR37599eCDD6pr1658RymAaDZd0s2B2XmvkrTTObfJd6ijkZOTo/T0dM4vBYAoVNx3AMCXrKwsde/eXVOmTNGAAQP0/PPP83UwACKWmb0rqbGkCmb2vaTHJcVKknPudUmzJLWQ9JWkPZJ6+El67NasWaNff/2VYgoAUYhiiiJp27Ztatu2rVJSUvTyyy+rX79+viMBwHFxznX7i8edpLtDFCcomPgIAKIXxRRFzoYNG5SQkKB169Zp/Pjx6tKli+9IAIACWLJkicqWLatzzz3XdxQAQCGjmKJIWbFihRISErRnzx599NFHatSoke9IAIACWrp0qf7+97/LzHxHAQAUMk6oQ5Exf/58NWzYUDExMUpOTqaUAkAE+fXXX/Xpp59yfikARCmKKYqEcePGKSEhQVWrVlVaWpouuugi35EAAEchPT1dzjnOLwWAKEUxRVRzzmnQoEG68cYbVa9ePS1evFiVK1f2HQsAcJQOTHx05ZVXek4CAAgGiimiVk5Ojvr166d//etf6tKli+bMmaNTTz3VdywAwDFYsmSJzj//fJUrV853FABAEFBMEZWysrLUuXNnvfrqq7rvvvv07rvvqmTJkr5jAQCOgXPu94mPAADRiVl5EXU2b96sNm3aaOnSpXxHKQBEgQ0bNuinn35i4iMAiGIUU0SVr776SgkJCfr+++81adIkdejQwXckAMBxOnB+KXtMASB6UUwRNZYsWaLWrVvLOacFCxaobt26viMBAArBkiVLdMIJJ+jiiy/2HQUAECScY4qo8P777+uaa67RKaecorS0NEopAESRpUuXqk6dOoqNjfUdBQAQJBRTRLz//ve/6tixoy699FKlpaXp3HPP9R0JAFBIsrOztWzZMg7jBYAoRzFFxMrNzVX//v3Vt29ftW3bVgsWLFBcXJzvWACAQrRixQplZ2cz8REARDnOMUVEysrK0k033aTJkyfr3nvv1csvv6yYmBjfsQAAhYyJjwCgaKCYIuJs2bJFbdu2VVpaml566SX169dPZuY7FgAgCJYsWaJKlSqpcuXKvqMAAIKIYoqIsm7dOiUkJOjbb7/VxIkT1alTJ9+RAABBtHTpUvaWAkARwDmmiBgff/yx6tatq61bt2r+/PmUUgCIclu3btW6desopgBQBFBMERGmT5+uxo0bq0yZMkpLS1P9+vV9RwIABFlGRoYk6YorrvCcBAAQbBRThL2hQ4eqffv2uuiii5SWlqbzzjvPdyQAQAikp6fLzHT55Zf7jgIACDKKKcJWbm6u/vWvf+mee+5Rq1atlJiYqIoVK/qOBQAIkfT0dJ1//vk6+eSTfUcBAAQZxRRhKSsrS926ddOgQYN09913a8qUKTrxxBN9xwIAhFB6ejqH8QJAEcGsvAg7mzdvVtu2bbVkyRINGjRI/fv35+tgAKCI2bhxozZt2kQxBYAigmKKsLJ27Vq1bNlSmzZt0qRJk9ShQwffkQAAHqSnp0ti4iMAKCoopggbiYmJ6tChg0qUKKGFCxfqyiuv9B0JAOBJenq6ihcvrssuu8x3FABACHCOKcLCW2+9pWbNmunMM8/U0qVLKaUAUMSlp6fr4osvVqlSpXxHAQCEAMUUXjnn9Oijj6pHjx5q3LixUlJSVK1aNd+xAAAeOeeUkZHBYbwAUIRwKC+82bt3r2677Ta9++67uv322/Xaa68pNjbWdywAgGfr1q3T9u3bKaYAUIRQTOHFli1b1K5dO6WkpOjZZ5/VAw88wMy7AABJTHwEAEURxRQh98UXX6hFixb6/vvvNWHCBHXp0sV3JABAGElPT9cJJ5ygmjVr+o4CAAgRiilCKikpSe3atVNMTIwSExNVt25d35EAAGEmPT1dtWrVUvHivE0BgKKCyY8QMmPGjFHTpk112mmnaenSpZRSAMCf7N+/X5mZmRzGCwBFDMUUQeec0+OPP66bb75ZDRo0UGpqqmrUqOE7FgAgDH322WfKysqimAJAERO0Ympmb5rZz2b2abC2gfCXlZWlrl276sknn9Stt96q2bNnq2zZsr5jAQDCFBMfAUDRFMw9pm9Jah7E9SPMbdq0SY0aNdJ7772n559/Xm+++aZKlCjhOxYAIIylp6frlFNO0TnnnOM7CgAghII2q4BzLsnMqgVr/Qhvn3zyiVq3bq0dO3bo/fffV9u2bX1HAgBEgPT0dNWpU0fFinG2EQAUJfzWR6GbMmWKGjRooGLFiik5OZlSCgAokL1792rlypUcxgsARZD3Ympmvc0sw8wyNm/e7DsOjoNzTs8884w6duyoiy++WB9//LEuu+wy37EAABFixYoV2r9/P8UUAIog78XUOTfcOVfHOVcnLi7Odxwco7179+rmm2/Www8/rG7dumnhwoU6/fTTfccCAEQQJj4CgKLLezFF5Pvpp5/UpEkTjR07Vk899ZTeeecdlSpVyncsAChSzKy5mX1uZl+Z2YOHePxWM9tsZssDl54+ch5Jenq6TjvtNFWuXNl3FABAiAVt8iMze1dSY0kVzOx7SY8750YFa3vwY9WqVWrVqpU2b96s9957T506dfIdCQCKHDOLkTRU0nWSvpeUbmbTnXOfHbToBOfcPSEPWEDp6em64oorZGa+owAAQiyYs/J2C9a6ER5mzpypbt266eSTT1ZSUpLq1KnjOxIAFFVXSvrKOfe1JJnZeEltJR1cTMPWr7/+qrVr16pr166+owAAPOBQXhw155yeffZZtWnTRuedd54+/vhjSikA+FVJ0nf5bn8fuO9gHc1spZlNMrOzQhOtYDIzM+Wc4/xSACiiKKY4Krt371a3bt300EMP6frrr9fixYtVqdKh3vsAAMLMDEnVnHOXSJor6e3DLehjxvyMjAxJ4oNOACiiKKYosA0bNqhBgwaaOHGinn/+eY0bN06lS5f2HQsAIG2UlH8PaOXAfb9zzm11zmUHbo6UdPnhVuZjxvyMjAxVrVpVzNAPAEVT0M4xRXRZtGiROnXqpH379mnmzJlq0aKF70gAgP8vXdK5ZlZdeYW0q6Qb8i9gZmc45zYFbraRtCa0EY8sMzNTl19+2K4MAIhy7DHFETnnNGzYMDVt2lQVKlTQxx9/TCkFgDDjnNsv6R5Jc5RXOCc651ab2ZNm1iawWF8zW21mKyT1lXSrn7R/tmPHDn311VccxgsARRh7THFY2dnZuueeezRy5Ei1bNlS77zzjk455RTfsQAAh+CcmyVp1kH3PZbv+kBJA0OdqyCWLVsmSewxBYAijD2mOKQff/xRTZo00ciRI/XQQw9p2rRplFIAQFBkZmZKopgCQFHGHlP8ydKlS9WxY0dt375dEydOVOfOnX1HAgBEsczMTFWrVk3ly5f3HQUA4Al7TPE755xef/11NWzYUCVKlFBKSgqlFAAQdBkZGewtBYAijmIKSdKePXt066236q677tJ1112nzMxMXXbZZb5jAQCi3I4dO7Ru3TqKKQAUcRRTaN26dapbt67GjBmjJ554QjNmzFDZsmV9xwIAFAFMfAQAkDjHtMibOXOmunfvrmLFiumDDz5QQkKC70gAgCIkIyNDEsUUAIo69pgWUTk5OXr00UfVunVrnX322crMzKSUAgBCjomPAAASe0yLpK1bt+qGG27QRx99pNtuu01Dhw5VqVKlfMcCABRBmZmZ7C0FALDHtKhJT0/X5ZdfroULF2rEiBEaNWoUpRQA4MX27du1bt061alTx3cUAIBnFNMiwjmnl19+WfXr15dzTsnJyerZs6fvWACAIoyJjwAAB3AobxGwbds23XrrrZoxY4batWunN998k1l3AQDeZWZmSpJq167tOQkAwDf2mEa5lJQUXXbZZZo9e7YGDx6sKVOmUEoBAGGBiY8AAAdQTKNUbm6unnvuOTVq1EixsbFKTU1V3759ZWa+owEAICmvmHJ+KQBAophGpZ9//lktWrTQwIED1bFjRy1btoyBHwAQVg5MfMT5pQAAiWIadWbPnq1LL71UCxcu1Ouvv67x48frlFNO8R0LAIA/YOIjAEB+FNMosWfPHt17771KSEhQ+fLl9fHHH+uOO+7g0F0AQFhi4iMAQH4U0yhw4FDdIUOG6L777lNGRoYuueQS37EAADisjIwMVa9enYmPAACSKKYRLScnR88995yuuuoq7dy5U3PnztVLL72kUqVK+Y4GAMARZWZmchgvAOB3FNMItX79el1zzTUaOHCg2rVrp1WrVqlp06a+YwEA8Je2b9+ur7/+mmIKAPgdxTTCOOc0atQoXXrppVq+fLlGjx6tCRMmqFy5cr6jAQBQIAcmPmLGeADAARTTCLJhwwY1b95cPXv2VK1atbRy5UrddNNNTHAEAIgoGRkZkpj4CADw/1FMI0Bubq5ee+01XXTRRUpNTdWwYcO0YMECVatWzXc0AACOWmZmpqpXr87RPgCA3xX3HQBH9vXXX+v222/XwoULdd1112nEiBGqWrWq71gAABwzJj4CAByMPaZh6rffftOzzz6rmjVratmyZRoxYoTmzJlDKQUARLQDEx9xfikAID/2mIahRYsW6a677tKaNWvUsWNHvfLKK6pcubLvWAAAHLfMzExJYo8pAOAP2GMaRjZv3qxbb71VjRs3VlZWlj744ANNmjSJUgoAiBoHiikTHwEA8qOYhoH9+/dr2LBhuuCCC/TOO+9o4MCBWr16tVq0aOE7GgAAhYqJjwAAh8KhvJ7NmTNH999/vz777DM1btxYQ4YMUc2aNX3HAgAgKDIyMji/FADwJ+wx9WTNmjVq0aKFmjdvruzsbL3//vtasGABpRQAELW2bdumb775hvNLAQB/QjENsU2bNunuu+/WxRdfrNTUVL344otavXq12rVrJzPzHQ8AgKBZtmyZJCY+AgD8GYfyhsjmzZv1/PPPa+jQodq3b5969+6tJ554QnFxcb6jAQAQEszICwA4HIppkG3fvl0vvviiBg8erKysLHXv3l2PPfaYzj77bN/RAAAIqYyMDNWoUUNly5b1HQUAEGYopkHy448/asiQIRoyZIh27typ66+/Xv/+9791wQUX+I4GAIAXmZmZTHwEADgkzjEtZJ999pl69uypqlWr6plnntG1116rFStWaPz48ZRSAECRxcRHAIAjYY9pIXDOKTExUS+99JI++OADnXDCCerZs6f69eunc88913c8AAC8OzDxEXtMAQCHQjE9Dps3b9bo0aM1fPhwffHFF4qLi9OTTz6pu+66SxUqVPAdDwCAsJGRkSFJql27tuckAIBwRDE9Snv37tWHH36ocePGadq0adq3b5/q1aungQMH6vrrr9cJJ5zgOyIAAGEnMzOTiY8AAIdFMS2ArKwsLViwQFOmTNHkyZO1c+dOxcXFqU+fPurVq5dq1qzpOyIAAGEtMzNTV1xxhe8YAIAwRTE9BOecvvnmGy1YsEAzZszQvHnztGfPHpUpU0YdOnTQDTfcoCZNmqh4cZ4+AAD+ytatW/XNN9/ozjvv9B0FABCmgtqszKy5pMGSYiSNdM49F8ztHavs7Gx9+umnyszMVHJyshYuXKjvvvtOklSlShX16NFDrVu3VuPGjVWyZEnPaQEA+LO/GnPNrKSk0ZIul7RV0vXOufWhyHZg4iNm5AUAHE7QiqmZxUgaKuk6Sd9LSjez6c65z4K1zSPJysrSzz//rE2bNmndunX68ssv9dVXX2n16tX69NNPtX//fklSXFycGjdurAcffFCNGzfW3/72N5mZj8gAABRIAcfc2yVtd86dY2ZdJT0v6fpQ5MvMzJTExEcAgMML5h7TKyV95Zz7WpLMbLyktpKCWkz79u2rtLQ0ZWdnKzs7W3v37tW2bdu0a9euPyxnZqpSpYrOP/98DRgwQLVr19bll1+u6tWrU0QBAJGmIGNuW0n/DlyfJGmImZlzzgU7XGZmps4++2wmPgIAHFYwi2klSd/lu/29pL8fvJCZ9ZbUW8o7bPZ4lS5dWnFxcSpZsqRKlSqlkiVLqly5cqpYsaLi4uJ02mmn6eyzz1b16tVVqlSp494eAABhoCBj7u/LOOf2m9lOSeUlbTl4ZYU9NlerVk3VqlU77vUAAKKX99l7nHPDJQ2XpDp16hz3p7bPPReWp7ECABAxCntsHjRo0HFnAgBEt2JBXPdGSWflu105cB8AAChcBRlzf1/GzIpLOkV5kyABAOBdMItpuqRzzay6mZWQ1FXS9CBuDwCAoqogY+50SbcErneStCAU55cCAFAQQTuUN3D+yj2S5ihv6vo3nXOrg7U9AACKqsONuWb2pKQM59x0SaMkjTGzryRtU155BQAgLAT1HFPn3CxJs4K5DQAAcOgx1zn3WL7reyV1DnUuAAAKIpiH8gIAAAAA8JcopgAAAAAAryimAAAAAACvKKYAAAAAAK8opgAAAAAAryimAAAAAACvKKYAAAAAAK8opgAAAAAAryimAAAAAACvzDnnO8PvzGyzpA2FsKoKkrYUwnpCjdyhRe7Qi9Ts5A6twspd1TkXVwjrKdIYm8kdYuQOvUjNTu7QCvrYHFbFtLCYWYZzro7vHEeL3KFF7tCL1OzkDq1IzY0ji9R/V3KHFrlDL1Kzkzu0QpGbQ3kBAAAAAF5RTAEAAAAAXkVrMR3uO8AxIndokTv0IjU7uUMrUnPjyCL135XcoUXu0IvU7OQOraDnjspzTAEAAAAAkSNa95gCAAAAACJExBVTM2tuZp+b2Vdm9uAhHi9pZhMCjy81s2r5HhsYuP9zM4sPs9z3m9lnZrbSzOabWdV8j+WY2fLAZXqY5b7VzDbny9cz32O3mNmXgcstYZb75XyZvzCzHfke8/l8v2lmP5vZp4d53Mzs1cDfa6WZ1c73mM/n+69y3xjIu8rMUs3s0nyPrQ/cv9zMMkKXukC5G5vZznyvh8fyPXbE11gwFSD3P/Nl/jTwmi4XeMzn832WmSUGftetNrN/HGKZsHyN48gYmxmbCyk3Y3MhYmxmbC5g7vAZm51zEXORFCNpnaQakkpIWiHpwoOW6SPp9cD1rpImBK5fGFi+pKTqgfXEhFHuaySVDly/60DuwO1dYfx83yppyCF+tpykrwN/lg1cLxsuuQ9a/l5Jb/p+vgPbvlpSbUmfHubxFpI+lGSSrpK01PfzXcDc9Q7kkZRwIHfg9npJFcL0+W4saebxvsZCnfugZVtLWhAmz/cZkmoHrpeR9MUhfqeE5WucyxH/XRmbw+/5vlWMzYWZnbE5vJ7vxmJsLszcYTM2R9oe0yslfeWc+9o595uk8ZLaHrRMW0lvB65PknStmVng/vHOuWzn3DeSvgqsLyxyO+cSnXN7AjeXSKocomxHUpDn+3DiJc11zm1zzm2XNFdS8yDlPNjR5u4m6d2QJPsLzrkkSduOsEhbSaNdniWSTjWzM+T3+f7L3M651EAuKXxe3wV5vg/neP5vHLejzB1Or+9Nzrllgeu/SlojqdJBi4XlaxxHxNgcWozNIcbYHFqMzaEVTmNzpBXTSpK+y3f7e/35ift9Gefcfkk7JZUv4M8Gy9Fu+3blfSpxQCkzyzCzJWbWLgj5DqeguTsGdutPMrOzjvJng6HA2w4cllVd0oJ8d/t6vgvicH83n8/30Tr49e0kfWRmmWbW21OmI6lrZivM7EMzqxm4LyKebzMrrbwBYnK+u8Pi+ba8QzlrSVp60EPR8BovahibGZsLgrE5vH9vMTaHCGPz4RU/1h9EcJhZd0l1JDXKd3dV59xGM6shaYGZrXLOrfOT8E9mSHrXOZdtZnco7xPxJp4zHY2ukiY553Ly3RfOz3dEM7NrlDf4Nch3d4PA811R0lwzWxv41DEcLFPe62GXmbWQNFXSuX4jHZXWklKcc/k/wfX+fJvZScobkPs5534J5baBY8HYHHKMzSHE2BxyjM2HEWl7TDdKOivf7cqB+w65jJkVl3SKpK0F/NlgKdC2zayppIcltXHOZR+43zm3MfDn15IWKu+TjFD4y9zOua35so6UdHlBfzaIjmbbXXXQoRQen++CONzfzefzXSBmdonyXiNtnXNbD9yf7/n+WdL7Ct1hfH/JOfeLc25X4PosSbFmVkER8HwHHOn17eX5NrNY5Q187zjnphxikYh9jRdhjM2MzQXB2ByGv7cYm71gbD4c5+Ek22O9KG8P79fKO7zjwEnNNQ9a5m79cYKFiYHrNfXHCRa+VugmWChI7lrKO2H73IPuLyupZOB6BUlfKkQnchcw9xn5rreXtMT9/5OhvwnkLxu4Xi5ccgeWu0B5J5tbODzf+TJU0+FP+G+pP558/rHv57uAuaso79yxegfdf6KkMvmup0pqHka5Tz/w+lDeIPFt4Lkv0GvMV+7A46co71yXE8Pl+Q48d6MlvXKEZcL2Nc7lsP9mjM2MzYWSO7AcY3PocjM2hzB34HHG5iNlCeU/ViE9eS2UN1vUOkkPB+57UnmfZEpSKUnvBf6jfSypRr6ffTjwc59LSgiz3PMk/SRpeeAyPXB/PUmrAv+5Vkm6PcxyPytpdSBfoqQL8v3sbYF/h68k9Qin3IHb/5b03EE/5/v5flfSJkn7lHec/u2S7pR0Z+BxkzQ08PdaJalOmDzff5V7pKTt+V7fGYH7awSe6xWB19HDYZb7nnyv7yXKN3gf6jUWLrkDy9yqvEll8v+c7+e7gfLOo1mZ77XQIhJe41z+8t+WsTm8cjM2F25uxubwys3YXLi5w2ZsPvBpAwAAAAAAXkTaOaYAAAAAgChDMQUAAAAAeEUxBQAAAAB4RTEFAAAAAHhFMQUAAAAAeEUxBQAAAAB4RTEFAAAAAHhFMQUiiJlVNrPrfecAAAB5GJuBwkExBSLLtZJq+w4BAAB+x9gMFAJzzvnOAKAAzKyBpGmSdkj6VVIH59zXXkMBAFCEMTYDhYdiCkQQM5staYBz7lPfWQAAAGMzUFg4lBeILOdLWus7BAAA+B1jM1AIKKZAhDCzCpJ2Ouf2+84CAAAYm4HCRDEFIkc1ST/4DgEAAH5XTYzNQKGgmAKRY62kCmb2qZnV8x0GAAAwNgOFhcmPAAAAAABesccUAAAAAOAVxRQAAAAA4BXFFAAAAADgFcUUAAAAAOAVxRQAAAAA4BXFFAAAAADgFcUUAAAAAOAVxRQAAAAA4NX/A1wPBN2Kz0MEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(16, 6))\n",
    "axs[0].plot(t, X_num[:,0], label=\"numerical solution\", c=\"black\")\n",
    "axs[0].set_ylabel(\"$x$\")\n",
    "axs[0].set_xlabel(\"$t$\")\n",
    "axs[0].set_title(\"position\")\n",
    "\n",
    "axs[1].plot(t, X_num[:,1], label=\"numerical solution\", c=\"black\")\n",
    "axs[1].set_title(\"velocity\")\n",
    "axs[1].set_ylabel(\"$v$\")\n",
    "axs[1].set_xlabel(\"$t$\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 5
}