2023-06-06_Meeting_Fenics_solvers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e798b8a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fenics as fn\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a4e401e",
   "metadata": {},
   "source": [
    "**Number of elements**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9befc164",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ne = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84fdb8fa",
   "metadata": {},
   "source": [
    "**Mesh definition**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e9dbb0c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = fn.UnitIntervalMesh(Ne)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "006f0e54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0cefb95a60>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAABACAYAAAAkuq3OAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAFwElEQVR4nO3dT4hVdRjG8e/jP5goGmkscsymgsyKRJuiLMJsMWkLC1wVBNKiFkUrUVtU0MZoExElISJtskUiBJVEUkZm4uCf0cKYDM1xYVZDYIKNvi3ujcZxLp47c875de95PjBw77nHOe/DHZ57+J2ZoyICMzMr35TUA5iZVZUL2MwsERewmVkiLmAzs0RcwGZmiUxrZueurq7o6ekpaBQzs/bU399/OiJmjd3eVAH39PSwd+/e/KYyM6sAScfG2+4lCDOzRFzAZmaJuIDNzBJxAZuZJeICNjNLxAVsZpaIC9jMLBEXsJlZIi5gM7NEXMBmZom4gM3MEnEBm5kl0tTNeCZi274h3th+hJPDZ5nd2cHqvnk8vrC76MMm5czO3I6qlheKz1xoAW/bN8S6rQOc/fs8AEPDZ1m3dQCgbd84Z3bmdsxctbxQTmY1878i9/b2RjO3o3xg/Q6Ghs9esn3G1CksnNuZ+fu0kn3Hhzl3/sIl2525vVQtc9XyQuPM3Z0dfLN2aVPfS1J/RPSO3V7oGvDJccoXGDdUu2iUzZnbS9UyVy0vNM7WqNcmotAliNmdHeOeAXd3dvDhs/cXeehkGp31O3N7qVrmquWFxplnd3bkdoxCz4BX982jY/rUi7Z1TJ/K6r55RR42KWeuceb2UrW8UE7mQs+A/12ortKVU2d25nZUtbxQTuZCL8KZmVmii3BmZtaYC9jMLBEXsJlZIi5gM7NEXMBmZom4gM3MEnEBm5kl4gI2M0vEBWxmlogL2MwsERewmVkiLmAzs0SauhmPpF+BYxM8VhdweoL/tlU5czVULXPV8sLkM98YEbPGbmyqgCdD0t7x7gbUzpy5GqqWuWp5objMXoIwM0vEBWxmlkiZBfxeicf6v3Dmaqha5qrlhYIyl7YGbGZmF/MShJlZIi5gM7NEci9gSY9KOiJpUNLacV6XpLfqrx+UtCjvGcqUIe9T9ZwHJe2StCDFnHm6XOZR+90j6byklWXOV4QsmSUtkbRf0mFJX5U9Y94y/GxfLeljSQfqmVelmDMvkjZJOiXpUIPX8++uiMjtC5gK/ATcDMwADgC3j9lnOfApIOA+4Ls8ZyjzK2PexcDM+uNlrZw3a+ZR++0APgFWpp67hPe5E/gemFt/fm3quUvI/BLwev3xLOB3YEbq2SeR+SFgEXCoweu5d1feZ8D3AoMRcTQizgFbgBVj9lkBvB81u4FOSdfnPEdZLps3InZFxB/1p7uBOSXPmLcs7zHAC8BHwKkyhytIlsxPAlsj4jhARLR67iyZA7hKkoArqRXwSLlj5icidlLL0Eju3ZV3AXcDv4x6fqK+rdl9WkWzWZ6h9gnayi6bWVI38ASwocS5ipTlfb4VmCnpS0n9kp4ubbpiZMn8NjAfOAkMAC9GxIVyxksi9+6aNqlxLqVxto39Pbcs+7SKzFkkPUytgB8sdKLiZcn8JrAmIs7XTo5aXpbM04C7gUeADuBbSbsj4seihytIlsx9wH5gKXAL8LmkryPiz4JnSyX37sq7gE8AN4x6Pofap2Oz+7SKTFkk3QVsBJZFxG8lzVaULJl7gS318u0ClksaiYhtpUyYv6w/16cj4gxwRtJOYAHQqgWcJfMqYH3UFkgHJf0M3AbsKWfE0uXfXTkvYk8DjgI38d/C/R1j9nmMixey96RefC8471xgEFicet6yMo/ZfzOtfxEuy/s8H/iivu8VwCHgztSzF5z5XeDV+uPrgCGgK/Xsk8zdQ+OLcLl3V65nwBExIul5YDu1q6ibIuKwpOfqr2+gdlV8ObVS+ovap2hLypj3ZeAa4J36GeFItPCdpDJmbitZMkfED5I+Aw4CF4CNETHurzO1gozv82vAZkkD1EppTUS07G0qJX0ALAG6JJ0AXgGmQ3Hd5T9FNjNLxH8JZ2aWiAvYzCwRF7CZWSIuYDOzRFzAZmaJuIDNzBJxAZuZJfIP/SG5Y35hu8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fn.plot(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ba2a229",
   "metadata": {},
   "source": [
    "**Finite element space**\n",
    "$$\\mathbb{U}_n := \\{\\phi_1,\\phi_2,\\dots,\\phi_n\\}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c725ecd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "Un = fn.FunctionSpace(mesh,'CG',1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dff8497",
   "metadata": {},
   "source": [
    "$$-\\Delta u = f \\text{   en } \\Omega$$\n",
    "$$u = 0 \\text{    sobre } \\partial\\Omega$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d40d5589",
   "metadata": {},
   "source": [
    "$$\\int \\nabla u \\cdot \\nabla v = \\int f v$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3840ba7d",
   "metadata": {},
   "source": [
    "**Trial and test functions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f1685745",
   "metadata": {},
   "outputs": [],
   "source": [
    "uh = fn.TrialFunction(Un)\n",
    "vh = fn.TestFunction(Un)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a507253",
   "metadata": {},
   "source": [
    "**Implementation of the variational formulation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "878cfe39",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = fn.inner(fn.grad(uh),fn.grad(vh))*(fn.dx)\n",
    "l = fn.Constant(2.)*vh*(fn.dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c9aa733",
   "metadata": {},
   "source": [
    "**Boundary condition**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6ece6dc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = fn.Constant(0.)\n",
    "bc = fn.DirichletBC (Un, g, 'on_boundary')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da2e8768",
   "metadata": {},
   "source": [
    "**Finite element function**\n",
    "$$u_n(x) = \\sum_{i=1}^{n} c_i \\phi_{i}(x)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9e51fb02",
   "metadata": {},
   "outputs": [],
   "source": [
    "un = fn.Function(Un)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d43a78ce",
   "metadata": {},
   "source": [
    "## Solver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad28afaf",
   "metadata": {},
   "source": [
    "**Option 1: Solve directly with a Fenics function**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "49334ccd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving linear variational problem.\n"
     ]
    }
   ],
   "source": [
    "fn.solve(b == l, un, bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88b3026a",
   "metadata": {},
   "source": [
    "Access to $\\{c_i\\}_{i=1}^{n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "002f452f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.16 0.24 0.24 0.16 0.  ]\n"
     ]
    }
   ],
   "source": [
    "print(un.vector()[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "053d26c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0ced2be100>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnVElEQVR4nO3deXxV9Z3/8dcnG5AQEpawhpAFRBEBIQRECNiOVbQVW5fiwq5Aq21/nem0/tpp++u004cz0047nbElyCa0ilq12nGrbS1BWULYd0zCFrawhrBk//7+SOikGMyF5N5zl/fz8eBh7r3nJJ8vhLeHk3vex5xziIhI+IryegAREfEvBb2ISJhT0IuIhDkFvYhImFPQi4iEOQW9iEiY8ynozexOM9ttZkVm9lQzrz9iZlsaf60ys6FNXttnZlvNbJOZFbbl8CIi0jJr6X30ZhYN7AFuB0qBdcBDzrkdTbYZA+x0zp02s4nA/3POjWp8bR+Q7Zw74Z8liIjIJ4nxYZscoMg5VwJgZsuBScBfg945t6rJ9muA1NYM1a1bN5eent6aTyEiElHWr19/wjmX0txrvgR9H+Bgk8elwKhP2H4W8HaTxw74g5k5IM85N7+lL5ienk5hoc7yiIj4ysz2X+k1X4Lemnmu2fM9ZnYbDUE/tsnTtzrnDptZd+A9M9vlnMtvZt/ZwGyAtLQ0H8YSERFf+PLD2FKgb5PHqcDhyzcysyHAAmCSc+7kpeedc4cb/1sGvEbDqaCPcc7Nd85lO+eyU1Ka/deHiIhcA1+Cfh0wwMwyzCwOmAy80XQDM0sDXgWmOOf2NHk+wcwSL30MfAbY1lbDi4hIy1o8deOcqzWzJ4F3gWhgkXNuu5nNbXx9HvA9oCvwSzMDqHXOZQM9gNcan4sBnnfOveOXlYiISLNafHulF7Kzs51+GCsi4jszW994gP0xujJWRCTMKehFRMKcL2+vFAkaF6vr+M3a/Zy9WOP1KAGVHB/Ho6P7ERejYzO5egp6CRk1dfU88fwG/ryrDGvu6o4w5hysP3CaX0y+meioCFu8tJqCXkJCfb3jH1/ezJ93lfHjz9/Ew6Mi66K6+fnF/PitXSR1iOVf7h2MRdr/6aRVFPQS9Jxz/PP/7OB3mw7zj3cMjLiQB5idm8XpCzX86i/FdImP4xt3DPR6JAkhCnoJer/4UxFLVu3jsbEZfHlCltfjeOabdwzkzIVq/vv9IpLjY3lsXKbXI0mIUNBLUHtu1T5+9sc93D8ile/cfUNEn7IwM350702cuVDDj97cSef4OO4b0aqiWIkQ+hG+BK3XNx3i+29s5/ZBPXj6CzdFdMhfEh1l/HzyMMb278Y3X9nCezuOeT2ShAAFvQSl93eV8Q8vbWZURhf+66GbiYnWt+ol7WKiyZsygsF9knji+Q2sKTnZ8k4S0fS3R4JO4b5TfOk367m+VyILpmXTPjba65GCTkK7GJZMH0lal3gee66QbYfKvR5JgpiCXoLKziNnmblkHb2TOrBkRg6J7WO9HilodU6IY9msHJI6xDJtUQElx895PZIEKQW9BI39J88zdVEBCe1iWDorh24d23k9UtDrldSBZbMabvEwZWEBR8ovejyRBCMFvQSFsrOVTFlYQG1dPctm5ZDaOd7rkUJGZkpHnpuZQ/nFGqYuLOD0+WqvR5Igo6AXz5VfqGHqogJOnKti8Ywc+ndP9HqkkDO4TxILpmWz/9QFpi9Zx/mqWq9HkiCioBdPXaiuZeZz6yg5fp75U7IZ1jfZ65FC1ujMrjzz8HC2HSpnzrL1VNXWeT2SBAkFvXimuraeL/16AxsPnOY/Jw9j7IBuXo8U8m4f1IN/vW8IHxSd4OsvbqKuPvhuLCSBpytjxRP19Y5/eHkzK/Yc5+kv3MTEm3p5PVLYuH9EKmcuVPOjN3eS1GErP/68LjaLdAp6CTjnHN9/Yzu/33yYb915PZNzIq+kzN8eG5fJ6QvVPPN+McnxcXzrzuu9Hkk8pKCXgPvZHz9i2Zr9zMnN5EsRXFLmb9/4zMC/Nl52jo9ldq5+ryOVgl4CavGHe/nFnz7iwexUnpqoo0x/MjN+OGkw5Rdr+PFbu0iOj+PB7L5ejyUeUNBLwLy2sZQf/H4Hd9zYQ+eNAyQ6yvjZg8M4e7GGp17ZQlKHWO64safXY0mA6V03EhB/2nmMb7y8hTFZXfnPySopC6S4mCjmPTqCoX2T+crzG1lVfMLrkSTA9LdN/K5g7ym+/JsN3Ni7E/OnqqTMCwntYlg8fSTp3eJ5/LlCtpSe8XokCSAFvfjV9sPlzFqyjj6dO7B4+kg6ttPZQq8kx8exdOYoOifEMX3xOorKVIIWKRT04jd7T5xn2qICEtvH8OtZo+iqkjLP9Uxqz7JZo4gymLpwLYfPqAQtEijoxS+OllcyZeFa6h0snTWK3skdvB5JGmV0S+C5mTlUVNYyZeFaTqkELewp6KXNnblQzdRFazl9vpolM0bSv3tHr0eSy9zYO4mF00dSevoi0xcXcE4laGFNQS9t6kJ1LTOWrGPfiQs8Oy2bIanJXo8kV5CT0YVfPjKc7YfPMntpIZU1KkELVwp6aTPVtfXMWbaezQfP8F8P38yYLJWUBbtP39CDnzwwhFXFJ/na8o3U1tV7PZL4gYJe2kRdvePrL21i5UcnePq+IbooJ4R8/uZUvv+5Qby7/RjfeW0bzqnxMtzovW7Sas45vvv6Nt7ccoRv33W9LrMPQTNuzeD0+Wp+8ecikhNi+b8Tb/B6JGlDPh3Rm9mdZrbbzIrM7KlmXn/EzLY0/lplZkN93VdC30//sIfn1x7gSxOyVJwVwr5++3VMGd2PvBUlzFtR7PU40oZaPKI3s2jgGeB2oBRYZ2ZvOOd2NNlsLzDeOXfazCYC84FRPu4rIWzByhL++/0iHsrpyzfvGOj1ONIKZsYP7rmRMxdrePrtXSR3iFWFdJjw5dRNDlDknCsBMLPlwCTgr2HtnFvVZPs1QKqv+0roemV9KT96cycTB/fkR/eqpCwcREUZP31gKGcv1vDt17aSHB/LnYN1U5hQ58upmz7AwSaPSxufu5JZwNvXuK+EiPd2HOObr2xhbP9u/HzyMKKjFPLhIi4mil89Opyb0zrz1Rc28WGRStBCnS9B39zf4GZ/LG9mt9EQ9N+6hn1nm1mhmRUeP37ch7HEK2tKTvLE8xsY3CeJvCkjaBejkrJwEx8Xw6JpI8nolsDspYVsPnjG65GkFXwJ+lKg6dsoUoHDl29kZkOABcAk59zJq9kXwDk33zmX7ZzLTklJ8WV28cC2Q+U89lwhaV3iWTJ9JAkqKQtbSfGxLJuVQ5eOcUxfXEBRWYXXI8k18iXo1wEDzCzDzOKAycAbTTcwszTgVWCKc27P1ewroaPk+DmmLSogqUNDAHROiPN6JPGz7p3a8+tZo4iOimLKwgIOqQQtJLUY9M65WuBJ4F1gJ/CSc267mc01s7mNm30P6Ar80sw2mVnhJ+3rh3WInx0pv8iUhQUALJuVQ68klZRFin5dE1g2K4dzVbVMWbCWk+eqvB5JrpIF41Vw2dnZrrCw0OsxpNHp89U8mLeaI+WVLJ89msF9krweSTywbt8ppixcS//uHXnh8dEkto/1eiRpwszWO+eym3tNFQjyic5X1TJ9yTr2n7rAgmnZCvkINjK9C796ZAS7jlTwuErQQoqCXq6oqraOOcvWs+1QOc88PJzRmV29Hkk8dtv13fnpg0NZU3KKr7ygErRQoaCXZtXVO77+4iY+KDrBv943hNsH9fB6JAkSk4b14Qf33Mh7O47x1KtbVYIWAvTeOPkY5xz/9LutvLX1KP909w3cPyK15Z0kokwbk87pC9X8/I8f0Tk+lm/fdYOujA5iCnr5mH9/dzcvFBzkiduyeGxcptfjSJD62qcHcPp8Nc+u3EvnhDi+PKG/1yPJFSjo5W88m1/CL/9SzMOj0vjGZ1RSJldmZnz/cw0laP/2zm6SO8Tx8CiVoAUjBb381UuFB/mXt3Zy95Be/HDSYP1TXFoUFWX8pLEE7Tu/ayhBu+smlaAFG/0wVgB4d/tRnnplC+MGdONnD6qkTHwXGx3FLx8ZwYi0znxt+UZWfqSuqmCjoBdWFZ/gKy9sZGjfZOY9OoK4GH1byNXpEBfNwukjyUrpyJxl69l44LTXI0kT+hsd4baWljN76XrSu8azWCVl0gpJHWJZOiuHlMR2zFiyjj3HVIIWLBT0Eaz4+DmmLS4gOT6WpTNHkRyvkjJpne6J7Vk2cxRx0VFMWbiWg6cueD2SoKCPWIfPXGTKgrVEGSybNYqeSe29HknCRFrXeJbOyuFidR1TFxVwvEIlaF5T0EegU+ermbJwLRWVtSyZkUNGtwSvR5Iwc33PTiyeMZIj5ReZvriAs5U1Xo8U0RT0EeZcVS3TFxdQevqiSsrEr0b068K8R0ew+2gFjz2nEjQvKegjSGVNHbOXFrL98FmeeXg4o1RSJn42YWB3/uOLw1i37xRPPr+BGpWgeUJBHyFq6+r52vKNrCo+yb/fP4S/U0mZBMg9Q3vzz5MG88edZXzrt1uor1cJWqDpvXQRwDnHd17bxrvbj/G9zw7iC8NVUiaBNWV0P06fr+Y/3ttDcnwc3/2sStACSUEfAZ5+ZxcvFh7kq5/qz8yxGV6PIxHqK5/qz+kL1Sz6cC9dEmJ58lMDvB4pYijow9y8FcXkrShhyuh+fP3267weRyKYmfHduwdx5kINP/nDHpLi45gyup/XY0UEBX0YW15wgKff3sXnhvbmB/fcqH8qi+eioox/u38IZy/W8L3Xt5HcIZbPDe3t9VhhTz+MDVPvbDvCt1/byvjrUvjpA0OJUkmZBInY6CieeWQ4I/t14e9f2sSKPSpB8zcFfRj6sOgEX31hE8P6JvOrR4erpEyCTvvYaBZMz6Z/90TmLlvP+v0qQfMnJUCY2XzwDLOXFpLRLYFF00cSH6ezcxKcOrWPZenMHHp0asfMJevYfVQlaP6ioA8jRWUVTF9cQJeOcSydlaOSMgl6KYntWDZrFO1jVYLmTwr6MHHozEWmLCwgOiqKZTNH0aOTSsokNPTtEs/SmaOoqq3n0YVrKauo9HqksKOgDwMnz1UxZcFazlXVsnRmDukqKZMQM7BnIotnjKTsbBXTFq2j/KJK0NqSgj7EVVTWMG1xAYfLL7Jo+kgG9e7k9Ugi12R4WmfypoygqKyCx55bx8VqlaC1FQV9CKurd8xeup5dRyr41SMjGJnexeuRRFol97oUfvbFYRTuP82Tz2/AOfXitAUFfQh7Z9tRVpec5If3Dua267t7PY5Im/jskN58e+IN/GlXGfkfnfB6nLCgoA9RzjnmrSgmo1sCD2b39XockTY1bUw6PTq1I29FsdejhAUFfYhaXXySrYfKeXxcJtG66lXCTFxMFLPGZrCq+CRbSs94PU7IU9CHqHn5JXTr2I4vDO/j9SgifvFQThqJ7WPIW1Hi9SghT0EfgnYcPkv+nuPMuDWd9rHRXo8j4heJ7WN5dHQ/3t52hH0nzns9TkjzKejN7E4z221mRWb2VDOvX29mq82sysy+cdlr+8xsq5ltMrPCtho8ks3PLyYhLppHR6niVcLbjDHpxERFseADHdW3RotBb2bRwDPARGAQ8JCZDbpss1PAV4GfXOHT3OacG+acy27NsAIHT13g91uO8FBOGknxsV6PI+JX3Tu15wvD+/ByYSknzlV5PU7I8uWIPgcocs6VOOeqgeXApKYbOOfKnHPrAF3O5mcLP9iLAbPG6U5REhkez82kuq6e51bt83qUkOVL0PcBDjZ5XNr4nK8c8AczW29ms6+0kZnNNrNCMys8flz91M05fb6aF9cdZNKwPvRK6uD1OCIBkZXSkc8M6sHS1fs5X1Xr9TghyZegb+69e1dzudqtzrnhNJz6ecLMcpvbyDk33zmX7ZzLTklJuYpPHzmWrdnPxZo6Zudmej2KSEDNGZ9F+cUaXlx3sOWN5WN8CfpSoOkVOanAYV+/gHPucON/y4DXaDgVJFfpYnUdS1bt41PXd2dgz0SvxxEJqOFpnclJ78LCD/ZSU1fv9Tghx5egXwcMMLMMM4sDJgNv+PLJzSzBzBIvfQx8Bth2rcNGst+uP8ip89XMHZ/l9Sginpg7IZNDZy7yP1t8Ps6URi3efsg5V2tmTwLvAtHAIufcdjOb2/j6PDPrCRQCnYB6M/s/NLxDpxvwWuNNqWOA551z7/hlJWGstq6eZ1fu5ea0ZEamd/Z6HBFPTLiuO9f16EjeihLuHdZHN7u/Cj7dZ8459xbw1mXPzWvy8VEaTulc7iwwtDUDCryz/SgHTl3g23fdoG9uiVhRUcbs3Cy+8fJmVuw5zoSBKvLzla6MDXKXyssyuyVw+6AeXo8j4ql7hvamV1J75qns7Koo6IPcquKTbDt0ltm5Ki8TuVR2tqbkFJsOnvF6nJChoA9y81YUk5LYjntvVnmZCMDkxrKz+fk6qveVgj6IbT9czsqPTqi8TKSJju1imDK6H29vO6qyMx8p6INY3ooSOraL4RGVl4n8jem3phMbHcX8lSo784WCPkgdPHWBN7ce4eFRaSR1UHmZSFPdE9tz3/BUfru+lOMVKjtriYI+SC38YC9RBjNuTfd6FJGg9Pi4DGpUduYTBX0QOnW+muXrDqi8TOQTZKZ05I5BPVm6ep/KzlqgoA9CS1fvo7KmnjkqLxP5RHPGZ3K2spYXCg54PUpQU9AHmYvVdTy3ah9/d0N3BvRQeZnIJ7k5rTOjMlR21hIFfZB5ef1BTl+oYY7Ky0R8Mnd8FkfKK/n9ZpWdXYmCPog0lJeVMDwtmex+Ki8T8cWEgSkM7JFI3ooSnLuaW2VEDgV9EHlr21EOnrrInPFZKi8T8ZGZMTs3k93HKvjLbt2drjkK+iDhnCNvRTGZKQncfoPKy0Suxj3DetNbZWdXpKAPEh8WnWT74bPMyc0kSuVlIlclNjqKmWMzWLv3FBsPnPZ6nKCjoA8SefkqLxNpjck5aXRqH8P8fNUiXE5BHwS2HWooL5t5awbtYlReJnItOraLYcot/Xhn+1FKjp/zepygoqAPAnn5jeVlo9O8HkUkpE0fk0FsdBTPrtzr9ShBRUHvsYOnLvDmlsM8MiqNTu1VXibSGimJ7bh/RCqvbCilrKLS63GChoLeYwtWlhAdZcy4NcPrUUTCwuPjMlV2dhkFvYdOnqvixcKD3DusDz2T2ns9jkhYyOiWwJ039mTZ6v2cU9kZoKD31NLV+xvKy8arvEykLc0dn8XZylqWq+wMUNB75kJ1LUtX7+PvbuhB/+4qLxNpS0P7JjM6s6HsrLpWZWcKeo+8XFjK6Qs1zNXRvIhfzFHZ2V8p6D1wqbxsRL/OZKd38XockbA04boUru+ZSF5+MfX1kV12pqD3wJtbj1B6+iJzVUUs4jdmxpzxmew5do6/7CnzehxPKegDrKG8rISslAQ+fX13r8cRCWufHXKp7CyyaxEU9AH2QdEJdhw5y5zcLJWXifhZbHQUs8ZlUrD3FBsiuOxMQR9g81YU06NTOybd3NvrUUQiwuSRfUnqEEteBFcYK+gDaGtpOR8WnVR5mUgAJbSLYeot/fjDjmMUR2jZmYI+gPLyi0lsF8NDo1ReJhJI08akExsdxYKVkXmuXkEfIAdOXuCtrUd4eLTKy0QCrVvHdjwwIpVX1h+KyLIzn4LezO40s91mVmRmTzXz+vVmttrMqszsG1ezb6R4dmUJMVFRzFR5mYgnHh+XSW19PYs/3Of1KAHXYtCbWTTwDDARGAQ8ZGaDLtvsFPBV4CfXsG/YO3muipcKD/L5m/vQo5PKy0S8kN4tgYmDe/HrNfupqKzxepyA8uWIPgcocs6VOOeqgeXApKYbOOfKnHPrgMt/91rcNxI8t3o/VbX1PJ6rugMRL83OzaSispblBQe9HiWgfAn6PkDT35XSxud80Zp9w8Kl8rLbB/Wgf/eOXo8jEtGG9k3mlsyuEVd25kvQN3dVj6/FET7va2azzazQzAqPHz/u46cPfi+uO8gZlZeJBI054zM5eraS1zcd8nqUgPEl6EuBvk0epwK+1sH5vK9zbr5zLts5l52SkuLjpw9uNXX1LFi5l5HpnRnRT+VlIsFgfGPZ2fz8kogpO/Ml6NcBA8wsw8zigMnAGz5+/tbsG/Le2nqEQ2cuMidX5WUiwcLMmDs+i4/KzvH+7sgoO2sx6J1ztcCTwLvATuAl59x2M5trZnMBzKynmZUCfw/8k5mVmlmnK+3rr8UEE+cc81aU0L97Rz6l8jKRoHL3kF70Se5AXoSUncX4spFz7i3grcuem9fk46M0nJbxad9IkP/RCXYeOcu/3T9E5WUiQSY2OopZYzP45//Zwfr9pxnRr7PXI/mVroz1k7zG8rJ7h0XUm4xEQsbknL4kx0dG2ZmC3g+2lJ5hVfFJZo3NIC5Gv8UiwSg+Loapo/vx3s5jFJWFd9mZUsgP8vJLGsrLclReJhLMpo5JJy4Cys4U9G1s/8nzvL31CI+M7keiystEglq3ju14IDuVVzccouxs+JadKejb2P+Wl6V7PYqI+OBS2dmiMC47U9C3oRPnqni5sJQvDO9Dd5WXiYSEfl0TmHhTL34TxmVnCvo2tHTVPqrrVF4mEmrm5GZSUVXLCwUHvB7FLxT0beR8VS3Prd7P7Tf0ICtF5WUioWRIajJjshrKzqpq67wep80p6NvIi+sOUn6xhrkTVHcgEormjs/i2NkqXt/ka5VX6FDQt4GaunoWfrCXnPQuDE8L7yvsRMLVuAHduKFXp7AsO1PQt4E3tzSWl6mKWCRkNZSdZVJUdo4/7wqvsjMFfSs1lJcVM6B7R24bqPIykVB2900NZWfzwqwWQUHfSiv2HGfX0QrmjM9SeZlIiIuJjuLxcRkU7j9N4b5TXo/TZhT0rZS3ooSendpzz9DeXo8iIm3gwZGNZWf54VOLoKBvhc0Hz7C6ROVlIuEkPi6Gqbek896OYxSVVXg9TptQOrVCXn4xie1jmJzTt+WNRSRkTLulH+1jo5gfJkf1CvprtO/Eed7edpQpKi8TCTtdO7bjwey+vLbxEMfCoOxMQX+Nnl1ZQmxUFNNVXiYSlh4bm0ldvWPRh3u9HqXVFPTX4HhFFS+vL+W+EX3onqjyMpFwlNY1nrtu6sXzaw5wNsTLzhT01+C5Vfuoqavn8XG6QEoknM0dn0VFVS3Prw3tsjMF/VU6X1XL0tX7uGNQTzJVXiYS1gb3SWJs/24sCvGyMwX9VVq+7iBnK2tVdyASIeaMz6SsoorXN4Zu2ZmC/irU1NWzcGUJORlduFnlZSIRYWz/bgzq1Yl5+cUhW3amoL8Kv998mMPllczV0bxIxDAz5ozPpOT4ef6485jX41wTBb2PnHPkrShhYI9ElZeJRJi7b+pFaucOIVuLoKD30V/2HGf3sQpm52ZipvIykUjSUHaWyfoQLTtT0Pto3l+K6ZXUns+pvEwkIj2QnUrn+NiQrDBW0Ptg44HTrN17SuVlIhHsUtnZH3eW8dGx0Co7U2r5YH5+CZ3axzA5J83rUUTEQ9PGpIdk2ZmCvgV7T5znne1HmXJLPzq2i/F6HBHxUJeEOL6Y3ZffbTrE0fLQKTtT0Ldgfn4JsdFRTBuT7vUoIhIEHhsXemVnCvpPUFZRySsbSrlveKrKy0QEgL5d4rl7SG+eX3uA8ouhUXbmU9Cb2Z1mttvMiszsqWZeNzP7RePrW8xseJPX9pnZVjPbZGaFbTm8v10qL5udqwukROR/zcnN5FwIlZ21GPRmFg08A0wEBgEPmdmgyzabCAxo/DUb+NVlr9/mnBvmnMtu/ciBca6qlmWr93PnjT3J6Jbg9TgiEkQG90li3IBuLPowNMrOfDmizwGKnHMlzrlqYDkw6bJtJgFLXYM1QLKZ9WrjWQNqecEBzlbW6mheRJo1JzeL4xVVvLbhkNejtMiXoO8DHGzyuLTxOV+3ccAfzGy9mc2+1kEDqbq2noUf7GWUystE5Apu7d+VG3t3Yn5+SdCXnfkS9M1d73/5qj5pm1udc8NpOL3zhJnlNvtFzGabWaGZFR4/ftyHsfzn95sPc6S8krkTsjydQ0SCl5kxd3wWJSfO816Ql535EvSlQN8mj1OBy4uZr7iNc+7Sf8uA12g4FfQxzrn5zrls51x2SkqKb9P7gXOOvPxiBvZIZMJ13s0hIsFv4uCe9O3SgXkrinEueI/qfQn6dcAAM8swszhgMvDGZdu8AUxtfPfNaKDcOXfEzBLMLBHAzBKAzwDb2nD+Nvf+7jL2HDvHnPEqLxORT3ap7GzjgTOs23fa63GuqMWgd87VAk8C7wI7gZecc9vNbK6ZzW3c7C2gBCgCngW+3Ph8D+ADM9sMFABvOufeaeM1tKl5K0rorfIyEfHRAyP60iUhjrwgLjvz6Zp+59xbNIR50+fmNfnYAU80s18JMLSVMwbMhgOnKdh7iu9+dhCx0bqWTERa1iEummm3pPOzP+5hz7EKruuR6PVIH6M0a2L+ihKSOsQyeWTfljcWEWk09ZZ+dIiNDtqyMwV9o+Lj53h3x1GmjO5HgsrLROQqdE6I44sj+/L6pkMcKb/o9Tgfo6BvtGClystE5NrNGptBvYNFHwRf2ZmCnsbysvWHeGBEKimJ7bweR0RCUN8u8Xx2SK+gLDtT0ANLPtxHTX09j49T3YGIXLvZuZmcr67jN2v3ez3K34j4oK+orGHZmv1MHNyTdJWXiUgr3Ni7sezsg31U1gRP2VnEB/3ygoNUVNYyJ1d1ByLSenPHZ3HiXBWvbQyesrOIDvpL5WW3ZHZlaN9kr8cRkTAwJqsrN/VJ4tn8EuqCpOwsooP+jc2HOXq2kjnjdW5eRNqGmTFnfGZD2dmO4Cg7i9igr6935K0o5vqeiYxXeZmItKE7b+xJWpf4oCk7i9igf393GR+VqbxMRNpeQ9lZBpsOnqFg7ymvx4ncoM9bUUKf5A58dojKy0Sk7T2Q3ZeuCXHkBUEtQkQG/fr9pynYd4pZYzNUXiYiftE+NpppY9L5864ydh+t8HSWiEy5vBXFJHWI5YsqLxMRP5oyuqHsLC/f2wrjiAv6orJzvLfzGFNvUXmZiPjXpbKzNzYd5vAZ78rOIi7oF6wsIU7lZSISII+Ny8DhbdlZRAV92dlKXt1wiAeyU+nWUeVlIuJ/qZ3j+dyQXrxQcIDyC96UnUVU0C/6cB+19fU8NlYXSIlI4MzOzeJ8dR2/9qjsLGKCvqKyht+s2c/Ewb1UXiYiATWodydyr0th8Yd7PSk7i5igf6HgABVVtao7EBFPzB2fyYlz1by6IfBlZxER9JfKy8ZkdWVIarLX44hIBLolsytDUpN4dmXgy84iIuh/t+kQx85WMWe8qohFxBtmxpzcLPaeOM8fth8N6NcO+6Cvr3fMzy/hhl6dyB3QzetxRCSC3Tm4J/26Br7sLOyD/s+7yigqO8dclZeJiMeio4zHx2WyubSctQEsOwv7oM/LL6ZPcgfuuqmX16OIiHD/iNSGsrMVgatFCOugX7//FOv2neaxcSovE5Hg0D42mulj0nl/93F2HT0bkK8Z1uk3b0UJyfEqLxOR4DLlln7Ex0Uzf0VgKozDNuiLys7x3o5jTL0lnfg4lZeJSPBIjo9j8sg03th8mEMBKDsL26B/Nr+EdjFRTLuln9ejiIh8zKwAlp2FZdAfO1vJaxsP8WB2X7qqvExEglCf5A7cM7Q3LxQc4MyFar9+rbAM+kUf7qW2vp7Hx6nuQESC1+zcTC5U1/HrNf4tOwu7oD9bWcPzaw5w1029SOsa7/U4IiJXdEOvTkwYmMKSVfv8WnbmU9Cb2Z1mttvMiszsqWZeNzP7RePrW8xsuK/7trUX1jaWl+Wq7kBEgt+c3CxOnKvmlQ2lfvsaLQa9mUUDzwATgUHAQ2Y26LLNJgIDGn/NBn51Ffu2maraOhZ+sJdb+3flptQkf30ZEZE2MzqzC0NTk3g2339lZ74c0ecARc65EudcNbAcmHTZNpOApa7BGiDZzHr5uG+beX3jYcoqqnQ0LyIhw8yYMz6LfScv8K6fys58Cfo+wMEmj0sbn/NlG1/2bRP19Y68/GIG9erEOJWXiUgIuePGnqR3jSfPT2VnvlxJ1FwT2OWTXGkbX/Zt+ARms2k47UNaWpoPY/2tCzV1jEzvwrgBKSovE5GQEh1lPHFbf7aUllNVW0/72Og2/fy+BH0p0LRDIBU47OM2cT7sC4Bzbj4wHyA7O/uq/5fWsV0MT9835Gp3ExEJCg9k9+WBbP/Utfhy6mYdMMDMMswsDpgMvHHZNm8AUxvffTMaKHfOHfFxXxER8aMWj+idc7Vm9iTwLhANLHLObTezuY2vzwPeAu4CioALwIxP2tcvKxERkWZZIO9y4qvs7GxXWFjo9RgiIiHDzNY757Kbey3srowVEZG/paAXEQlzCnoRkTCnoBcRCXMKehGRMBeU77oxs+PAtRY0dwNOtOE4oUBrDn+Rtl7Qmq9WP+dcSnMvBGXQt4aZFV7pLUbhSmsOf5G2XtCa25JO3YiIhDkFvYhImAvHoJ/v9QAe0JrDX6StF7TmNhN25+hFRORvheMRvYiINBGSQd+am5WHKh/W/EjjWreY2SozG+rFnG3J1xvLm9lIM6szs/sDOZ8/+LJmM5tgZpvMbLuZrQj0jG3Nh+/tJDP7vZltblzzDC/mbCtmtsjMysxs2xVeb/v8cs6F1C8a6o6LgUwabmyyGRh02TZ3AW/TcIer0cBar+cOwJrHAJ0bP54YCWtust2faajKvt/ruQPw55wM7ADSGh9393ruAKz528C/Nn6cApwC4ryevRVrzgWGA9uu8Hqb51coHtG35mbloarFNTvnVjnnTjc+XEPD3bxCma83lv8K8ApQFsjh/MSXNT8MvOqcOwDgnAv1dfuyZgckWsM9QjvSEPS1gR2z7Tjn8mlYw5W0eX6FYtC35mbloepq1zOLhiOCUNbims2sD/B5YF4A5/InX/6crwM6m9lfzGy9mU0N2HT+4cua/xu4gYbbkG4Fvuacqw/MeJ5o8/zy5Z6xwaY1NysPVVdzk/XbaAj6sX6dyP98WfPPgW855+rC5Ibwvqw5BhgBfBroAKw2szXOuT3+Hs5PfFnzHcAm4FNAFvCema10zp3182xeafP8CsWgb83NykOVT+sxsyHAAmCic+5kgGbzF1/WnA0sbwz5bsBdZlbrnPtdQCZse75+b59wzp0HzptZPjAUCNWg92XNM4CnXcMJ7CIz2wtcDxQEZsSAa/P8CsVTN625WXmoanHNZpYGvApMCeGju6ZaXLNzLsM5l+6cSwd+C3w5hEMefPvefh0YZ2YxZhYPjAJ2BnjOtuTLmg/Q8C8YzKwHMBAoCeiUgdXm+RVyR/SuFTcrD1U+rvl7QFfgl41HuLUuhAuhfFxzWPFlzc65nWb2DrAFqAcWOOeafZteKPDxz/mHwBIz20rDaY1vOedCttXSzF4AJgDdzKwU+D4QC/7LL10ZKyIS5kLx1I2IiFwFBb2ISJhT0IuIhDkFvYhImFPQi4iEOQW9iEiYU9CLiIQ5Bb2ISJj7/zuLmZsyv11JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fn.plot(un)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9699b812",
   "metadata": {},
   "source": [
    "**Option 2: Get the matrices from Fenics and solve using an external solver**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c06ada99",
   "metadata": {},
   "outputs": [],
   "source": [
    "B = fn.assemble(b)\n",
    "L = fn.assemble(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97211346",
   "metadata": {},
   "source": [
    "Apply boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "579ab1dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "bc.apply(B)\n",
    "bc.apply(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b25ec73d",
   "metadata": {},
   "source": [
    "Transform to a numpy matrix/vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "94ccef42",
   "metadata": {},
   "outputs": [],
   "source": [
    "B_np = fn.as_backend_type(B).mat()[:,:]\n",
    "L_np = L[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d328959",
   "metadata": {},
   "source": [
    "solve the linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "eca3afa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "vec = np.linalg.solve(B_np,L_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba9c82f",
   "metadata": {},
   "source": [
    "create a $\\mathbb{U}_n$ function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "64d61462",
   "metadata": {},
   "outputs": [],
   "source": [
    "un = fn.Function(Un)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afa4a947",
   "metadata": {},
   "source": [
    "assign the solution of the linear system to $\\{c_i\\}_{i=1}^{n}$ in \n",
    "$$u_n(x) = \\sum_{i=1}^{n} c_i \\phi_{i}(x)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5a4664d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-8.8817842e-17  1.6000000e-01  2.4000000e-01  2.4000000e-01\n",
      "  1.6000000e-01  0.0000000e+00]\n"
     ]
    }
   ],
   "source": [
    "un.vector()[:] = vec\n",
    "print(un.vector()[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0dff4ba1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0cecdde6a0>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnVElEQVR4nO3deXxV9Z3/8dcnG5AQEpawhpAFRBEBIQRECNiOVbQVW5fiwq5Aq21/nem0/tpp++u004cz0047nbElyCa0ilq12nGrbS1BWULYd0zCFrawhrBk//7+SOikGMyF5N5zl/fz8eBh7r3nJJ8vhLeHk3vex5xziIhI+IryegAREfEvBb2ISJhT0IuIhDkFvYhImFPQi4iEOQW9iEiY8ynozexOM9ttZkVm9lQzrz9iZlsaf60ys6FNXttnZlvNbJOZFbbl8CIi0jJr6X30ZhYN7AFuB0qBdcBDzrkdTbYZA+x0zp02s4nA/3POjWp8bR+Q7Zw74Z8liIjIJ4nxYZscoMg5VwJgZsuBScBfg945t6rJ9muA1NYM1a1bN5eent6aTyEiElHWr19/wjmX0txrvgR9H+Bgk8elwKhP2H4W8HaTxw74g5k5IM85N7+lL5ienk5hoc7yiIj4ysz2X+k1X4Lemnmu2fM9ZnYbDUE/tsnTtzrnDptZd+A9M9vlnMtvZt/ZwGyAtLQ0H8YSERFf+PLD2FKgb5PHqcDhyzcysyHAAmCSc+7kpeedc4cb/1sGvEbDqaCPcc7Nd85lO+eyU1Ka/deHiIhcA1+Cfh0wwMwyzCwOmAy80XQDM0sDXgWmOOf2NHk+wcwSL30MfAbY1lbDi4hIy1o8deOcqzWzJ4F3gWhgkXNuu5nNbXx9HvA9oCvwSzMDqHXOZQM9gNcan4sBnnfOveOXlYiISLNafHulF7Kzs51+GCsi4jszW994gP0xujJWRCTMKehFRMKcL2+vFAkaF6vr+M3a/Zy9WOP1KAGVHB/Ho6P7ERejYzO5egp6CRk1dfU88fwG/ryrDGvu6o4w5hysP3CaX0y+meioCFu8tJqCXkJCfb3jH1/ezJ93lfHjz9/Ew6Mi66K6+fnF/PitXSR1iOVf7h2MRdr/6aRVFPQS9Jxz/PP/7OB3mw7zj3cMjLiQB5idm8XpCzX86i/FdImP4xt3DPR6JAkhCnoJer/4UxFLVu3jsbEZfHlCltfjeOabdwzkzIVq/vv9IpLjY3lsXKbXI0mIUNBLUHtu1T5+9sc93D8ile/cfUNEn7IwM350702cuVDDj97cSef4OO4b0aqiWIkQ+hG+BK3XNx3i+29s5/ZBPXj6CzdFdMhfEh1l/HzyMMb278Y3X9nCezuOeT2ShAAFvQSl93eV8Q8vbWZURhf+66GbiYnWt+ol7WKiyZsygsF9knji+Q2sKTnZ8k4S0fS3R4JO4b5TfOk367m+VyILpmXTPjba65GCTkK7GJZMH0lal3gee66QbYfKvR5JgpiCXoLKziNnmblkHb2TOrBkRg6J7WO9HilodU6IY9msHJI6xDJtUQElx895PZIEKQW9BI39J88zdVEBCe1iWDorh24d23k9UtDrldSBZbMabvEwZWEBR8ovejyRBCMFvQSFsrOVTFlYQG1dPctm5ZDaOd7rkUJGZkpHnpuZQ/nFGqYuLOD0+WqvR5Igo6AXz5VfqGHqogJOnKti8Ywc+ndP9HqkkDO4TxILpmWz/9QFpi9Zx/mqWq9HkiCioBdPXaiuZeZz6yg5fp75U7IZ1jfZ65FC1ujMrjzz8HC2HSpnzrL1VNXWeT2SBAkFvXimuraeL/16AxsPnOY/Jw9j7IBuXo8U8m4f1IN/vW8IHxSd4OsvbqKuPvhuLCSBpytjxRP19Y5/eHkzK/Yc5+kv3MTEm3p5PVLYuH9EKmcuVPOjN3eS1GErP/68LjaLdAp6CTjnHN9/Yzu/33yYb915PZNzIq+kzN8eG5fJ6QvVPPN+McnxcXzrzuu9Hkk8pKCXgPvZHz9i2Zr9zMnN5EsRXFLmb9/4zMC/Nl52jo9ldq5+ryOVgl4CavGHe/nFnz7iwexUnpqoo0x/MjN+OGkw5Rdr+PFbu0iOj+PB7L5ejyUeUNBLwLy2sZQf/H4Hd9zYQ+eNAyQ6yvjZg8M4e7GGp17ZQlKHWO64safXY0mA6V03EhB/2nmMb7y8hTFZXfnPySopC6S4mCjmPTqCoX2T+crzG1lVfMLrkSTA9LdN/K5g7ym+/JsN3Ni7E/OnqqTMCwntYlg8fSTp3eJ5/LlCtpSe8XokCSAFvfjV9sPlzFqyjj6dO7B4+kg6ttPZQq8kx8exdOYoOifEMX3xOorKVIIWKRT04jd7T5xn2qICEtvH8OtZo+iqkjLP9Uxqz7JZo4gymLpwLYfPqAQtEijoxS+OllcyZeFa6h0snTWK3skdvB5JGmV0S+C5mTlUVNYyZeFaTqkELewp6KXNnblQzdRFazl9vpolM0bSv3tHr0eSy9zYO4mF00dSevoi0xcXcE4laGFNQS9t6kJ1LTOWrGPfiQs8Oy2bIanJXo8kV5CT0YVfPjKc7YfPMntpIZU1KkELVwp6aTPVtfXMWbaezQfP8F8P38yYLJWUBbtP39CDnzwwhFXFJ/na8o3U1tV7PZL4gYJe2kRdvePrL21i5UcnePq+IbooJ4R8/uZUvv+5Qby7/RjfeW0bzqnxMtzovW7Sas45vvv6Nt7ccoRv33W9LrMPQTNuzeD0+Wp+8ecikhNi+b8Tb/B6JGlDPh3Rm9mdZrbbzIrM7KlmXn/EzLY0/lplZkN93VdC30//sIfn1x7gSxOyVJwVwr5++3VMGd2PvBUlzFtR7PU40oZaPKI3s2jgGeB2oBRYZ2ZvOOd2NNlsLzDeOXfazCYC84FRPu4rIWzByhL++/0iHsrpyzfvGOj1ONIKZsYP7rmRMxdrePrtXSR3iFWFdJjw5dRNDlDknCsBMLPlwCTgr2HtnFvVZPs1QKqv+0roemV9KT96cycTB/fkR/eqpCwcREUZP31gKGcv1vDt17aSHB/LnYN1U5hQ58upmz7AwSaPSxufu5JZwNvXuK+EiPd2HOObr2xhbP9u/HzyMKKjFPLhIi4mil89Opyb0zrz1Rc28WGRStBCnS9B39zf4GZ/LG9mt9EQ9N+6hn1nm1mhmRUeP37ch7HEK2tKTvLE8xsY3CeJvCkjaBejkrJwEx8Xw6JpI8nolsDspYVsPnjG65GkFXwJ+lKg6dsoUoHDl29kZkOABcAk59zJq9kXwDk33zmX7ZzLTklJ8WV28cC2Q+U89lwhaV3iWTJ9JAkqKQtbSfGxLJuVQ5eOcUxfXEBRWYXXI8k18iXo1wEDzCzDzOKAycAbTTcwszTgVWCKc27P1ewroaPk+DmmLSogqUNDAHROiPN6JPGz7p3a8+tZo4iOimLKwgIOqQQtJLUY9M65WuBJ4F1gJ/CSc267mc01s7mNm30P6Ar80sw2mVnhJ+3rh3WInx0pv8iUhQUALJuVQ68klZRFin5dE1g2K4dzVbVMWbCWk+eqvB5JrpIF41Vw2dnZrrCw0OsxpNHp89U8mLeaI+WVLJ89msF9krweSTywbt8ppixcS//uHXnh8dEkto/1eiRpwszWO+eym3tNFQjyic5X1TJ9yTr2n7rAgmnZCvkINjK9C796ZAS7jlTwuErQQoqCXq6oqraOOcvWs+1QOc88PJzRmV29Hkk8dtv13fnpg0NZU3KKr7ygErRQoaCXZtXVO77+4iY+KDrBv943hNsH9fB6JAkSk4b14Qf33Mh7O47x1KtbVYIWAvTeOPkY5xz/9LutvLX1KP909w3cPyK15Z0kokwbk87pC9X8/I8f0Tk+lm/fdYOujA5iCnr5mH9/dzcvFBzkiduyeGxcptfjSJD62qcHcPp8Nc+u3EvnhDi+PKG/1yPJFSjo5W88m1/CL/9SzMOj0vjGZ1RSJldmZnz/cw0laP/2zm6SO8Tx8CiVoAUjBb381UuFB/mXt3Zy95Be/HDSYP1TXFoUFWX8pLEE7Tu/ayhBu+smlaAFG/0wVgB4d/tRnnplC+MGdONnD6qkTHwXGx3FLx8ZwYi0znxt+UZWfqSuqmCjoBdWFZ/gKy9sZGjfZOY9OoK4GH1byNXpEBfNwukjyUrpyJxl69l44LTXI0kT+hsd4baWljN76XrSu8azWCVl0gpJHWJZOiuHlMR2zFiyjj3HVIIWLBT0Eaz4+DmmLS4gOT6WpTNHkRyvkjJpne6J7Vk2cxRx0VFMWbiWg6cueD2SoKCPWIfPXGTKgrVEGSybNYqeSe29HknCRFrXeJbOyuFidR1TFxVwvEIlaF5T0EegU+ermbJwLRWVtSyZkUNGtwSvR5Iwc33PTiyeMZIj5ReZvriAs5U1Xo8U0RT0EeZcVS3TFxdQevqiSsrEr0b068K8R0ew+2gFjz2nEjQvKegjSGVNHbOXFrL98FmeeXg4o1RSJn42YWB3/uOLw1i37xRPPr+BGpWgeUJBHyFq6+r52vKNrCo+yb/fP4S/U0mZBMg9Q3vzz5MG88edZXzrt1uor1cJWqDpvXQRwDnHd17bxrvbj/G9zw7iC8NVUiaBNWV0P06fr+Y/3ttDcnwc3/2sStACSUEfAZ5+ZxcvFh7kq5/qz8yxGV6PIxHqK5/qz+kL1Sz6cC9dEmJ58lMDvB4pYijow9y8FcXkrShhyuh+fP3267weRyKYmfHduwdx5kINP/nDHpLi45gyup/XY0UEBX0YW15wgKff3sXnhvbmB/fcqH8qi+eioox/u38IZy/W8L3Xt5HcIZbPDe3t9VhhTz+MDVPvbDvCt1/byvjrUvjpA0OJUkmZBInY6CieeWQ4I/t14e9f2sSKPSpB8zcFfRj6sOgEX31hE8P6JvOrR4erpEyCTvvYaBZMz6Z/90TmLlvP+v0qQfMnJUCY2XzwDLOXFpLRLYFF00cSH6ezcxKcOrWPZenMHHp0asfMJevYfVQlaP6ioA8jRWUVTF9cQJeOcSydlaOSMgl6KYntWDZrFO1jVYLmTwr6MHHozEWmLCwgOiqKZTNH0aOTSsokNPTtEs/SmaOoqq3n0YVrKauo9HqksKOgDwMnz1UxZcFazlXVsnRmDukqKZMQM7BnIotnjKTsbBXTFq2j/KJK0NqSgj7EVVTWMG1xAYfLL7Jo+kgG9e7k9Ugi12R4WmfypoygqKyCx55bx8VqlaC1FQV9CKurd8xeup5dRyr41SMjGJnexeuRRFol97oUfvbFYRTuP82Tz2/AOfXitAUFfQh7Z9tRVpec5If3Dua267t7PY5Im/jskN58e+IN/GlXGfkfnfB6nLCgoA9RzjnmrSgmo1sCD2b39XockTY1bUw6PTq1I29FsdejhAUFfYhaXXySrYfKeXxcJtG66lXCTFxMFLPGZrCq+CRbSs94PU7IU9CHqHn5JXTr2I4vDO/j9SgifvFQThqJ7WPIW1Hi9SghT0EfgnYcPkv+nuPMuDWd9rHRXo8j4heJ7WN5dHQ/3t52hH0nzns9TkjzKejN7E4z221mRWb2VDOvX29mq82sysy+cdlr+8xsq5ltMrPCtho8ks3PLyYhLppHR6niVcLbjDHpxERFseADHdW3RotBb2bRwDPARGAQ8JCZDbpss1PAV4GfXOHT3OacG+acy27NsAIHT13g91uO8FBOGknxsV6PI+JX3Tu15wvD+/ByYSknzlV5PU7I8uWIPgcocs6VOOeqgeXApKYbOOfKnHPrAF3O5mcLP9iLAbPG6U5REhkez82kuq6e51bt83qUkOVL0PcBDjZ5XNr4nK8c8AczW29ms6+0kZnNNrNCMys8flz91M05fb6aF9cdZNKwPvRK6uD1OCIBkZXSkc8M6sHS1fs5X1Xr9TghyZegb+69e1dzudqtzrnhNJz6ecLMcpvbyDk33zmX7ZzLTklJuYpPHzmWrdnPxZo6Zudmej2KSEDNGZ9F+cUaXlx3sOWN5WN8CfpSoOkVOanAYV+/gHPucON/y4DXaDgVJFfpYnUdS1bt41PXd2dgz0SvxxEJqOFpnclJ78LCD/ZSU1fv9Tghx5egXwcMMLMMM4sDJgNv+PLJzSzBzBIvfQx8Bth2rcNGst+uP8ip89XMHZ/l9Sginpg7IZNDZy7yP1t8Ps6URi3efsg5V2tmTwLvAtHAIufcdjOb2/j6PDPrCRQCnYB6M/s/NLxDpxvwWuNNqWOA551z7/hlJWGstq6eZ1fu5ea0ZEamd/Z6HBFPTLiuO9f16EjeihLuHdZHN7u/Cj7dZ8459xbw1mXPzWvy8VEaTulc7iwwtDUDCryz/SgHTl3g23fdoG9uiVhRUcbs3Cy+8fJmVuw5zoSBKvLzla6MDXKXyssyuyVw+6AeXo8j4ql7hvamV1J75qns7Koo6IPcquKTbDt0ltm5Ki8TuVR2tqbkFJsOnvF6nJChoA9y81YUk5LYjntvVnmZCMDkxrKz+fk6qveVgj6IbT9czsqPTqi8TKSJju1imDK6H29vO6qyMx8p6INY3ooSOraL4RGVl4n8jem3phMbHcX8lSo784WCPkgdPHWBN7ce4eFRaSR1UHmZSFPdE9tz3/BUfru+lOMVKjtriYI+SC38YC9RBjNuTfd6FJGg9Pi4DGpUduYTBX0QOnW+muXrDqi8TOQTZKZ05I5BPVm6ep/KzlqgoA9CS1fvo7KmnjkqLxP5RHPGZ3K2spYXCg54PUpQU9AHmYvVdTy3ah9/d0N3BvRQeZnIJ7k5rTOjMlR21hIFfZB5ef1BTl+oYY7Ky0R8Mnd8FkfKK/n9ZpWdXYmCPog0lJeVMDwtmex+Ki8T8cWEgSkM7JFI3ooSnLuaW2VEDgV9EHlr21EOnrrInPFZKi8T8ZGZMTs3k93HKvjLbt2drjkK+iDhnCNvRTGZKQncfoPKy0Suxj3DetNbZWdXpKAPEh8WnWT74bPMyc0kSuVlIlclNjqKmWMzWLv3FBsPnPZ6nKCjoA8SefkqLxNpjck5aXRqH8P8fNUiXE5BHwS2HWooL5t5awbtYlReJnItOraLYcot/Xhn+1FKjp/zepygoqAPAnn5jeVlo9O8HkUkpE0fk0FsdBTPrtzr9ShBRUHvsYOnLvDmlsM8MiqNTu1VXibSGimJ7bh/RCqvbCilrKLS63GChoLeYwtWlhAdZcy4NcPrUUTCwuPjMlV2dhkFvYdOnqvixcKD3DusDz2T2ns9jkhYyOiWwJ039mTZ6v2cU9kZoKD31NLV+xvKy8arvEykLc0dn8XZylqWq+wMUNB75kJ1LUtX7+PvbuhB/+4qLxNpS0P7JjM6s6HsrLpWZWcKeo+8XFjK6Qs1zNXRvIhfzFHZ2V8p6D1wqbxsRL/OZKd38XockbA04boUru+ZSF5+MfX1kV12pqD3wJtbj1B6+iJzVUUs4jdmxpzxmew5do6/7CnzehxPKegDrKG8rISslAQ+fX13r8cRCWufHXKp7CyyaxEU9AH2QdEJdhw5y5zcLJWXifhZbHQUs8ZlUrD3FBsiuOxMQR9g81YU06NTOybd3NvrUUQiwuSRfUnqEEteBFcYK+gDaGtpOR8WnVR5mUgAJbSLYeot/fjDjmMUR2jZmYI+gPLyi0lsF8NDo1ReJhJI08akExsdxYKVkXmuXkEfIAdOXuCtrUd4eLTKy0QCrVvHdjwwIpVX1h+KyLIzn4LezO40s91mVmRmTzXz+vVmttrMqszsG1ezb6R4dmUJMVFRzFR5mYgnHh+XSW19PYs/3Of1KAHXYtCbWTTwDDARGAQ8ZGaDLtvsFPBV4CfXsG/YO3muipcKD/L5m/vQo5PKy0S8kN4tgYmDe/HrNfupqKzxepyA8uWIPgcocs6VOOeqgeXApKYbOOfKnHPrgMt/91rcNxI8t3o/VbX1PJ6rugMRL83OzaSispblBQe9HiWgfAn6PkDT35XSxud80Zp9w8Kl8rLbB/Wgf/eOXo8jEtGG9k3mlsyuEVd25kvQN3dVj6/FET7va2azzazQzAqPHz/u46cPfi+uO8gZlZeJBI054zM5eraS1zcd8nqUgPEl6EuBvk0epwK+1sH5vK9zbr5zLts5l52SkuLjpw9uNXX1LFi5l5HpnRnRT+VlIsFgfGPZ2fz8kogpO/Ml6NcBA8wsw8zigMnAGz5+/tbsG/Le2nqEQ2cuMidX5WUiwcLMmDs+i4/KzvH+7sgoO2sx6J1ztcCTwLvATuAl59x2M5trZnMBzKynmZUCfw/8k5mVmlmnK+3rr8UEE+cc81aU0L97Rz6l8jKRoHL3kF70Se5AXoSUncX4spFz7i3grcuem9fk46M0nJbxad9IkP/RCXYeOcu/3T9E5WUiQSY2OopZYzP45//Zwfr9pxnRr7PXI/mVroz1k7zG8rJ7h0XUm4xEQsbknL4kx0dG2ZmC3g+2lJ5hVfFJZo3NIC5Gv8UiwSg+Loapo/vx3s5jFJWFd9mZUsgP8vJLGsrLclReJhLMpo5JJy4Cys4U9G1s/8nzvL31CI+M7keiystEglq3ju14IDuVVzccouxs+JadKejb2P+Wl6V7PYqI+OBS2dmiMC47U9C3oRPnqni5sJQvDO9Dd5WXiYSEfl0TmHhTL34TxmVnCvo2tHTVPqrrVF4mEmrm5GZSUVXLCwUHvB7FLxT0beR8VS3Prd7P7Tf0ICtF5WUioWRIajJjshrKzqpq67wep80p6NvIi+sOUn6xhrkTVHcgEormjs/i2NkqXt/ka5VX6FDQt4GaunoWfrCXnPQuDE8L7yvsRMLVuAHduKFXp7AsO1PQt4E3tzSWl6mKWCRkNZSdZVJUdo4/7wqvsjMFfSs1lJcVM6B7R24bqPIykVB2900NZWfzwqwWQUHfSiv2HGfX0QrmjM9SeZlIiIuJjuLxcRkU7j9N4b5TXo/TZhT0rZS3ooSendpzz9DeXo8iIm3gwZGNZWf54VOLoKBvhc0Hz7C6ROVlIuEkPi6Gqbek896OYxSVVXg9TptQOrVCXn4xie1jmJzTt+WNRSRkTLulH+1jo5gfJkf1CvprtO/Eed7edpQpKi8TCTtdO7bjwey+vLbxEMfCoOxMQX+Nnl1ZQmxUFNNVXiYSlh4bm0ldvWPRh3u9HqXVFPTX4HhFFS+vL+W+EX3onqjyMpFwlNY1nrtu6sXzaw5wNsTLzhT01+C5Vfuoqavn8XG6QEoknM0dn0VFVS3Prw3tsjMF/VU6X1XL0tX7uGNQTzJVXiYS1gb3SWJs/24sCvGyMwX9VVq+7iBnK2tVdyASIeaMz6SsoorXN4Zu2ZmC/irU1NWzcGUJORlduFnlZSIRYWz/bgzq1Yl5+cUhW3amoL8Kv998mMPllczV0bxIxDAz5ozPpOT4ef6485jX41wTBb2PnHPkrShhYI9ElZeJRJi7b+pFaucOIVuLoKD30V/2HGf3sQpm52ZipvIykUjSUHaWyfoQLTtT0Pto3l+K6ZXUns+pvEwkIj2QnUrn+NiQrDBW0Ptg44HTrN17SuVlIhHsUtnZH3eW8dGx0Co7U2r5YH5+CZ3axzA5J83rUUTEQ9PGpIdk2ZmCvgV7T5znne1HmXJLPzq2i/F6HBHxUJeEOL6Y3ZffbTrE0fLQKTtT0Ldgfn4JsdFRTBuT7vUoIhIEHhsXemVnCvpPUFZRySsbSrlveKrKy0QEgL5d4rl7SG+eX3uA8ouhUXbmU9Cb2Z1mttvMiszsqWZeNzP7RePrW8xseJPX9pnZVjPbZGaFbTm8v10qL5udqwukROR/zcnN5FwIlZ21GPRmFg08A0wEBgEPmdmgyzabCAxo/DUb+NVlr9/mnBvmnMtu/ciBca6qlmWr93PnjT3J6Jbg9TgiEkQG90li3IBuLPowNMrOfDmizwGKnHMlzrlqYDkw6bJtJgFLXYM1QLKZ9WrjWQNqecEBzlbW6mheRJo1JzeL4xVVvLbhkNejtMiXoO8DHGzyuLTxOV+3ccAfzGy9mc2+1kEDqbq2noUf7GWUystE5Apu7d+VG3t3Yn5+SdCXnfkS9M1d73/5qj5pm1udc8NpOL3zhJnlNvtFzGabWaGZFR4/ftyHsfzn95sPc6S8krkTsjydQ0SCl5kxd3wWJSfO816Ql535EvSlQN8mj1OBy4uZr7iNc+7Sf8uA12g4FfQxzrn5zrls51x2SkqKb9P7gXOOvPxiBvZIZMJ13s0hIsFv4uCe9O3SgXkrinEueI/qfQn6dcAAM8swszhgMvDGZdu8AUxtfPfNaKDcOXfEzBLMLBHAzBKAzwDb2nD+Nvf+7jL2HDvHnPEqLxORT3ap7GzjgTOs23fa63GuqMWgd87VAk8C7wI7gZecc9vNbK6ZzW3c7C2gBCgCngW+3Ph8D+ADM9sMFABvOufeaeM1tKl5K0rorfIyEfHRAyP60iUhjrwgLjvz6Zp+59xbNIR50+fmNfnYAU80s18JMLSVMwbMhgOnKdh7iu9+dhCx0bqWTERa1iEummm3pPOzP+5hz7EKruuR6PVIH6M0a2L+ihKSOsQyeWTfljcWEWk09ZZ+dIiNDtqyMwV9o+Lj53h3x1GmjO5HgsrLROQqdE6I44sj+/L6pkMcKb/o9Tgfo6BvtGClystE5NrNGptBvYNFHwRf2ZmCnsbysvWHeGBEKimJ7bweR0RCUN8u8Xx2SK+gLDtT0ANLPtxHTX09j49T3YGIXLvZuZmcr67jN2v3ez3K34j4oK+orGHZmv1MHNyTdJWXiUgr3Ni7sezsg31U1gRP2VnEB/3ygoNUVNYyJ1d1ByLSenPHZ3HiXBWvbQyesrOIDvpL5WW3ZHZlaN9kr8cRkTAwJqsrN/VJ4tn8EuqCpOwsooP+jc2HOXq2kjnjdW5eRNqGmTFnfGZD2dmO4Cg7i9igr6935K0o5vqeiYxXeZmItKE7b+xJWpf4oCk7i9igf393GR+VqbxMRNpeQ9lZBpsOnqFg7ymvx4ncoM9bUUKf5A58dojKy0Sk7T2Q3ZeuCXHkBUEtQkQG/fr9pynYd4pZYzNUXiYiftE+NpppY9L5864ydh+t8HSWiEy5vBXFJHWI5YsqLxMRP5oyuqHsLC/f2wrjiAv6orJzvLfzGFNvUXmZiPjXpbKzNzYd5vAZ78rOIi7oF6wsIU7lZSISII+Ny8DhbdlZRAV92dlKXt1wiAeyU+nWUeVlIuJ/qZ3j+dyQXrxQcIDyC96UnUVU0C/6cB+19fU8NlYXSIlI4MzOzeJ8dR2/9qjsLGKCvqKyht+s2c/Ewb1UXiYiATWodydyr0th8Yd7PSk7i5igf6HgABVVtao7EBFPzB2fyYlz1by6IfBlZxER9JfKy8ZkdWVIarLX44hIBLolsytDUpN4dmXgy84iIuh/t+kQx85WMWe8qohFxBtmxpzcLPaeOM8fth8N6NcO+6Cvr3fMzy/hhl6dyB3QzetxRCSC3Tm4J/26Br7sLOyD/s+7yigqO8dclZeJiMeio4zHx2WyubSctQEsOwv7oM/LL6ZPcgfuuqmX16OIiHD/iNSGsrMVgatFCOugX7//FOv2neaxcSovE5Hg0D42mulj0nl/93F2HT0bkK8Z1uk3b0UJyfEqLxOR4DLlln7Ex0Uzf0VgKozDNuiLys7x3o5jTL0lnfg4lZeJSPBIjo9j8sg03th8mEMBKDsL26B/Nr+EdjFRTLuln9ejiIh8zKwAlp2FZdAfO1vJaxsP8WB2X7qqvExEglCf5A7cM7Q3LxQc4MyFar9+rbAM+kUf7qW2vp7Hx6nuQESC1+zcTC5U1/HrNf4tOwu7oD9bWcPzaw5w1029SOsa7/U4IiJXdEOvTkwYmMKSVfv8WnbmU9Cb2Z1mttvMiszsqWZeNzP7RePrW8xsuK/7trUX1jaWl+Wq7kBEgt+c3CxOnKvmlQ2lfvsaLQa9mUUDzwATgUHAQ2Y26LLNJgIDGn/NBn51Ffu2maraOhZ+sJdb+3flptQkf30ZEZE2MzqzC0NTk3g2339lZ74c0ecARc65EudcNbAcmHTZNpOApa7BGiDZzHr5uG+beX3jYcoqqnQ0LyIhw8yYMz6LfScv8K6fys58Cfo+wMEmj0sbn/NlG1/2bRP19Y68/GIG9erEOJWXiUgIuePGnqR3jSfPT2VnvlxJ1FwT2OWTXGkbX/Zt+ARms2k47UNaWpoPY/2tCzV1jEzvwrgBKSovE5GQEh1lPHFbf7aUllNVW0/72Og2/fy+BH0p0LRDIBU47OM2cT7sC4Bzbj4wHyA7O/uq/5fWsV0MT9835Gp3ExEJCg9k9+WBbP/Utfhy6mYdMMDMMswsDpgMvHHZNm8AUxvffTMaKHfOHfFxXxER8aMWj+idc7Vm9iTwLhANLHLObTezuY2vzwPeAu4CioALwIxP2tcvKxERkWZZIO9y4qvs7GxXWFjo9RgiIiHDzNY757Kbey3srowVEZG/paAXEQlzCnoRkTCnoBcRCXMKehGRMBeU77oxs+PAtRY0dwNOtOE4oUBrDn+Rtl7Qmq9WP+dcSnMvBGXQt4aZFV7pLUbhSmsOf5G2XtCa25JO3YiIhDkFvYhImAvHoJ/v9QAe0JrDX6StF7TmNhN25+hFRORvheMRvYiINBGSQd+am5WHKh/W/EjjWreY2SozG+rFnG3J1xvLm9lIM6szs/sDOZ8/+LJmM5tgZpvMbLuZrQj0jG3Nh+/tJDP7vZltblzzDC/mbCtmtsjMysxs2xVeb/v8cs6F1C8a6o6LgUwabmyyGRh02TZ3AW/TcIer0cBar+cOwJrHAJ0bP54YCWtust2faajKvt/ruQPw55wM7ADSGh9393ruAKz528C/Nn6cApwC4ryevRVrzgWGA9uu8Hqb51coHtG35mbloarFNTvnVjnnTjc+XEPD3bxCma83lv8K8ApQFsjh/MSXNT8MvOqcOwDgnAv1dfuyZgckWsM9QjvSEPS1gR2z7Tjn8mlYw5W0eX6FYtC35mbloepq1zOLhiOCUNbims2sD/B5YF4A5/InX/6crwM6m9lfzGy9mU0N2HT+4cua/xu4gYbbkG4Fvuacqw/MeJ5o8/zy5Z6xwaY1NysPVVdzk/XbaAj6sX6dyP98WfPPgW855+rC5Ibwvqw5BhgBfBroAKw2szXOuT3+Hs5PfFnzHcAm4FNAFvCema10zp3182xeafP8CsWgb83NykOVT+sxsyHAAmCic+5kgGbzF1/WnA0sbwz5bsBdZlbrnPtdQCZse75+b59wzp0HzptZPjAUCNWg92XNM4CnXcMJ7CIz2wtcDxQEZsSAa/P8CsVTN625WXmoanHNZpYGvApMCeGju6ZaXLNzLsM5l+6cSwd+C3w5hEMefPvefh0YZ2YxZhYPjAJ2BnjOtuTLmg/Q8C8YzKwHMBAoCeiUgdXm+RVyR/SuFTcrD1U+rvl7QFfgl41HuLUuhAuhfFxzWPFlzc65nWb2DrAFqAcWOOeafZteKPDxz/mHwBIz20rDaY1vOedCttXSzF4AJgDdzKwU+D4QC/7LL10ZKyIS5kLx1I2IiFwFBb2ISJhT0IuIhDkFvYhImFPQi4iEOQW9iEiYU9CLiIQ5Bb2ISJj7/zuLmZsyv11JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fn.plot(un)"
   ]
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}