{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "49a9311e-18d1-489a-9b58-f90a1e40b333",
   "metadata": {},
   "source": [
    "# Practical Numerical Analysis Lecture 3: Inbuilt Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e7a65b0-dce2-4221-8191-ea5441b4c131",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import exp\n",
    "from numpy.polynomial.polynomial import polyroots\n",
    "from scipy.optimize import newton, root\n",
    "from scipy.integrate import solve_bvp, quad\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "489f3933-c7fb-4f8e-ae87-128dd68cf587",
   "metadata": {},
   "source": [
    "## 1D Rootfinding\n",
    "\n",
    "Use scipy.optimise.newton to find the real root of $f(x)=(x^2+2)(x-1)=x^3-x^2+2x-2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fd033a2d-5de1-4f0b-bd84-c98e9373570d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    \"\"\"Define f(x).\"\"\"\n",
    "    return x**3-x**2+2*x-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "58d3ea76-bf9a-4943-a12b-b923366d9a90",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fp(x):\n",
    "    \"\"\"Derivative of f(x).\"\"\"\n",
    "    return 3*x**2-2*x+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2e9afb84-265e-4e34-958b-61f8acdfca6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fpp(x):\n",
    "    \"\"\"Second derivative of f(x).\"\"\"\n",
    "    return 6*x-2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8bb0d27-ad14-4c6f-961b-3e8df4f1949e",
   "metadata": {},
   "source": [
    "Use the secant method (no derivatives), Newton's method (first derivative), and Halley's method (second derivative)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "01defb9d-bf71-4f5e-89a3-c9a2850a2de5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Secant root=1.0000000000000935, computed in 9 iterations\n",
      "Newton root=1.0, computed in 7 iterations\n",
      "Halley root=1.0, computed in 4 iterations\n"
     ]
    }
   ],
   "source": [
    "(s, rs) = newton(f, 3, full_output=True)\n",
    "(n, rn) = newton(f, 3, fprime=fp, full_output=True)\n",
    "(h, rh) = newton(f, 3, fprime=fp, fprime2=fpp, full_output=True)\n",
    "print(f\"Secant root={s}, computed in {rs.iterations} iterations\")\n",
    "print(f\"Newton root={n}, computed in {rn.iterations} iterations\")\n",
    "print(f\"Halley root={h}, computed in {rh.iterations} iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351cc8b2-d393-42c9-81ac-732bdd3ca4cf",
   "metadata": {},
   "source": [
    "## 1D Polynomial Rootfinding\n",
    "Use numpy.polynomial.polynomial.polyroots to find all roots of $f(x)$ using companion matrix approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ce7d9380-547f-4610-ab03-f0f1a2c299af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.60822483e-16-1.41421356j 3.60822483e-16+1.41421356j\n",
      " 1.00000000e+00+0.j        ]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([-2, 2, -1, 1])\n",
    "r1 = polyroots(a)\n",
    "print(r1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "764410ad-b1a4-4356-b2bc-88663bc84301",
   "metadata": {},
   "source": [
    "## Higher Dimensional Rootfinding.\n",
    "\n",
    "Find the intersection between the circle with centre at the origin and radius 1, and a straight line, i.e. find roots of\n",
    "\\begin{eqnarray*}\n",
    "g_1(x,y) & = & x^2+y^2-1=0 \\\\\n",
    "g_2(x,y) & = & x+0.5-y=0\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ed237563-8ba4-4b67-9a38-37fe9a397549",
   "metadata": {},
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    \"\"\"2D function defining cube roots of unity.\"\"\"\n",
    "    g1 = x[0]**2+x[1]**2-1\n",
    "    g2 = x[0]+0.5-x[1]\n",
    "    return [g1, g2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bccf6cd8-504b-4847-8267-4aa7258fc66e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Roots found: [0.41143783 0.91143783], [-0.91143783 -0.41143783]\n"
     ]
    }
   ],
   "source": [
    "r1 = root(g, [1, 1])\n",
    "r2 = root(g, [-1, -1])\n",
    "print(f\"Roots found: {r1.x}, {r2.x}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c403ae11-70eb-4807-931c-113de36e1a1a",
   "metadata": {},
   "source": [
    "## Integration\n",
    "\n",
    "Use scipy.integrate.quad to compute $\\int_0^1 f(x) {\\rm d}x$ (where $f(x)$ is the cubic above) which has exact value $-13/12$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "322b0521-3d8a-4032-8733-b1c50c5c33e9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical approximation I=-1.0833333333333333, error: 0.0\n"
     ]
    }
   ],
   "source": [
    "v1 = quad(f, 0, 1)\n",
    "# Exact value\n",
    "v = -13/12\n",
    "print(f\"Numerical approximation I={v1[0]}, error: {abs(v-v1[0])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30a2f124-cacc-4ea8-b4ad-3f09c0cf9065",
   "metadata": {},
   "source": [
    "## Solution of Elliptic BVPs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cff0dbbb-b12e-4ea9-a3dc-8df7e720b07e",
   "metadata": {},
   "source": [
    "Use scipy.integrate.solve_bvp to solve\n",
    "\\begin{eqnarray*}\n",
    "y^{\\prime \\prime} -3 y^{\\prime}+2y & = & 4, \\\\\n",
    "y(0) & = & 1, \\\\\n",
    "y(1) & = & 2,\n",
    "\\end{eqnarray*}\n",
    "which has exact solution\n",
    "\\begin{eqnarray*}\n",
    "y(x) & = & \\frac{{\\rm e}^{2x}-{\\rm e}^{1+x}}{{\\rm e}-1}+2.\n",
    "\\end{eqnarray*}\n",
    "\n",
    "We set $y=y_1$ and re-write the problem as\n",
    "\\begin{eqnarray*}\n",
    "y_1^{\\prime} & = & y_2, \\\\\n",
    "y_2^{\\prime} & = & 3y_2-2 y_1+4,\n",
    "\\end{eqnarray*}\n",
    "with boundary conditions\n",
    "\\begin{eqnarray*}\n",
    "y_1(0) & = & 1,\\\\\n",
    "y_1(1) & = & 2.\n",
    "\\end{eqnarray*}\n",
    "\n",
    "First define the right hand side function and boundary condition residual function, then solve and plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "240963f0-27c9-4d32-a182-79a5a9731111",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rhsfun(x, y):\n",
    "    \"\"\"Right hand side for BVP.\"\"\"\n",
    "    return np.vstack((y[1], 3*y[1]-2*y[0]+4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ee4b4c59-92a5-4348-af84-56f2a37b4503",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bcfun(ya, yb):\n",
    "    \"\"\"Boundary condition residuals for BVP.\"\"\"\n",
    "    return np.array([ya[0]-1, yb[0]-2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9c13a44e-290c-437c-b6bf-ff72bf8a835f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial mesh with 5 nodes\n",
    "x = np.linspace(0, 1, 5)\n",
    "# Initial guess for y\n",
    "yinit = np.zeros((2, x.size))\n",
    "# Solve!\n",
    "yapprox = solve_bvp(rhsfun, bcfun, x, yinit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "84153a12-1c01-4cce-a2e7-56a31127e27a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGyCAYAAAA/E2SwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAb/NJREFUeJzt3Xd4FNX+BvA3jQTSQ2hJqKEX6RKKgPSiovxUsCCgKNjgil69XEUvtqiIDcEOAjZEEFREBCUqXXqT3luoyRLSk/n9cThbkt3NzO5sfz/Pk2d2Z3dmTwLZffM9ZYIURVFARERE5CWCPd0AIiIiInMMJ0RERORVGE6IiIjIqzCcEBERkVdhOCEiIiKvwnBCREREXoXhhIiIiLwKwwkRERF5FYYTL/fUU0/hl19+8dnzu/s1PfH9EBGRvkI93QBHTZ48GXl5ebjpppvQs2dPi8dOnz6Nt956C0OGDMENN9zgmQbqZNq0aYiIiMCAAQN88vx6v2ZxcTH+85//oF+/fujXr5+u5yYi3yY/F6yJi4vDc8895+YWkaN8NpxMnz4d2dnZWLx4Mf755x+EhYUZHzt37hymTZuGlJQUnw8nU6dORefOnT3dDK9RXFyMadOmISoqymo44c+LKHBNnz4d1atXx0MPPVTusejoaA+0iBzls+EEAGrUqIFDhw5h5syZmDBhgqeb4xJPPfWUp5vgU/jzIgpsderU4fuAH/DpcJKWloa8vDy89NJLGDVqFGJjY20+9+zZs3jzzTdx9913o127dhaPvfjii2jSpAmGDRsGAMjPz8dzzz2HwYMHIy0tDfPnz8fhw4fRsWNH3HzzzcbnzJ8/HwcOHEDz5s0xbNgwhISElHtdg8GAJUuWYP/+/QgPD0ePHj0sqjnmr9W5c2csXLgQu3fvxtChQ9GhQwc89dRT6NOnT7luipycHPzwww/Yt28foqOjceONN6J9+/YW35PBYAAAhIWFISkpCQMGDECjRo00/pSF3NxcfP/99zh48CCqVKmCdu3aoXfv3uWed/LkSSxZsgSnTp1C1apVMXDgQDRv3tzuudX+25w/fx6vvvoqAODXX39FTk4OAKBt27a45557AMDmz6uidpn/O3Tt2hXffvst9u7di3r16uHuu+9GlSpVtP/QiMgr2XvfDQ4OxldffYWnnnoKiqLg22+/xalTp/C///0PVapUQWFhIZYsWYKdO3ciJCQEHTt2xIABAxAcbBrC+emnn+LKlSt44oknsHr1aqxatQq1atXCmDFjPPhd+xafHxA7depUXL58Ga+99prd5124cAHTpk3Dnj17yj02c+ZMLFu2zHg/Pz8f06ZNw++//45bbrkF27Ztg8FgwB133IGJEyfi8uXLGDhwIDZv3oyrV6/igQcewMiRI8udNyMjAw0aNMC0adNQUlKCc+fOYfDgwbj33ntRWlpq8Vq//fYbBg0ahDVr1uDChQvGdk6bNg2rV6+2OO9vv/2G+vXr46WXXkJOTg7Onz+Phx9+GM8//7zxOdWqVUPNmjVRs2ZNREREYNmyZWjevDnmzp2r/od7zfHjx9G4cWNMnToVxcXFyMrKwhtvvFGuW+XLL79Ew4YNMX/+fISGhmLNmjVo1aoVXnjhBbvnV/tvExISgho1agAAoqKijN+feSi19vNS0y7575CRkYEhQ4Zg/fr1KCoqwqRJk5CWlobCwkJtPzQi8lr23nf37NmDadOmYdmyZRg6dCjOnDmDDRs2IDc3F0eOHEHLli0xceJE5Obm4vz587jnnnvQrVs3XL582Xj+7777DnPmzMELL7yAqVOnIjc3t9z7ElVA8VGxsbHKkCFDFEVRlJEjRyoRERHK8ePHFUVRlK1btyoAlLffftv4/J07dyoAlHnz5pU7V40aNZSRI0ca71++fFkBoNSoUUPZt2+fcf+7776rBAcHK0OGDFH27Nlj3D9jxgwFgLJz507jvszMTCU2Nla55ZZblKKiIuP+LVu2KCEhIcoHH3xg8VqJiYnK9u3bjc/LyspSFEVRACjPPvuscf/p06eV6OhoZeDAgUp+fr7F97Fjxw67P7NnnnlGiYqKUq5evWrcV/b81jz77LNKTEyMkpeXZ7F/27ZtxtsHDx5UKlWqpNx1111KaWmpcf/LL7+sAFCWLl1q8zW1/Nvk5eUpAJQXXnjBalvLnlttu+S/Q82aNS3+bVevXq0AUD777DNbPx4i8hKxsbFKnTp1lCeffLLc11dffWV8nr333Xnz5ikAlBtuuEHJyclRFEVR8vPzlfz8fKVr165KjRo1lLNnzxqP2b17txIeHq6MGDHCuK9///5KTEyM8r///c/i3KSez1dOAODll18GIEZq6+nGG29E48aNjfcHDx6M0tJShIWFoVmzZsb9gwYNAgBs2LDBuG/u3LnIzs7GSy+9hNBQU+9Z27ZtceONN5arYHTp0gXXXXed8b6tLqo5c+bgypUrSE9PR3h4uMVjrVq1Mt5WFAWrVq3C66+/jmeeeQZPPfUUDh48iJycHOzdu1fLjwGFhYUoLCzEoUOHLPa3bt3aePurr75CYWEhJk+ejKCgIOP+J598EtHR0Zg1a5am19SL1nb16NHD4t+2a9euqFatGtavX++2NhPpISsrC19//bXbptbn5eXh888/L/d17Ngxt7y+FB4ebqyqmn9Ze0+19747YsQIREZGGs95/PhxrFmzBmPHjjVWcAGgefPmuPPOO/HNN98gNzfXuD8vLw9PPPGE1XNTxXx6zImUkpKCCRMmYOrUqRb/GZxlHkwAIDEx0er+atWqAQDOnDlj3Ldjxw4AwBdffIHg4GAoigJFUQAAmZmZyMzMtDiH+QeiPXv27EFoaChatGhh8zn5+fkYOHAgtm7diqFDh6JBgwaIiIgwjkG5dOmSqteSHnnkESxYsACtW7dGWloaevbsib59+6JHjx7G5+zfvx+hoaHlfjYRERFo2LAh9u/fr+k19aK1XU2aNCl3jho1auD06dMubSeRXgoKCvDoo4/i559/RmhoKBo0aOCWqfXZ2dkYPXo0hgwZgri4OOP+li1bom7dui5/fUnLgFh777tlH5PvFdbee1u2bImioiIcPXrUOJYtKSkJMTExaptNZfhFOAGASZMm4bPPPsPTTz+N119/vdzjcqpxcXFxucfkwMqyKleubHFfDni1tb/suUNCQlC9evVy573vvvsQERFhsc/8l9meoKAgi6BjzaxZs5CRkYG1a9daTKv94osv8Mknn6h6HXP16tXDP//8g+XLlyMjIwNLly7FK6+8gu7du2P58uWIiIiw267i4mKLqkVZjvzbqKW1XWX/bQHx72itbUTeqKioCGlpaXj33XcxYcIEHDx40OZzjx49io0bNyI2NhbXX3894uPjnX79l19+GS1btnT6PO5g73237GPyvUKOFzQn3x/M30/UvqeTdX4TTmJjY/Hcc8/hX//6F/r06VPu8eTkZACW1Q0AOHLkCK5evap7e9q0aYN58+Zh0KBBFc5W0XreuXPnYsuWLejUqZPV5+zfvx9hYWHl1vtwpmsiIiICQ4YMwZAhQwAAn3zyCR566CH88MMPuPPOO9G8eXOUlJRg9+7dFt09OTk5OHz4MG655Rab59bybyNHxNsLZ+acaReRL4qKilI1K+Q///kPPv74Y/To0QM5OTnYunUrvvzyS/Tv39+p11+9ejX27NmDhg0bom3btnb/MPEl8n18+/btGD58uMVj27dvR0REBBo0aOCJpvklvxhzIj3yyCNITU01Tjc1FxUVhXbt2mHBggUoKioCAJSUlOD111839ivqafTo0ahevTomTJhQ7q//8+fP488//3TovCNHjkS1atXw5JNPIjs727i/uLgYf/31FwDRNVFUVISNGzcaH9+yZYvFjCQtVq5cafFaAFCrVi0ApqrHvffeiypVqmDy5MnGny8ATJkyBbm5uRg7dqzN82v5t6lUqRKqVq2Ks2fPqmq7M+0i8lcLFizAvHnzsGfPHnz//fdYsWIFpk6dipEjR6KgoMDh8wYHB+OLL77AV199hX79+uGGG25Q/bvq7erVq4d+/frh448/xpEjR4z7N2zYgEWLFmH06NHlxgGS4/ymcgKID8pXX33VuF5JWW+++SYGDRqETp06oWPHjti6dSv++9//YvHixbq3JT4+HitWrMDw4cPRsGFD9OrVC7GxsTh06BAOHDhgNUCpkZCQgGXLluGOO+5As2bN0KtXL4SEhGD9+vUYO3YsbrjhBowePRpz5sxB3759cfvtt8NgMODw4cN44YUXrE55rsj+/fvxwAMPoEWLFmjQoAEuXryIpUuXYvjw4cbKQ0pKCubPn48RI0agbdu26Nq1K3bt2oVNmzbhnXfesRifYo2Wf5tRo0Zh5syZKCwsRNWqVS3WOSnL2XYR+aM5c+agSZMmWLlypbHbs6CgAJmZmdi7dy9at26NnJwcfPfdd3bPk5qaaly3qUqVKtiwYQM6dOgAQIxt69atGx599FEsXLjQ5d+TdODAAZtjTl599VVUqlTJ4XN//vnnGDJkCNq3b4/BgwejoKAAP/30EwYMGICpU6c6fF4qz2fDycsvv4zatWuX23/nnXfi/PnzyMvLQ/fu3S0eu/HGG7F//34sX74cpaWl+O9//4u6devizJkzqF+/vvF5lStXxtSpU8sdHx4ejqlTp6Jbt24W+8PCwjB16lR06dLFYv91112HXbt2Ye3atdi1axeCgoJwxx13oGvXrsaEbeu1JGvLsbdv3x779u1DRkYG9u3bh4SEBLz00kuoU6cOANEFs3btWixfvhyHDx9GSkoKBgwYgMzMTEydOtVicKia5d4feeQRjBw5En/99RcOHTqEuLg4TJkypdwg05tuuglHjx7FihUrcOrUKdxwww3o3bu3scpi7zXV/tvI4wcPHozdu3cjPz/fYhS8tXOraZe9f4eJEye6pLpG5CnHjx9HcHAwVq5cabF/5MiRxmpofn4+MjIy7J6nsLDQGE5iYmKMwQQQf0g9+uijePrpp1FaWmqxSJmrvPzyy8jPz7f5uOxisvf73r59e0ydOrXc+xYgKsYbNmzAn3/+iV27diEkJARPP/20xfcNAA8++KBLhgsEkiBFbec9ERH5lDFjxuDgwYPlQsYNN9yA1NRUfP755y59/dmzZ+P+++9HTk4OAz5p4ldjToiIqGKDBg3CokWLyk2RP3r0qMPnPHHiRLl9CxYsQMuWLRlMSDOf7dYhIiLrvvvuO+Tk5ODAgQPIzMw0VkhGjRoFAHjiiSfwyy+/oGPHjnjooYcQFxeHzZs3Y8eOHdi2bZtDr7l06VJ88cUXGDBgAKKjo7FkyRJs3brVJWP6yP+xW4eIyM889dRTuHDhQrn95t04JSUlWLhwIVavXo2QkBB06NABd955p3HMiSO2bduGxYsX4+LFi2jUqBHuueceVK1a1eHzUeBiOCEiIiKvwjEnRERE5FUYToiIiMir+MSA2NLSUpw+fRrR0dF+sxQy+TZFUXDlyhUkJSW5Zf0GoorwfZK8kaPvlT4RTk6fPm11wTUiTztx4gRSUlI83Qwivk+SV9P6XukT4SQ6OhqA+OZ4CWryBgaDAbVr1zb+3yTyNL5Pkjdy9L3SJ8KJLFHGxMTwl468Csvn5C34PkneTOt7JTvLiYiIyKswnBAREZFXYTghIiIir8JwQkRERF6F4YSIiIi8CsMJEREReRWGEyIiIvIqDCdERETkVRwKJwaDAZcvX9Z0TH5+Pk6cOIHCwkJHXpKIiIgChKZwsnDhQrRv3x61a9dGgwYN0LBhQ/z0008VHvfCCy8gISEBbdq0QWJiIqZPn+5wg4mIiMi/aQonq1atwqeffoqsrCxcunQJo0ePxu233479+/fbPGbevHmYOnUqVq5ciYsXL+KLL77AE088gRUrVjjdeCIiIvI/QYqiKI4eXFJSgsqVK2PmzJkYM2aM1ed06dIFqampmDdvnnHfjTfeiPj4eCxatEjV6xgMBsTGxiI7O5vXjCCvwP+T5G34f5K8kaP/L50aEHv48GEUFRUhOTnZ6uOlpaXYsmULOnfubLG/W7du2LRpkzMvTeRSa9cCr70GZGR4uiVE5MvWr+f7iCMcvipxcXExHnzwQbRu3Rp9+/a1+pwrV66goKAAVatWtdifmJiICxcu2Dx3QUEBCgoKjPcNBoOjzSRyyPLlwIsvAg89BPTs6enWEJEvKioC+vUDCgqACxeA6GhPt8h3OFQ5KS0txciRI3HgwAEsWrQIoaHWM05wsDh9UVGRxf7CwkKEhITYPH96ejpiY2ONX7Vr13akmUQOO3FCbPlfj4gcdfw4cOUKUFgInD7t3LlOngSWLgUcH4jhWzSHk9LSUowePRqrVq3CqlWr0KBBA5vPjY6ORmxsLM6ePWux/+zZs0hJSbF53KRJk5CdnW38OiE/KYjchOGEiJx15Ijptp3OAlU6dQJuuglYuNC58/gKTeGktLQU999/P1asWIFVq1ahcePG5Z5z6dIlHDH7F+nZsyeWL19u8Zxly5ahp51aeXh4OGJiYiy+iNyJ4YSInGUeTs6fd+5csvLyxRfOncdXaAonDz30EBYuXIhPP/0UISEhOHjwIA4ePIhLly4Zn/Pee++hbdu2xvvPPvss/vzzTzz//PP4+++/8eijj+LEiRN48skn9fsuiHSkKAwnROQ8PcOJ5Gz3kK/QFE7WrVuHGjVqYPz48RgwYIDx6wuzKJeQkGDR1dOxY0f8+uuv+PvvvzFq1CicOnUKf/zxBxo2bKjfd0Gko8uXgdxccdtO7yMR+amNG4GBA4H33nPuPAwnjtM0W2f37t0VPmf8+PEYP368xb4ePXqgR48e2lpG5CGyapKYCFSu7Nm2EJF7zZ0LjBwpbu/aBZT5ONNEr3BiftWXU6dEdTcoyPHz+QJe+I+oDHbpEPme7duBv/5y/jwffGC6ffKkc+fSa0BsdrblfWcH1/oChhOiMhhOiHxLaalYT6RXL+DYMefOVfaD39Fr1V69Cpw7Z7rvTOWkbDj55x/Hz+UrGE6IymA4IfItZ86IIFBcDPz+u3PnKhtOLl927DxHj1re97ZwcuqU+Ll5K4YTojIYToh8y8GDptt//OH4eYqKgKwsy31mk1E1kV06cr1RbwonOTlisH+DBkBJiXPnchWGE6IyGE6IfIte4eTiRbENDgbq1RO3nQ0n110ntnqOOXE2nKxbJ7b5+SKoeCOGE6IyGE6I3Oeff4Ayi4hrduCA6fbRo2LZeEfIAJGQAFSrJm47Gk7kYNp27cQ2L0+MQ3GEDCeyCuPsdOKtW0238/KcO5erMJwQmSktNb2pMJwQudaePUCbNsDgwc6dx7xyAjhePZHhJDFRBBTA8XBy5YrYJicDERHitqNdOzKc1Kolts4Gig0bTLfz8507l6swnBCZOX9ejM4PChJvKkTkOrNni9+3LVscryoApnDSooXYrl7t2HlcEU6io01VGGfDSc2aYutMOFEUU7eOs+dyJYYTIjOyHJyUBISFebYtRP6spAT46ivT/f37HTuPopjCyYABYutot4erwkliorjtDeHk5EnLWToMJ0Q+QIaTOnU82w4if/f775YhwtFBnpmZouoSHAy0by/2ORooXF05cXRQrMEgtnqEk/XrLe8znBD5AIYTIvf4+mvL+3v3OnYeWTWpU8c0JkPOutHK27t1atQQ2/x8UTFyxKZNlve9NZxourYOkb9jOCG9ZGVlYeXKlcjLy0O3bt1Qv359lxzjbufPi69mzZy7vou8VFvXrsCaNY5XTmQ4adgQqFpV3HY2nFSrpl84iYrSf8wJABQUmAbaalF2UTlvDSesnBCZYTghPWzatAkNGzbEW2+9hW+//RYtWrTAJ598ovsxntC3rxh42rFj+b/CtZBT9vv3F1tHKydyJdb69S0DRWmp9nPJ8GBeOXF0hVjzyomzQcdaOHE0VJQ9jrN1iHwAwwnpYfTo0ejfvz/Wrl2LH3/8EdOmTcPjjz+O03ZGajpyjLudOycusAcAmzcDU6Y4dp7CQtPaJn37iu3+/Y6tViqrHdWrmyonpaWmcRqOnEuPbh25uFl0NBAZKW47GihkOElMFGNrnDlX2eNYOSHyAQwn5Kxdu3Zh165dePTRR437Ro8ejdDQUHz//fe6HeMIR8cpSDKYSPv2OXaeU6dEW8LDgQ4dRPdEYaHlVXzVMg8UERFAlSriviOhQq9woiiWlZPKlcXt3Fzt5wJM4SQ21nQuhhOiAJGXZ7qKKMMJOWr3tcEUzZs3N+6LiIhAgwYNjI/pcUxBQQEMBoPFly0HDgB33w0MG6b527GwbZvYdu4stocPi+vRaCW7dFJSgNBQoEkTcd+Rrh05vkRO15WhwpFxJ9bCSVaW9opOXp6pWyk62hSY9AwnjnbHlD3O2XAyZgwwcKBYq0ZPDCdE18iVYSMjgfh4z7aFfJcMCbGxsRb74+PjbQYIR45JT09HbGys8au2nSWNS0rE7JiFC8tfLVcLWTkZOFB84JaUOHa+shXKxo3F9tAh7ecyDxSA44Nic3NN4SEx0fQeoCjlr21TEVk1AcT7iTPhpKjIdFxsrGkQrLOVExm+nA0nf/0F/PKL/tfoYTghusb8DdOZWQgU2Cpf+9M2p8y7tcFgQBX5KaXDMZMmTUJ2drbx64QsR1jRtKkY21FaCnzwgepvpRxZOWnXTsyOASyva6NW2etXySmysnKphQwnMpQ42h0jw0xYmKh2yK0j5zKfqRMcbAonjgQB82waE6Nft44MX86GE/PuKz0xnBBdw/EmpIdGjRoBAI6YDaAoKSnB8ePH0VB+outwTHh4OGJiYiy+7HnsMbH99FPHPpDy8kzdLm3amKodjqzsWjacODrVVlH0q5yYTyOWf5w4GnTKfmA7UzmRVZsqVURg0iucyO/N2dk6DCdELsZwQnro0KEDatWqha/M1mZftmwZLl++jJtuusm479tvv8Uf165Qp/YYZwweDNStKz5ov/tO+/G7d4tunMREcXmHa3nKocpJ2d81R8NJbq7pw7VsONEaKMyvSCw5G06iosTWmQGx5uNNzM/lDd06paWm7hz5veqFi7ARXXPsmNgynJAzQkJC8P7772P48OEwGAxITEzEjBkzMGHCBIsBry+++CLS0tLQo0cP1cc41y7g9tuBadOAv/8GRozQdrwcb9KmjagseEPlRFZHKlUyfTg6OiBWdp/ExZn2ORpOzKcRA/pUTsqGE2cHxOoRTsy/H1ZOiFxEhpO6dT3bDvJ9Q4cOxcaNGxEfH4+cnBx8/vnnePvtty2eM2zYMPTs2VPTMc6SOceRWTFyJdamTcVWz8pJ9epiq3XMiXmXjuyKcbRbR4YT894xedt8gKsatrp1nBlzItui14BYPcacyO/TfFyNXlg5IbqG4YT01KZNG7Rp08bm45MnT9Z8jLNksHBkqXj5+1GvntjKcHL8uPhrXO1S6jk5Ynou4HzlpOxgWPPbWqsd1sKJXDxN60wUPcecyGNkW7xpQKx595XekwhYOSGC6DuVpWZ265C/kuHk5Ent1QA5ZViGk2rVRFeDomibAix/z2JjTUFAhpOsLG3rppQdDAs43q0ju0/Mw4nsKnI2nMhAUVysfV0YGR5k+HMmnBQXiy9An24dVw2GBRhOiACIpbQLC0V5MjnZ060hco2EBFMXitaxIjKcyMpiUJApyJ86pf488rnmv2cJCaZl2WXgUKPsAmyA85UT86VmZLXi6lVt57JVOQG0V09keJChxJlwYn6MHrN1XDUYFmA4IQJg6gNPThbT9Yj8lSNdO/n5wJkz4rasnABi1g4AaLn8T2am2Mq1TQARTGSo0NK1Y61youeYE70qJ5UqmcKXXuHEkVBhfoye3TqsnBC5CMebUKBo1kxstQyKlV0xkZGW4zv0CieAY4Ni7XXrZGWZujDUsBdOnK2cBAU5PijWFZWT8HDnBulKDCdELsZwQoFCVk60hBPzLh3zgY/OhJOaNS33OzIo1tqAWPN1SuTAWzVcMSDWvLvD0UGxZcOJM7N1zM/l7MBagOGEyOUYTihQOBNOzLt0AKBWLbGVXT5qnD0rtmUrJ86EE/PKSWioadyIlq4dPbt1yq5zAji+EJsrKicREc5PSQYYTohcjuGEAoUMJ/v3q7/arq1wome3jiPhxNqAWMAUTuxcqLkca7N19BoQC+hXOdEjnOhVOeGAWCIXYzihQFG7tqguFBWpDxW2fj9cEU6cHXMCmAKGlqsJu3JALOB4OJGDWPUcEGseTpyZrcPKCZELKQrDCQWOkBDTFGBZEalIRZWTM2fE75EatsacyAGxaisn1i76JzlSObE2ldgV4cQbBsSWrZyo/bcri+GEyIWysky/ZFyAjQJB/fpia3YRZLtshRNZ/SgqUje+o7TUVBlxtlsnLw8oKBC3zQfBAtorJ4pif0CsHt06zo45kWNE9B4QW1qqfWE4ieGEyIVk1aR6df2vD0HkjWTIUBNOCgtN3TZlK4uVKplChZqunYsXTeNcZKVE0hpO5CJroaHlxzxorZzk5ooPacD7unVcNSBWnsfRcwHWB/7qheGEAp6tvwqJ/JWsnKjp1jl5UlQWwsPLBwpA27gT2aWTkFB+sUOtY05kOElIKH9dF62VExliyl7AzrxyorbrQ1GsDxT1pnBSubIIlvLn5mg4sTZlWi8MJxTwGE4o0Gjp1jG/grC1i7vJ6cRqwomtacSAadxIVpa6WUTm4aQsrZUT85k65t+j/NBVFPUf4Lm5piDjyjEnzg6IDQpyfjoxu3WIXIjhhAKN/L+upnJS0WBxRyonZQfDAqbl1BVF3eJply9bHmdOVk7UhhNr400AyyqK2q4d8/Ep5sd70zonegQdgOGEyKXKXtCMyN/JysmJExUPhjSvnFhjPmOnIramEQOim0F+yKkZXKumcqK1W6dsODHv5lE7KFY+r3Jl0/V0AO/r1nH2XADDCZFLsXJCgaZGDTGGpLTUdN0cW2Q4qahyoubKxPbCCaDtisL2womjlRPzacSS1kGxMnzI8SqSNy1fL8/hbDjhImxELiTL1gwnFCiCg9V37cjfD1uVEy1L2NsbcwKYgoa3VE4A7dOJ5fPKzvzTexE2R9Yn0bNyYj7wl5UTIp1lZZn6t9mtQ4FE7aDYirp1tIQTtZUTNeHEHWNOAP0rJ3oNiAXENG8t7AUdrcxnMDGcEOlM/lVYrVr5NxMif6ZmrRNFqbhbR4aTs2cr/kteVk7kMWVpCSdqunXUVk6sXVdH0nplYluVE0cGxJaWmhaasxZOnA06znQRyfEmZadf64XhhAIax5tQoEpNFdtDh2w/58IF0wdXSor158iZN2pWiZXhxNpsHUC/MSdapxKrqZxo7dbRY8yJ+SwaGSTCwkwDbT25FL75eBNrU8ydxXBCAY3hhAKVmnAiK4u1aokBtNZUqmQKFfa6doqLTQus2aqc6DXmRIaMK1dMK7/a485uHS3hxDw0yCDhzPoktgbEOjKV2JULsAEMJxTgOI2YAlWDBmJ7+LDt51Q03kRSM+7k/HnR7RMcbFoNtiy9xpzIyon5oE17vHVArAwToaHiS3K04qFn5cSV04gBhhMKcLK/XQ4OJAoUMpxcvGh7bIbaq3WrCSeyS6daNXFlZGv0GnMSHm5aHl/NuBN3TCV2JAiUDRNlz6W14qHngFiGEyIXYjihQBUdbbpWjq2uHbm/ot8PNeFEPmarSwdQP+akqMgUKKyFk6AgbeNO5HOsfdDqNSDWmcqJrXDiLWNOXIHhhAKWojCcUGCraNzJgQNi26iR/fNoqZzYGgwLqB9zYr68fVyc9edomU6spnLiiQGxrg4nzszWKTuLSG8MJxSwLl0ylSY55oQCkQwntsadqA0napawr2gaMaC+W0eON4mJsRyLYU7LQmzeOiC2bDeM5Gw4kaGkUiWxregSBtbINVbkOfTGcEIBSw6GrVnTdemfyJvZq5wUFprGnOhROZGP2aucyHBy9arpL3Nr7I03kRypnLhjQGxRkZi5pIatyokMF1rHnJQ9nwwWWhdzMz9Gju3RG8MJBSx26VCgk4NirYWTw4fFNNyoKPuBAtDWrWOvchIba1rDw964EzXhRG3lRFHcUzlxZPE0W+HE0VBRthIjg4U3Vk5sFMRsy8rKwty5c7FlyxaMHTsWnTt3rvCYbdu2YdGiRTh37hySkpIwbNgwNGnSxKEGE+mF4YQCnb3KiezSadiw4kW2zMOJolh/vprKSXCwmBp88aL4shVkZLeOHpWT/HxTJcOVK8TKaod8jppZLmW7YSS55oy96pK98/lCONFUOfn222/RvHlz7N+/H3PmzMEhe6v3XPPNN9+gQ4cOuHz5Mtq0aYPDhw+jZcuW+O233xxuNJEeGE4o0DVsKLYnTpTvIjh40PI59sgQkZtrOwyoqZwA6sadyMqJtTVOJLWVE9neoCDrl7DQa0BsUJD2KcB6Vk5KSkwhRI9uHXkurwgnHTt2xMGDB/H++++rPuajjz7CHXfcgenTp2PcuHH4/PPP0bVrV3z22WeaG0ukJ4YTCnTVq4swUFoK7N5t+ZjawbCAqBLIqoOtrh01s3UAddOJ1YQTtZUT82nEwVY+EfXq1gG0jxWxFU4cqZyYdyWZL4UP+EHlpH79+qii8Qo/devWxVn5vxJAcXExzp8/j3pcL5w8jOGEAl1QENCmjbi9fbvlY1rCCWCasXPqVPnHrlwxVRQqCidqphPLxxITbT9Ha+XEWpcOoF+3DqB9lk1F4URLxcPWdXoA58KJzw6Iffvtt5GSkoI2bdrgrrvuQqtWrdC7d29MnjzZ5jEFBQUwGAwWX0R6UhTTTASGEwpkMpxs22a5X2s4kUvcyyXvzcm/T6OiKl60S1ZOLlyw/Zzz58XWXjhRe2ViteFEbbeOOyonslqhpXIiXzMszLRCrx6zdbyicuKIbdu2YeXKlejcuTNuuOEGdOzYEUuWLMG+fftsHpOeno7Y2FjjV+3atV3dTAowZ8+KX9aQENtXWyVyVElJCV577TW0adMGTZo0wbhx43CpgmVP+/fvj5SUFIuvp59+2uVttRZODAZTyGjcWN155FpBMvSbUzMYVpLX3bEXTuRjtq7RA6hfIVZLOFEU++eSzzM/zpzWcGJrnRNHunWshQlv7tbRPFtHqwceeADDhw/H22+/DQB45JFHMHz4cDz66KNYs2aN1WMmTZqEiRMnGu8bDAYGFNKVXHSqTh3XlSUpcD377LOYNWsWZs+ejcTEREyYMAE33XQT1qxZgyAbU18yMzMxevRojB071rgv0tonnM5atxbb7dtNM23Wrxe369c3LXFfEXvhRO6r6AKCgClwyOqINTKcuLNbp7RUhIGyM2fKcke3jiMVDxlkzK8u7Uw4cfWAWJeGk5KSEhw/fhwtW7a02N+qVStkZGTYPC48PBzhtq7PTaQDGU7YpUN6MxgMeOedd/D+++9j8ODBAIDPP/8czZo1w4oVK9CvXz+bx8bGxiLFzaW8pk3FB4zBIBYmrF8fWL1aPNatm/rz2AsncsFDNb9vMpycO2f7OZ4IJ4AIHvbCiaJ474BYa5WOgOrWWbx4MR555BEAQEhICFq2bIkffvgBpaWlAICioiIsXboUbdu21fuliVST4UQuQkWklw0bNqCgoAB9+/Y17mvatCnq1q2LP//80+6x77//PlJTU9G9e3dMmzYNxWqXEnVCWBgg/36UXTsynHTtqv489sKJHHyuZh5ERZUTRVE35kSvcBIaavoArmjciXlFxFrlRO8xJ56snHjVgNg9e/Zg1KhRGDVqFADg448/xqhRozB79mzjc7Zt24avvvrKeP+jjz7Cpk2b0KpVKwwbNgzNmjVDZmYm3nnnHV2+ASJHMJyQq5w8eRIAULPMAIsaNWrglLWpLNe0bNkSU6dOxdKlS/Hwww9j2rRpGDFihM3n6zlxQHbtrFsnPqg2bBD3HamcnDghukDMycqJmnAiu5FshRPzpe3VjDlxNpwA6gfFmj+uZ7eOHouw+fWYk/j4ePTs2RMAjFsAaGw2YurWW2+1uN+pUyccPHgQGzZsQGZmJh5//HF06tQJYezoJw9iOCEtevbsiYNyVTIrEhISsGPHDgAwVolD5JSIa8LCwlBSUmLzHPPmzTOOR2natCliYmJw00034YUXXkDTpk3LPT89PR1TpkzR/L1YM3AgMHs2MG8ecOutomsiPh5o1kz9OZKSxADzoiIxADY52fSYI906tsKJ7NKJiLAeACQZTgoKxJetkQJqw8nlyxWHE9mlExFhmhFjztsqJ97craMpnNSqVctYNbGlTZs2aCOHf19TuXJlizBD5GkMJ6TF/PnzUWTnz0vzIFLt2qfrxYsXUaNGDeP+CxcuoFOnTjbPUXagbFpaGgDgn3/+sRpO9Jw4cOutYibN2bPAY4+JfV26WF+UzJbQUDHz7dgx8SXDSXGxaeaPlm6dvDwRBsqO3TAfb2JvWf3oaPG4oojqia2Bva6onNgKTd425iRgB8QSeaO8POD0aXGb4YTUMA8ZFenQoQOCgoKwZs0aDB06FICYiXPgwAFcf/31qs9z+FqCTrBxARk9Jw6EhQEPPQS8+CKwdavYd9992s9Tt64pnHTpIvadOiWWTq9UqeKl6wGxDkp4uPjgPX++fDhRM94EEMEqOlqED3eFE3uDYQH/mq3jcwNiibydLDHHxNi/cBiRI5KSknDbbbfhxRdfxPnz51FQUICnn34aSUlJuOWWW4zP69Onj3Edk02bNuGtt94yjhs5cuQIHn30UTRr1gxdtYxKdcJDD5k+rN55B7jzTu3nsDYoVv6+1a2rrhITFGS/a0fNGieSmoXYtIQTGT5s8bXKiR7dOl4xIJbIH5hPI67oaqtEjvj0009Rt25dJCUlIS4uDtu2bcNPP/2EymafMmfPnjUuzNasWTNcvHgRDRo0QEJCAlq2bImGDRtixYoVCA11T4E7ORn49Vfgl1+ACRMcO4e1cKJlpo5kbzqxmmnEkppBsWrCiQwbart1bFVOHF2EreyA2EConLBbhwIOx5uQq8XHx2PJkiW4evUqCgoKrHbN/Pbbb6h07Z09MjISr7zyCl555RVcunQJ8fHxNhdrcyVnhwbKcCKrJea3HQkn9ion7gwnnurW4QqxRAGE4YTcJTIy0uYqr7bGsdgaY+IL5OyebdtMq806coFNe+FE7ZgTwP3hRO9uHVvhxBtm67h6QCy7dSjgHDoktgwnRPpq21ZMoT171nR1Yvn7pqVyYm+tE1+unDg65iQQ1zlhOKGAI98sU1M92w4if1Olimm12Y0bxYfwpk3ifrt26s+j14DYisJJYaEpKLijciIrIN425qS4WN1FDc1xQCyRjkpLTd06DCdE+uvYUWz//lusNFtQINZQUXt1Y8B9Y06uXDHdjo62fR6t4aSiyomaMSeKYjuc6D1bB9BePWHlhEhHZ86IX/iQENPgPSLSj1zKZeNG4I8/xO0ePbTNjHPXmBPZpVOlilhEzhZPdOuYV0Vcvc4J4H3hhANiKaDILp06dVxXjiQKZDKcbNokFl8DtM8CshVOSkuBa7OvdQ0n9rp0AM9065hXV1w95gTQHk44IJZIRxxvQuRaLVqID2GDwVQ5cTScZGZajoU4d04ElOBg7w4nUVHWH9fSrSMDTFBQ+T+kZDjRs3KidcYOu3WIdMRwQuRaoaFA796m+zVqAE2aaDtHUpLY5uYCWVmm/fIaPUlJ6iqfasOJvfEmgPpwkpNj+fyytHTrmI83KdslJgOBs5WToCBTd5aj3TocEEukA4YTItebPx947z1xfZ0XXtC+EnNkpKl6Yr7arAwndeqoO09F4eTyZbGtaGkZvSsnWsJJ2fEmgGXlRO0sG2uVE8Cx6cQlJaYuO1ZOiHTAcELkelWqAI8/DqxZAzz8sGPnsLYUvt7hRI5fiY+3fx69KidaVoi1tcYJ4NgsG1vdMI4MrjV/TYYTIh0wnBD5BneEE1k50SucuKJyYi2cmFc/1Hbt6Fk5YTgh0lFWlukvJa4OS+TdrF2nRwYVtcsAyHCSn2+9MiDfD/Tq1nHVmJOyzAOB2oqHrcqJI+HE/DU55oTISbJqUr16xQPgiMiz5HL3zlROzGfhWKueuLtyIrt1CgrErCN77I05CQkRX/JcatiqnDjSrSOfGxxsaofeGE4oYBw8KLaNGnm2HURUMT26dUJCTEHBWjjRWjkpLBRLvduidoVYoOJQYW/MCaA9VLiicuKqLh2A4YQCCMMJke8oG05yc01L16sNJwAQFye25lOSJa2VE8B+9UR261Q05gSouGvHXrcOoH0hNj3HnDCcEOlIhpOGDT3bDiKqmAwnFy6IQHDihLgfHW0aS6JG1apie/Fi+cfUVk7Cw0UXBmA7nJhXVWxVTsLCTN0gFc3YqSicaK2cyHCi52wdhhMiHRw4ILYMJ0TeLy7OFEKOHbPs0tGybopcSVZWXcyprZwEBVU87kRWTQDb4QRQPyhW78qJDB+snBB5GVZOiHyLedeO1vEmkq3KiaKor5wAFYcTub9SJfszWNSGE1lZsTYgVr4O4HzlxJlw4srrkzGcUEC4ckVcpwPgGidEvsI8nMgpxVqvJm6rcpKba/pArqhyAqivnNirmgDqF2JzV+XEmdk6rJwQOUlOI05MNA2QIyLvJq/Js2YNsHKluN2qlbZz2AonsmoSFlZxoADUV05sDYaV9OrW8YbKCcMJkZPkeBPO1CHyHUOHiu3ChcD69WLsx223aTuHrW4d8/Emasaw6FU58YcxJxwQS6QTjjch8j1paWIxNtkF0q0bUKuWtnNUVDlRM94E0K9yorZbx11jTtitQ+RBMpxwvAmR7wgKAu66y3T/jju0n8NWOFE7U0dSG04CoXLCcEKkE9mt07ixZ9tBRNrcfbfYBgWZunm0qKhbR23lRF7y4soV649XtACb5IkxJ4rimhViXTlbJ9R1pybyHvv3iy3HnBD5lpYtgY8/BqpUAZKTtR9vXjlRFNP4Etmto7ZyUlE4UVs5kd007qycmAcPX5mtw3BCfs9gME0jZjgh8j0PPuj4sTKcFBSIACErG56unDg75kSGDDWhwjzA6FE54YBYIh3I8SbVq2tb9pqIfF+VKqYPcvOuHU9VTvTu1lFTOTEPMBxzQuQlOI2YKHAFBVkfFKt1QKzayom7u3W0VE6Cg03X9pG8tVuH4YT8HsebEAU2a4NitU4lVls50atbxxWVk7JVE8B7B8QynJDf40wdosDmjsqJq7p1Khpzoiac2FrjBGC3DpHHsFuHKLBZCyenToltzZrqzuHuqcSysqLHVGJ7lRNHunU4IJZIB7Jbh5UTosBUtlunoAA4c0bcVnshQb2nErvzwn+snBB5mUuXTH3LXLqeKDCVrZzIqklEhOmxisTEiK3BYP1xvacS67kIm6vGnHCdEyIHyapJcrKYUkjkLiUlJVi+fDnOnTuHESNGIKTsNAkrcnJykJGRgby8PHTp0gXJjqw6RuWUDSfHjoltnTrqLvoHeK5y4uoxJ87M1uGAWCIHyXAiL71O5A7Tp09HamoqnnzySYwePRpFKv4s3b59Oxo2bIjnn38eH3/8MRo3boy5c+e6obX+T44rOXlSbI8fF1u1XTqAKZzk5wPFxeUfV1s5kX8k2QsniuK+MSfeWjlhOCG/tm+f2DKckDvFxsZi9erVeOmll1QfM3LkSHTv3h1btmzBihUrkJ6ejnHjxuHs2bMubGlgkOPN5PuBDCd16qg/hwwngPXqiZ6VE/PA4Y1jTjgglshJ8s2Ig2HJne677z6kpKSofv6ePXuwfft2jB8/3rjvwQcfRFBQEL7//ntXNDGgyJl6Fy+KL/NuHbUqVTJ9GNsLJxVVTtSEE/OZPN44W4djToicxG4d8gW7du0CALRo0cK4r3LlymjQoIHxsbIKCgpQYPZns8HWSE1CZCSQkiK6dfbtc6xbBxDVk4sXy4eTkhJT2NCjciLDSVCQ7QDg77N1GE7Ib5WWcgE20sf333+P7Oxsm49HRERg+PDhDp9fnjsuLs5if0JCgs3XTU9Px5QpUxx+zUDTpEn5cKKlcgLYDie5uabbaisn5seUZT7exNaAXUcqJ3qHE1cOiGU4Ib91/Lj4C6RSJaBePU+3hnzZxo0bcUYujGFFdHS0U+Gk8rVPrKtXryLK7NPtypUrxsfKmjRpEiZOnGi8bzAYULt2bYfb4O+aNAF++w3Yu9e5cAKUDyfyfkiI9a4Tc1oqJ7a6dADHKifs1iHyArJLp2HD8he7ItIiPT3dpedPTU0FABw9ehQtW7YEAJSWluLEiRO48847rR4THh6O8Io+CclIdu2uXSuCQVCQ6OrRwlY4ycoS27i4iqcm6xVOPFk54YBYIidwMCx5syVLlmDt2rUAgI4dO6JGjRqYP3++8fEVK1bgwoULGDx4sKea6FdkOFm9Wmxr1qy4ylGWmnBSES3hxNYaJ4B+lROOOSFyMw6GJU9Zt24d9u3bh02bNgEA5s2bh7CwMPTt29e4sNqzzz6LtLQ0dOnSBaGhoXjnnXcwYsQIXL16FYmJiXjnnXfw8MMPo1WrVp78VvxG2T9StA6GBUyrxOoRToqKxEBaa1XditY4AfSrnHhrtw4rJ+S3WDkhT9m/fz8yMjKQk5ODkSNHYs2aNcjIyMAleS0FALfeeiu6du1qvD98+HCsXr0aQUFBOHnyJGbMmIGZM2d6ovl+qez4ksce034OWTkpOzFKjlmOja34HOYrVduqnrhzzAkHxBK52d69Ytu0qWfbQYFn5MiRGDlypN3nvPzyy+X2derUCZ06dXJVswJaSAjQuzewahXwySfAPfdoP4ce3TrmgSMvz/rsHr3Diatm67Bbh0ijnBzgxAlxm906RAQAP/wgLgSqdSCspEc4CQ4WwaKgwHblRO9uHV+creNQt05RURGOHj2Kq3JJPJWuXLnChYLILeR4k2rVTJdLJ6LAVqWK48EE0CecABWvdSLDib2LlXqyciKf6zUX/jt9+jSeffZZ1K9fH/Xr11e9rPLOnTvRvXt31KhRA40bN8aQIUNw7tw5hxpMpAa7dIhIb3qHE1uVExla7IUT84qHoth/PV9cIVZTOPntt99QuXJl4wh0NU6cOIEePXqgZcuWuHjxIs6ePYtx48Zh+/btmhtLpJYcDMtwQkR6sRVOtAyIBfQJJ7JyoijWr5JsTu9r67ijcqJpzMmIESM0v8Crr76K+Ph4vPfeewgNFS83cOBAzech0oKVEyLSm7sqJ3K/vXVOzKsWhYX2g4Les3W8rlvHEcuWLcOQIUMQGhqKzMxMFGn5CRA5iOGEiPTmTd065kGjonEnasJJcXHF3UOS/BgPdeGUGpeHk5MnT6KwsBAtWrRAy5YtER0djf/7v//D+fPnbR5TUFAAg8Fg8UWkVkmJaUAswwkR6UWPRdgAfcJJaKiY+QNU3CWjZhE2QH31xC8qJ0FBQZg9ezZmzZqF8+fP48iRI9i3bx/GjRtn85j09HTExsYav3gxK9JCXvAvPNyxVSCJiKyR4USGEUlrOJGhw5lwApiChR6VEyDAwklycjIGDhxoXFioVq1aGDduHJYvXw7FRg1p0qRJyM7ONn6dkAtWEKnwzz9i26gRL/hHRPpJTBTbK1csA4EMJ3oNiFUz5gRQP51Yz3BSWiq+yh6rN917jLKysmAwGFDn2lrBvXr1KhcusrKyEBUVhSAbl2/k1TbJGTKcNGvm2XYQkX+JixN/8JSUABcuAMnJokorP/zd2a0DqJ9po2adE0BdODF/jtdUTvLz83H06FEcPXoUAHDhwgUcPXoUFy9eND7nnXfewXXXXWe8/5///AcbNmzAG2+8gT179mDBggV4++237XbrEDmD4YSIXCE42FQ9kcMm5TTioCDTgNmKVLQIm9pwokflJCjIVGFWE07Mpy17TTjZvHkzevbsiZ49e6Ju3bp455130LNnT0ydOtX4nLi4OGPVBAAaN26MjIwMrFmzBkOHDsXMmTORnp6O559/Xr/vgsiMDCfNm3u2HUTkf6pVE1sZTsy7dIJVfqLq1a2jtXJiq0NCy3Rid1VONHXrdO3a1Vg1seVf//oX/vWvf1nsa9euHZYsWaK1bUSaKQqwZ4+4zcoJEemtenWxlYucax0MC+jXraO1cmJrRdewMNE95U3hxOUDYoncKTNTvFkEBwONG3u6NUTkb+xVTtRy95gTe906gGOVk5AQ0SXkKgwn5Fdkl079+vav6ElE5Ahb4cQVlRO9ZutUdC0cR8KJK6smAMMJ+RkOhiUiVyobTuSAWC3hpKJ1TtRclRjQZ0AsYFrpleGEyEU43oSIXMmdlRO9pxJX1K1T0QUEAYYTIoewckJEruTqcKIo+g6IVRR1A2IBVk6IXEZWTlq08Gw7iMg/lQ0nly+LrSMDYq2tc2IeNPSYSmweOPQcEMtwQqTSpUvA2bPiNisnROQKMpzIqcTHjoltSor6c9irnJgHFj0GxJoHF4YTIg/YvVts69RRv1IjEZEWMpxkZYkP6kOHxP3UVPXnUBNOwsIqDgBqKifmwYXdOkQeILt0uDIsEblKQoJpfY9z54AjR8RtvcNJRVUTQFvlJDjYNCunLIYTIheSlROONyEiVwkJMV1fZ+tW8WEdFqZft47aacSAKZyoqZzYqpoAnEpM5FIcDEtE7iC7dtavF9v69U0Xz1NDTeVETTiRgcNe5aSiNU4AxyontqowemE4Ib8hKyfs1iEiVyobTrR06QD2F2HTEk60dOuoCSdc54RIZ+YzdRhOiMiVZDjZsEFstYYTWTkpLARKSiwf0zLmRMuAWHvdOhxzQuQiskundm3O1CEi10pLE9ucHLF1NJwA4mrA5hwZc+KJbh2GEyIVdu0S25YtPdsOIvJ/99xjOcbEmXBSdiE2R8ac2KucaOnWYTgh0hnDCRG5S82awMCBpvtaw0lwsClYlB13ovdUYnbrEHkQwwkRudOoUabb9etrP15WRspWTrR06/hz5cTFk4GIXE9RGE6IyL1uvhn4v/8D6tVTV+UoKypKrDIrx61Ies/W8dV1ThhOyOdlZgIXL4pVG3lNHSJyh0qVgO++c/x4OXD/yhXL/a4KJ75WOWG3Dvk8WTVp2NCxv2CIiNxNhhNblRO9phL7arcOwwn5PBlOWrXybDuIpPz8fMyZMwfdu3dH06ZNUWjv0+Oa2267DU2bNrX4mjx5shtaS54QFSW2ZSsnrppKrGZArDctwsZuHfJ5HG9C3ubmm29GrVq10K1bN6Snp6O0tLTCY44cOYIhQ4Zg9OjRxn1xcXEubCV5kh7dOp6onMgAw3BCVIGdO8WW19Qhb/HTTz8hPDwc32kclFCtWjU0bdrURa0ib1JRONF7KjG7dYjcqLTUVDm57jrPtoVICrf3SWDHJ598gtatW2PQoEH46KOPVFVcyDfZCid6TyX21XVOWDkhn3b4sPhLIzxcDIgl8lUNGzbE4MGD0apVK2zduhX//e9/sWHDBsyaNcvq8wsKClBg9iezwWBwV1NJB3LMiaunEqvp1uFUYiKd7dghti1auP4S3hS4Bg0ahMOHD9t8PD4+HuvWrXPqNebPn4+Qa2uid+jQAYmJiRg6dCgmTZqERo0alXt+eno6pkyZ4tRrkufoOZVYTeXE17p1+HZOPk2ON+FMHXKlGTNmWFQpygrVIRmHmF+sBUD37t0BALt377YaTiZNmoSJEyca7xsMBtSuXdvpdpB76DHmRHbV+OPy9Qwn5NNk5YTjTciV6juyPrmTjh8/DgCIjY21+nh4eLjDY1vI89y1CBvXOSHyAFk5YTghX3PLLbcY1zHZunUrPv74Y+N6KGfPnsX48eORmpqKrl27erKZ5CK2wkl2ttjayKQWzAfEKor15/jqOicMJ+Szrl4FDh4Ut9mtQ97kf//7H5o2bYoJEyYAAFq3bo2mTZti9erVxuccPnwYZ86cAQCkpqZi9+7dqF69OurUqYO6desiPj4eK1asQCV7nyrks6wNiFUUbeFEVkMUxXawcFXlxNVj/NitQz5rzx7xS1m9OlCjhqdbQ2Ty8MMPY/jw4eX2m48J+fHHH1H52sCCmJgYvPvuu3jrrbdw6tQp1KpVC2Gu/tOUPMpa5SQnRyyPAGirnAAihFj7L8MBsURutn272LJqQt6mRo0aqFFBYrY2jiUkJAR16tRxVbPIi1gLJ7JqEhqqbcwJIEJIZGT55/jqgFh265DPkuGkTRuPNoOISDN74SQ2VlxlvSIhIabn2ZpO7KvrnDCckM+S4aR1a8+2g4hIKznm5OpVU1eOlvEmgAgmFc3Y8dVuHYYT8kmKYppGzHBCRL5GVk4AEVAAICtLbNWGE6DicCIrJ+zWIXKDY8fEXxlhYQCvk0ZEvqZyZSD42iew7NqRlRMtF6Ou6Po6rJwQuZHs0mne3P5fBERE3igoqPy4E63dOoD6bh1WTojcgONNiMjX6RFOKqqcaFnnhIuwETlp2zaxZTghIl9VdiE2OeZES7cOB8QSeRFOIyYiX+eOygm7dYjcJDsbkFevZzghIl/ljjEnXOeEyE1kl06dOkBCgkebQkTkMD3DCWfrEHnY1q1i27atZ9tBROQMOeZEhhNHxpzI7hquc0LkYQwnROQPZOVEDoh1pnKSn1/+sdJSU5hg5YTIxRhOiMgf6NGtc+3C1lYrJ+ZdPWoqJ8XFYvVtexhOiKzIzwf27BG3GU6IyJfpEU4iIsQ2L6/8Y+bhRE3lBABKSuy/nlwLheGEyMzu3eKXJyEBSEnxdGuIiBwnw4nBIN7XDAZxX8uYExlOrHXrmFdT1FROgIq7dlg5IbJiyxaxbdtW3SXFiYi8VY0aYnv2rKl6AjhWObEWTmTlJDTUdB0fa+RUYsB+OCkpMXX7mB/jCgwn5FNkOGnf3rPtICJyVp06Ynv8uKlLJzzcfhdMWXLMibVuHTXTiAH1lRPzx1g5ITKzebPYMpwQka+T4eTkSeDSJXFbS9UEsF85kYFFBhhbQkJMtxlOiDQqKgJ27BC327XzbFuIiJxVq5YIBkVFwL59Yp+W8SaAunAin2NLUJC66cReHU4MBgNmzpyJMWPGYP369ZqO3bVrF8aMGYMPP/xQ68sSYc8eUaaMjQVSUz3dGiIi54SGAsnJ4vbOnWKrtXIiqyLWwoncV1HlBNAeTrxqzMl3332Hpk2bYseOHfjss89w8OBB1cfm5uZi+PDhWLp0KVauXKm5oUQcDEtE/qZuXbGVVWFHu3WsjTlR260DaAsnoaGufw/WlH3atm2L/fv3IyoqCh999JGmFxo/fjz69OmjKdAQmeN4EyLyN3LcyapVYluvnrbj7XXryH0VdesAlgux2eKuacSAxspJamoqouTFADSYP38+NmzYgNdee03zsUSSDCccb0JE/kKGk6tXxbZLF23H6zEgFtBWOXFHOHFxrxFw+PBhjB8/HitWrECEmvgGoKCgAAVmq8cY5Mo0FLCKi01XI2blhIj8hQwnUufO2o63N5VY7YBYwDSGxFvCiUtn6xQVFeGuu+7Cf/7zH1x33XWqj0tPT0dsbKzxq3bt2i5sJfmC3bvFXwYxMUCjRp5uDRGRPszDSXw80LixtuPVdOv4YuXEpeFk+fLl2LZtm3GWzpgxY7Bjxw5s3rwZY8aMwcWLF60eN2nSJGRnZxu/Tpw44cpmkg/YtEls27e3v9IhEZEvMQ8naWna39/YreOAVq1aYcaMGRb7/v77b4SFhSEtLQ3hNpatCw8Pt/kYBSYZTjp29Gw7iIj0ZB5OtI43AdRNJdYyINZvw8kPP/yAlStX4r333kPdunUxZswYi8cXL16MiIiIcvuJ7Pn7b7Ht0MGz7SAi0lNMjOjOuXxZ+3gTwDNTib2uW+eff/4xds8AwGeffYYxY8Zgzpw5xuds2bIFc+fO1beVFNAKCkxrADCcEJG/eeUV4IEHgB49tB+rxwqxgPeFE02Vk9jYWKSlpQGAcQuIKcbSLbfcYnG/rPHjxyPEfCF/ogrs3Cl+KapW1b4GABGRt3v4YceP9dcBsZrCSVJSUoXdMe3atUM7OwtR9OvXT8tLEmHjRrHt0IErwxIRmTOfSqwolu+Rjkwl9slF2Ig8QYaT66/3bDuIiLyNefAoLLR8LGAqJ0SeIMNJp06ebQeRFqdPn8bOnTuRkJCA1q1bo1KlShUeU1BQgHXr1iEvLw+dOnVCQkKCG1pKvsw8nOTnA+YTXX15QCzDCXm17Gxg715xm9OIyRdkZmZi3Lhx2LRpE5o3b45jx44hPz8f33zzjcVYvbL++ecfDBgwAJUrV0Z8fDx27dqFWbNm4Y477nBj68nXVKokunIURYQT8wsH+vJUYnbrkFfbtEn80tWrB1Sv7unWEFXsypUrGDlyJI4fP47ly5djz5496N69O4YNG2b3uJEjR6JNmzbYs2cP1q1bh8mTJ2P06NE4f/68m1pOvigoyPZ0YldeldjVGE7Iq3G8Cfmahg0b4tZbb0XQtZGJwcHBuP3223H8+HFcunTJ6jH79u3D33//jSeeeALB15YIffTRR1FaWorvv//ebW0n32Rrxk7ATCUmcrcNG8SW403Il/3++++oWbOmzTEkO64t5GN+DbLIyEg0aNDA+FhZvEAqSbbCiasGxKoYPuU0hhPyWopiCiccb0Ke9Ouvv9r98A8PD8fNN99s9bGVK1dixowZ+OSTT2wen52dDQCIi4uz2F+1alVkZWVZPSY9PR1Tpkyx33AKCLaWsNe7W0fOBmI4oYB2/Dhw9qzo32zf3tOtoUC2dOlSnDp1yubjMTExVsPJ+vXrMXToUDz99NMYNWqUzePltcTy8vIQGRlp3J+Tk4MIGzX5SZMmYeLEicb7BoOBV3APULbGnGgZEOtt65wwnJDXWr9ebFu3BqpU8WxbKLC9++67mo/ZsGED+vfvj0ceeQSvvPKK3ec2aNAAAHD8+HE0a9YMAKAoCk6ePIn/+7//s3oML5BKUkVjTnyxcsIBseS1ZDixM/uSyCv9/fff6N+/P8aNG4fXXnvN6nNWrFiBzZs3AwCuv/56VK1aFd99953x8YyMDJw7dw4DBw50S5vJd3FALJEbyXDiyJU6iTzl4MGD6NevHxo1aoSOHTtaBI5+/fohJiYGAPDEE08gLS0Nn376KcLCwvDmm2/ioYceQlFRERITE/H6669j1KhRaNu2rae+FfIR5kvYm9N7QCzHnFDAKygAtmwRt1k5IV+SnZ2N3r17AwC++eYbi8c6depkDCf9+vVD48aNjY+NGjUKderUwfz583H06FG8/PLLGDlypPsaTj7LWuWkpMQUJrhCLJFOtm4Vv1iJicC17ngin9C+fXuLaoktb731Vrl9vXr1Qq9evVzRLPJj1sKJ2Sxz3bp1OOaEAt66dWKblsYrERMR2WMtnJh38agJJzJweEvlhOGEvNLatWLbpYtn20FE5O2sjTmRt0ND1S03LwNH2Ssbm2PlhAKaogBr1ojbXbt6ti1ERN7OWuVEy2BYwBQ47IUTd64Qy3BCXufYMeDMGZH2O3TwdGuIiLybvW4dNV06gLpuHRlc2K1DAUl26bRty8XXiIgqYq1bR2vlRE23DisnFNDYpUNEpJ69yome3TqsnFBA42BYIiL13NWtw8oJBazsbEBeIZ6VEyKiiukxIFbLbB1WTijgrFsHlJaKhdeSkjzdGiIi72dvKjErJ0Q6WL1abLt182w7iIh8hbumErNyQgHrr7/E9oYbPNsOIiJfoceAWM7WIbKhoADYuFHcZjghIlJHhhNXd+uwckIBafNmkfyrVQPMLtZKRER2yOoIV4glcgHZpdOtGy/2R0Sklh6VE87WIbLhjz/Etnt3z7aDiMiXREaK7dWrpn2OVk44W4fITHGxaaZOjx6ebQsRkS+JjRXb7GzTPlZOiHSwbRtw5Yr4JbvuOk+3hojId8hwUlhoqpiwckKkA9mlc8MNQEiIZ9tCRORLoqNNt2X1JDdXbLnOCZET/vxTbNmlQ0SkTXCwKaDIcCK3sqpSERk4SkrEKt3WsHJCAaWkxDRTh+GEiEg7GUIMBrHNyhLbuDh1x5sHDmtdO4oixgYCrJxQgNi+Hbh8WST/tm093RoiIt9TdlCsM+HEWteOeWBh5YQCwqpVYtu9OxAa6tm2EBH5IlvhJD5e3fHm1RBr4cR8HysnFBBkOOnVy7PtICLyVWXDyeXLYqu2chISYlr80lq3DisnFFCKi02DYW+80bNtISLyVebhpKDAtM6J2nASFGR/xo75PndUuBlOyKM2bxbrm8THA61be7o1RES+KSZGbA0Gy8XY5H417C3EJisnYWHuubwIwwl51O+/i22PHmI6HBERaWdeOZHjTWJjta0bZW8hNhlY3NGlAzCckIf99pvY9u7t2XYQEfky83CidbyJZK9bx7xy4g4MJ+QxeXmm6+n06ePZthAR+TJrlROt4UQGD1ZOKKCtWSMGbiUnA02aeLo1RES+S49wwsoJEYCVK8W2Tx/3DLAiIvJXcuCreThRu8aJZG9ALCsnFDDMwwkRETnOfPl6Zysn9tY5YeWE/NqFC8CWLeI2F18jInKOqwfEsnJCAeG338SFpFq2BJKSPN0aIiLfpueAWG8Yc8IrmZBHLF8utv37e7YdRK6wa9cufPjhh9i+fTsSEhJwyy23YNSoUQixs+jEsGHDsG/fPot9Q4cOxfPPP+/q5pIfkOEkNxc4f17c1jrmxJvWOWE4IbdTFODXX8VthhPyN5s2bcIDDzyAsWPHYtiwYTh8+DD+/e9/Y/PmzZg5c6bN4/bt24devXrhvvvuM+6rWrWqO5pMfsB8Jdjjx8XWl2frMJyQ2+3ZA5w6BUREAN26ebo1RPpq1aoVtm3bhqBrU9BuuOEGFBcXY9y4cXjrrbcQERFh89ikpCS0adPGTS0lfxIWBlSuLNaPOnZM7OM6J0QayC6dHj3ELxORPwkPDzcGE/N9iqKgtLTU7rGzZ89Gp06dMHToUMyZMweKoriyqeRnZNeOs906Pls5KS4uxtmzZ5GQkIAqVaqoOiY/Px+lpaWqn0/+a9kysR0wwLPtIHKH/Px8vPHGGxgwYIDd97969eqhf//+aNWqFbZu3YonnngC69evxwcffGD1+QUFBSgoKDDeNxgMuredfEtsLHD2rOm+L8/W0RROzp49iw8++ACzZ8/GiRMnMG/ePNx77712j/npp5+Qnp6O7du3Q1EU1K9fH2+//Tb69u3rVMPJN+XkAH/+KW4PHOjZthCpddttt+HIkSM2H4+Pj8eqVavK7S8tLcV9992Hy5cv45dffrH7GgsWLEDYtT9Lu3Xrhpo1a+LOO+/EU089hdTU1HLPT09Px5QpUzR+J+TPZOVE0rNbx6srJ8uXL0dQUBDWrVuHlJQUVccsWbIEb775Jq6//noAwJQpUzBkyBDs3LnT6i8c+bfffxcJvH59oHFjT7eGSJ309HTk5+fbfDw0tPxbaWlpKUaNGoU1a9YgIyMDSRXMmQ8r865/4403AoDN98pJkyZh4sSJxvsGgwG1a9e2+xrk35wNJz5bORk5cqTmF/jkk08s7j///PN47bXXkJGRwXASgGSXzqBBXLKefEfTpk01Pb+0tBT3338/Vq5ciYyMDDRq1Ejza546dQoAEB0dbfXx8PBwhIeHaz4v+a+WLYEVK0z3o6K0He9NY07cPiD2+PHjKCoqQs2aNd390uRhimIKJ+zSIX+lKArGjBmDFStWICMjA41tlAhvv/12vPjiiwCA7du3Y+7cuSgpKQEAXLx4EU888QTq1q2LbpzSRiq9+CLQtau4nZys/Q/AgJ2tU1JSgrFjx6JFixbo16+fzecVFBTAYDBYfJHv271bTHELDwd69vR0a4hc49dff8Xs2bMBAHfeeSfatGlj/Dp06JDxeXv37sXxawtS1K1bF2vWrEHVqlXRuHFjJCcnIzQ0FL/++iurI6RaVBSQkQHMmAHMmaP9eG+qnLhtnRP518TOnTvx559/lutfNceBXv7pp5/EtndvIDLSs20hcpXOnTtj69atVh9LTk423l64cKFx9k5cXBw++ugjvPfeezh27BiSk5MRyV8SckBoKPDII44d67NjThylKArGjh2Ln3/+GatWrbJZ5pQ40Ms//fij2N58s2fbQeRKMTExqhZSa9KkSbl94eHhFb4/ErmKmm4dn62cGAwG5OTkGEemK4qCcePGYcmSJfj999/RvHnzCs/BgV7+58IFYN06cXvwYM+2hYiIylPTreOVY04KCgpw8uRJnDx5EgBw+fJlnDx5Epfl9ZkBvPXWWxYB5LHHHsNXX32FuXPnIjY21ng8x5EElmXLxIDYNm0AFsGIiLyPzw6I/fvvv5GWloa0tDQkJyfj9ddfR1paGtLT043PiYmJsehXXb58OWJjY/HAAw8Yj01LS8Onn36q33dBXm/JErG96SbPtoOIiKzz2QGx3bp1M1ZNbJk4caLFeJGDBw861jLyG/n5gFwcc8gQz7aFiIis86YBsbzwH7ncb78BV6+Kefft23u6NUREZI03LV/PcEIut3ix2A4ZwlVhiYi8FSsnFDBKSkxTiG+91aNNISIiO7xpzAnDCbnUunVAZqa4IFWPHp5uDRER2eKzs3WItFq4UGxvvtl9/6mJiEg7Vk4oICgKsGiRuP1//+fZthARkX0ynLByQn5t0ybg+HFxHZ3+/T3dGiIiskdWRVg5Ib8mu3QGDQIqV/ZsW4iIyD7O1iG/pyjAt9+K23fc4dm2EBFRxex167ByQn5h0ybgyBGgShVe6I+IyBfY69Zh5YT8gqya3HyzCChEROTdOFuH/Jp5l86wYZ5tCxERqcN1TsivrVsnZulERQEDBni6NUREpAYrJ+TXvvxSbIcO5SwdIiJfYS+c5OeLbUSEe9rCcEK6Kioydencfbdn20JEROqZd+soiuVjDCfk01auBC5cAKpXB3r39nRriIhILfPxJMXFptuKAhQUiNsMJ+STvvhCbIcNA0JDPdsWIiJSzzycmA+KlcEEYDghH2QwAN9/L26PGOHZthARkTbmg13Nx53ILh2A4YR80HffAXl5QNOmQIcOnm4NERFpUVE4CQribB3yQXPniu1994n/xERE5DvMw4d5V475YFh3vbcznJAujhwB/vhD/Me9915Pt4aIiBwhl38w78px90wdgOGEdDJ7ttj26QPUru3ZthARkWNkAGE4IZ9XUgJ8/rm4/cADHm0KERE5QVZO8vJM+xhOyCetXAmcOAHExwNDhni6NURE5ChWTshvfPqp2N57r3v/8xIRkb5YOSG/cPYssHixuD1mjEebQkRETmLlhPzC7NlimeO0NOC66zzdGiIicgYrJ+TzSkuBTz4Rt8eO9WxbiIjIeTKAeDqc8Oon5LDly8X6JrGxwJ13ero1RN7l6tWr2L59O0JCQtC8eXNER0dXeExxcTE2bdqEvLw8dOjQQdUxRHrylnVOGE7IYe+/L7b33w9UqeLZthB5k8mTJ+Ozzz5Dw4YNYTAYcOTIEbz55pt48MEHbR5z4MABDBw4ECUlJYiLi8Phw4cxd+5cDOEUOHIjVk7Ipx08CCxbJm4/8ohn20LkbVJSUnDkyBGEh4cDAKZPn45x48bhtttuQ2JiotVj7rvvPjRu3Bg//vgjQkJC8PLLL2PEiBE4cuQIqlat6s7mUwDzlsoJx5yQQ2bOBBQFGDgQaNjQ060h8i5jx441BhMAaNu2LUpLS5GVlWX1+QcOHMD69evx1FNPISQkBAAwYcIEFBYW4nt5qW8iN/CWAbGsnJBmBgPw2Wfi9mOPebYtRN5q//792L9/P86cOYN3330X48ePR0MbSX779u0AgNatWxv3RUdHIzU11fhYWQUFBSgwuzqbwWDQsfUUqLxlKjHDCWk2a5YIKE2aAAMGeLo1RK6XkZGBnJwcm49XqlQJ/fr1s9i3adMmfPnllzh27BgKCgrQp08fm8fLikp8fLzF/qpVq9qstqSnp2PKlCnqvgEilVg5IZ9UUgK89564/cQTQDA7BikAzJ8/HydOnLD5eGxsbLlwcvfdd+Puu+8GAMyePRu33XYbtm/fjhYtWpQ7vlKlSgCAvLw8REZGGvfn5uYaHytr0qRJmDhxovG+wWBAbV51k5zEygn5pEWLxPThhARgxAhPt4bIPT744AOnjh81ahQefvhh/PHHH1bDSf369QEAJ0+eRJMmTYz7T548iVtvvdXqOcPDwy3GtRDpwVsqJ/y7l1RTFOD118Xtxx7j9GEiawwGA0pKSiz2HThwAAUFBUhKSjLu+/PPP7Fjxw4AQKdOnRAfH49FixYZH//rr7+QmZmJAew7JTdi5YR8zqpVwObNIllzICyRdfv378e4ceMwfPhw1KtXD8eOHcP06dPRvXt33HTTTcbnPfLII0hLS8Onn36KSpUq4bXXXsPjjz8ORVGQmJiIl19+GXfddRc6dOjgwe+GAg0rJ+RzXntNbO+/H6hWzbNtIfJWHTp0wPz583H58mV8++23OHz4MKZOnYpVq1YhNNT092CPHj0sZuc89NBDWLx4Mfbt24eVK1di0qRJmDt3rie+BQpgXISNfMqGDcCKFUBoKPDUU55uDZF3S01NxSuvvGL3OTNmzCi3b+DAgRg4cKCrmkVUIS7CRj7lpZfEdsQIoF49jzaFiIhcxFsqJwwnVKEtW4ClS8W04UmTPN0aIiJyFVZOyGc8/7zY3n030KiRZ9tCRESuw8oJ+YT160XVJCQEeOEFT7eGiIhciZUT8nqKAjz3nLg9ciQv8EdE5O9YOSGvt2IF8NtvQKVKwOTJnm4NERG5Gisn5NVKS4FnnhG3H3mEM3SIiAKBDCBFReJaagDDCXmRL78Etm0DoqOBZ5/1dGuIiMgdZOUEEKFEUYCCAnGfi7CRR129CvznP+L2f/8LJCZ6tj1EROQe5gEkL09MhrD2mKuxckLlTJ0KnD4tunL+9S9Pt4aIiNwlJAQICxO38/Mtx54wnJDHHDliuvLwG2+49z8jERF5nvnF/2Q4CQoyhRZ30Nytk5OTg6+++gpbtmzB/fffj+uvv77CYy5evIg5c+bg2LFjaNSoEUaNGoWoqCiHGkyuNWGC+M/Yqxdw++2ebg0REblbRARgMIhwIgNJRIQIKO6iqXKyaNEiNG7cGBs3bsRHH32E/fv3V3jM2bNn0bZtW/zwww+oWbMmPv/8c6SlpeHKlSsON5pc44cfgB9/FBf3e/999/5HJCIi72A+ndgTM3UAjeGkVatW2Lt3Lz799FPVx7z00kuIiorCr7/+ikmTJuH333/HuXPn8O6772puLLmOwSCmDAPAk08CzZp5tj1EROQZ5gux+UQ4adSoEWJiYjS9wA8//IDbb78dlSpVAgDExMTglltuwQ8//KDpPLYUFYnBm+Sc//4XOHUKaNDAdC0dIiIKPD5XOdEqPz8fJ0+eRP369S32169fHwcPHrR5XEFBAQwGg8WXNYoCPPYY0L49sGmTrk0PKBkZwIwZ4vZHHwFVqni0OURE5EE+VznRKu/a4vzR0dEW+2NiYpCbm2vzuPT0dMTGxhq/ateubfV5V64A69YBZ88C3bsDixfr1vSAceUKMHq0uP3QQ0CfPp5tDxEReZbfV06ioqIQHByMrKwsi/2XL1+22z00adIkZGdnG79OnDhh9XkxMcDq1UD//iLhDR0KvPmmqKiQOhMmAEePAnXrip8dEREFNr+vnISFhaFRo0bYs2ePxf49e/agRYsWNo8LDw9HTEyMxZctMTHATz8B48aJUPLvfwMPPGBabpdsW7AAmD1bzMqZO1csVU9ERIHNLysnP/30EyZOnGi8P3z4cMyfPx8XLlwAABw5cgRLly7F8OHDdXvN0FBg5kzg3XeB4GDxgdurF5CZqdtL+J3Dh4EHHxS3//tf0S1GRETkc5WTvXv3Yty4cRg3bhwAYPbs2Rg3bhy++OIL43M2bdqEWbNmGe8//fTTSE1NRfv27TFs2DB07twZvXr1wgMPPKDTtyAEBQHjxwM//wzExgJr1wIdOgB//63ry/iF/HzgjjuA7GygSxfghRc83SIiIvIW3lA50bRCbHR0NNq0aQMA+OCDD4z7zQes3nTTTahbt67xfpUqVbBq1SpkZGTg+PHjePzxx9GtWzcnm21b//7Axo3AkCHA3r1At25iJsqYMS57SZ+iKMDDDwNbtogL+s2f794liYmIyLuZV07keqmRke5tg6ZwkpycbKya2NKhQwd06NDBYl9ISAh69+6tvXUOatwY2LABGDlSzOB58EExq+f99y0vBx2Ipk8HPv9cdH99/TWQkuLpFhERkTcxv7aOHL9ZrZp72+C3F/6LiQEWLgRefVV8EM+aBaSlAfv2ebplnvPTT8ATT4jbb77JacNERFSerJzk5wPnzonbDCc6Cg4GJk0Cfv0VqF4d2LFDLNg2d66nW+Z+mzcDw4cDpaViNtO//uXpFhERkTeSlZPcXOD8eXGb4cQFevcGtm0DbrwRuHpVdPfcfTdQZvkVv7V3LzBggPje+/QBPviAF/UjIiLrqlYV2wsXGE5crlYtYMUK4KWXgJAQMd7iuuuA33/3dMtc6+BBEUguXBCzlxYt4gBYIiKyLTlZbE+fNoWT6tXd24aACSeACCXPPSdWlW3QADhxQlRVxo8XVQV/c+CAqBadOiWuMrxsGRdaIyIi+5KSxPbUKVZO3CotDdi+HRg7VtyfPh1o1UpUVvzF9u3ADTcAJ0+KYLJqlZg6TEREZI8MJ5mZpj/cGU7cJCoK+PBDYPlyoE4d4MgRoF8/YMQI319ZduVKseJrZibQurUIJjVqeLpVRETkC6pXFz0N8jp1YWFiBqw7BWw4kfr1A3btAh5/XAwS/eILoGlTUU0pLvZ067RRFLGM/8CBgMEgKicZGQwmRESkXnCwGKcpVavm/kkUAR9OADEO4733xMJt7dqJWTzjx4sBsz//7BtXOc7JAUaPBh59VISqu+8W3VRxcZ5uGRER+RrZtQO4v0sHYDix0LGjWPr+ww/F+Ix//gEGDxaDSteu9XTrbNu4UazfMmeOSLxTp4oKUHi4p1tGRES+SM7YATwTTjQtXx8IQkLEQNlhw4BXXhHdO3/8AXTtKq7b8+yzorvEG+TkAFOmAG+9JRZXS0kRoaRHD0+3jCiwZWdn44MPPsDatWsREhKCzp07Y9y4cYix03E/evRoHDhwwGLfzTffjGeeecbVzSUqx7xy4u5pxADDiU1xcaIC8fjjYm2U2bPF4Nnly4HOnYGJE8XFBT2xZkhJCTBvnghKp0+LfffcI7qmEhLc3x4istSzZ0/ccsstGDduHAwGA1577TV8++232LhxI4KDrRest27dim7dumH48OHGfTVr1nRXk4kseLpbh+GkAnXqAJ98IpbBf/11cdG8deuAO+4QZa/77wdGjRLrprhaYaFYPO6118Sqr4B43ffeE91PROQdVq9ejUizy7jWqVMHXbt2xf79+9G0aVObx9WrV8+lV20nUsvT3Tocc6JSgwbARx8Bx46JhdyqVRML1Lz0EpCaKqopb70FHD6s/2v/84+oktSpI4LQ3r1AfDzwxhvA7t0MJkTeJrLM9eXXrl2LhIQEJJu/41vx5Zdfonfv3rjvvvuwcOFCVzaRyC5WTnxMzZoikDz3HLB4sbja8cqVwPr14uvJJ4FGjYBevcTYlI4dgYYNxUBVtS5dEoNcf/8dWLoU2LPH9FitWsCECcC4cUBsrO7fHhHpZNasWZg1axbOnDmD0NBQZGRkINrOEs21atVCr1690KpVK2zduhWjR4/G2rVrMW3aNKvPLygoQIG8nj0Ag8Gg+/dAgcvT4SRIUbx/oqzBYEBsbCyys7PtDijzlDNngO++A77/Hvjrr/Lro1SuDDRpAtSvL8JFQgIQGSkCS3ExcOWKuCz18ePAvn1iWX1zoaHiwn0jR3punAtZ8vb/k6Sve+65B8eOHbP5eFxcHH766SeLfcePH8exY8dw7NgxTJ06FQkJCVixYgVCQ63/TZifn48Iea16iCrKiBEjcPjwYdSrV6/c8//3v/9hypQp5fbz/yTpIStLVOgB8bnmaG+jo++VDCc6MxjEiqx//CHGpmzdCpj9caNaw4Zildc+fUQwkf9JyDv40v9Jct6WLVuQm5tr8/GwsDB06tTJ5uNnzpxBcnIyvv32W9x+++2qXvP8+fOoXr06Fi9ejCFDhpR73FrlpHbt2vw/SbpQFDExxGAADh1yfFylo++V7NbRWUyMqG7I95LiYjEOZd8+MV4lM1N02+Tliem/oaFiKf1q1cQApEaNgObNGUaIvEm7du2cOr5atWoICwvDuXPnVB+Tee06GlWqVLH6eHh4OMK5mBG5SFCQmKV67Jh7JnyUxXDiYqGhQOPG4ouI/N+OHTtw9OhR3HLLLQCA0tJSvPHGGygtLUWfPn2Mz7vvvvvQokULPPPMM9i1axf279+PoUOHAhB/bf773/9GUlISZ++Qx1z77+gRnK1DRKSjpKQkfPPNN6hevTo6dOiAWrVqYfbs2ViwYAEam/2VsmXLFuOia7Vq1cLChQtRrVo1tGvXDikpKcjOzsYvv/yCypUre+pbIfIYVk6IiHSUmJiIr776CtnZ2Th06BCqVauG5OTkcouvzZs3zzh7p2rVqvjyyy+NxyQnJ6MGr9hJAYzhhIjIBWJjY+2OVWnbtq3mY4gCBbt1iIiIyKswnBAREZFXYTghIiIir8JwQkRERF6F4YSIiIi8CsMJEREReRWGEyIiIvIqDCdERETkVRhOiIiIyKswnBAREZFX8Ynl6xVFASCu1EnkDeT/Rfl/k8jT+D5J3sjR90qfCCdXrlwBANSuXdvDLSGydOXKFcTGxnq6GUR8nySvpvW9MkjxgT/9SktLcfr0aURHRyMoKMjiMYPBgNq1a+PEiROIiYnxUAt9A39W2tj7eSmKgitXriApKanc1WaJPMHe+yTA338t+LNSr6KflaPvlT5ROQkODkZKSord58TExPA/kUr8WWlj6+fFigl5EzXvkwB//7Xgz0o9ez8rR94r+ScfEREReRWGEyIiIvIqPh9OwsPD8cILLyA8PNzTTfF6/Flpw58X+RP+f1aPPyv1XPWz8okBsURERBQ4fL5yQkRERP6F4YSIiIi8CsMJEREReRWfWOfk8OHD2LFjB6pXr45OnTohJCTEJcf4g8zMTGzYsAGRkZHo1q2b3UFK+fn5WLx4cbn93bp1U7Vegq8rKirCypUrUVhYiCFDhqg6xmAwYM2aNQCArl27cg0E8hqlpaXYuHEjzp49ixYtWqBRo0YuOcZf7NixA4cOHUL9+vXRpk0bu889cOAANm/eXG7/8OHDXdQ675KZmYk//vgDjRo1Qtu2bVUdc/DgQezatQs1atRAp06dtC9WqXi5559/XqlSpYrSp08fJSUlRWnfvr1y4cIF3Y/xB3PmzFGqVKmidO/eXWnatKlSt25dZd++fTaff+bMGQWAMmDAAGXYsGHGr/Xr17ux1Z7x4osvKrVr11bq1aun1KhRQ9Uxq1atUuLj45UOHTooHTt2VOLj45VVq1a5tqFEKmRnZytdunRRkpKSlL59+yqRkZHKk08+qfsx/qC4uFgZNmyYEhcXp/Tr109JSEhQbrvtNqWwsNDmMdOnT1eio6Mt3ieHDRvmxlZ7xqlTp5Thw4crycnJSmxsrDJhwgRVxz3zzDNKZGSk0rdvXyUpKUlJS0tTsrKyNL22V4eTv/76SwGg/P7774qiKIrBYFCaNm2qjBkzRtdj/MGJEyeU8PBw5YMPPlAURfwC9uvXT+nWrZvNY2Q42blzp7ua6TXeeecd5fTp08rUqVNVhZP8/HwlKSnJ4pfz0UcfVZKTk5WCggIXtpSoYhMmTFBSU1OVy5cvK4qiKGvXrlWCgoKUX375Rddj/MEHH3ygxMTEKAcPHlQURVGOHDmixMXFKe+8847NY6ZPn640adLEXU30Gnv27FG++uorpaCgQGnfvr2qcLJixQolKChIWbNmjaIoipKVlaWkpqYqjz32mKbX9upwMnbsWKVt27YW+9566y0lKipKKSoq0u0Yf/D2228rMTExFh+UP//8swJAOXLkiNVjZDiZNWuWsmTJEmX37t1uaq33UBtOli1bpgBQjh49atx3+PBhBYDfv5mT90tMTFReeeUVi31dunRR7r33Xl2P8Qddu3ZVRowYYbHv/vvvVzp06GDzmOnTpyv16tVTfv75Z+XXX39Vzpw54+pmeh214WTkyJFK586dLfa9+uqrSnx8vFJaWqr69bx6QOzOnTvRsmVLi32tWrVCTk4Ojh49qtsx/mDnzp1o0qQJKlWqZNzXqlUrAMCuXbtsHhcUFIQZM2Zg5syZ6NKlC3r37o3z58+7vL2+ZufOnYiOjkbdunWN++rXr4/IyEjs3LnTgy2jQHfmzBlcuHDB6vuerf+bjhzjL2x9RlT0fZ8/fx5vvfUWXnjhBdStWxeTJ092ZTN9lq2f7+XLl3Hq1CnV5/HqAbHZ2dlISEiw2Fe1alUAQFZWlm7H+ANHvu8qVapg9erV6NKlCwDg3Llz6NatGx577DHMnz/fpe31NdZ+voD4Gfvz/yvyftnZ2QBg9fff3vuk1mP8gaIoMBgMVr/vgoIC5OfnIyIiotxxnTp1wtGjR5GYmAgAWL58OQYNGoTrrrsOd9xxh1va7isq+ixSO9nCqysn4eHhyMnJsdgn71v7D+ToMf7Ake87JibGGEwAoHr16nj88cexdOlSKFw42IK1ny8gfsb+/P+KvJ+ckWft99/e+6TWY/xBUFAQKlWqZPX7lo9Z07FjR2MwAYD+/fujW7du+PHHH13aXl+k12ewV4eT1NRUHD9+3GLfsWPHEBwcjHr16ul2jD+w9X0DQIMGDVSfJyYmBlevXkVubq6u7fN1qampuHz5Mq5cuWLcZzAYkJWVpennS6S3lJQUVKpUyervv63/m44c4y9svVfWq1dP03TXmJgYdoFbYevnGxoaijp16qg+j1eHk0GDBuHPP//EuXPnjPvmz5+P7t27IyoqCoCYf/3NN98Yy5RqjvFHgwYNwokTJ7Bhwwbjvvnz5yMpKck4hz8vLw/ffPMNTp48CUD0O5e1aNEiNGvWDJGRkW5ptzdbuHAh9uzZAwDo06cPQkNDsWjRIuPjCxYsQGhoKPr06eOpJhIhLCwMffv2xYIFC4z7Ll26hJUrV2Lw4MHGfZs2bcLPP/+s6Rh/NGjQICxZsgSFhYUAxHpH33//vcX3vX//fnzzzTfG+2XfK8+dO4fVq1ejY8eO7mm0Fzt9+jS++eYbY3Vk0KBB+O2333Dp0iXjc7799lv06dPHZmXKKm3jdd2rsLBQ6dSpk3Ldddcp77//vjJq1CilcuXKyoYNG4zPWbFihcV0WDXH+Ku7775bSUlJUd5++23l3//+txIaGqrMnz/f+PiJEycUAMqCBQsURVGU9957T7nxxhuVqVOnKh9++KEycOBAJTo6WlmxYoWnvgW3ycjIUL7++mvlnnvuUWJjY5Wvv/5a+frrrxWDwWB8TmRkpPLSSy8Z77/44otKVFSU8sorryivvPKKEhUVZfE4kafs2LFDiYqKUu6++25lxowZSvv27ZXWrVsreXl5xueMHTvWYjqsmmP8UWZmplK7dm2lT58+ysyZM5X+/fsrtWrVUk6dOmV8zttvv62Yfzz27dtXeeCBB5QPP/xQeeONN5T69esr1113nXLp0iVPfAtuU1hYaHxvbNCggTJw4EDl66+/tpih+OOPPyoAlAMHDiiKIpZdaNeundKuXTtlxowZyj333KNERkYqW7du1fTaXl05CQsLw6pVqzBq1Chs3LgRVatWxebNm3H99dcbn1OzZk0MGzYMcXFxqo/xV/PmzcNLL72EHTt2IDc3F3/88QfuvPNO4+NVqlTBsGHDULt2bQDA448/jvT0dFy4cAGbNm1Cz549ceDAgYCoBKxbtw6LFy9GcXExBgwYgMWLF2Px4sUW3Ta33347WrRoYbw/efJkfPnllzh69CiOHj2KL7/8Es8995wnmk9koVWrVti6dStSUlKwYcMGDBs2DKtXr7bo4+/YsaNFdUDNMf6oevXq2LRpE7p164Z169ahU6dO2LJlC5KSkozPadKkCYYNG2a8v2zZMtx4443Yvn07Tp06hf/973/YvHkz4uPjPfEtuE1RUZHxvbFjx46IiYnB4sWLsWrVKuNzkpOTMWzYMERHRwMQY07+/PNP3HXXXdiwYQOSkpKwdevWClfhLStIUTjykYiIiLyHV1dOiIiIKPAwnBAREZFXYTghIiIir8JwQkRERF6F4YSIiIi8CsMJEREReRWGEyIiIvIqDCdERETkVRhOiIiIyKswnBAREZFXYTghIiIir8JwQkRERF7l/wGLlzQtuaIg5gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xplot = np.linspace(0, 1, 100)\n",
    "yplot = yapprox.sol(xplot)[0]\n",
    "yexact = (exp(2*xplot)-exp(1+xplot))/(exp(1)-1)+2\n",
    "\n",
    "fig, axes = plt.subplots(1, 2)\n",
    "axes[0].plot(xplot, yplot, 'b')\n",
    "axes[0].set_title('Numerical solution')\n",
    "axes[1].plot(xplot, yplot-yexact, 'b')\n",
    "fig.subplots_adjust(wspace=1);\n",
    "axes[1].set_title('Error');"
   ]
  }
 ],
 "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.14.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}