{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute numerically with numpy and scipy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Numerical Solving of Systems of Linear Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have been using numpy already in our plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_168/1288942697.py:2: UserWarning: get_systems() requires Cython profiling to be enabled, otherwise the results will be very unreliable. Rebuild Sage with the environment variable 'SAGE_PROFILE=yes' to enable profiling.\n",
      "  get_systems(\"point([0,0], size=50).show(figsize=1, axes=False)\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFkAAAA4CAYAAACWo1RQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAABg0lEQVR4nO3aMWoCQRjF8bdiYWe3CFtvqaW9nY21XsWj6BG0tBEUPIPYegE9gYq4KYZAIBvGje77ZpL3AyE6Uwx/RpmMJkVRQOrVsF7Af6DIBIpMoMgEikygyASKTKDIBM0Kc/VfS7nEN0E7mUCRCRSZQJEJFJlAkQkUmUCRCRSZQJEJFJlAkQmqXBCZu1yAw8H93e0CrZbtep4VxU6+XoHpFMgyoN93jyxzr91u1qvzSyr87sLkqvN+B0YjYL0uHx8OgdUKaNq9J+O/6lwsfg4MuLHlkree3wg+8nzunzOb1b+OVwQf+Xh8zxxLwUdut98zx1Lwkcdj/5zJpP51vCL408X5DPR6wOlUPt7pAPs9kKbcdX0R/+kiTYHtFsjz72N5Dmw2poGfEvxO/vR4uOPabueeDwbujNyw3ybenRxN5IDF/3HxFygygSITKDKBIhMoMoEiEygygSITKDKBIhMoMkGV73i9FyFSTjuZQJEJFJlAkQkUmUCRCRSZQJEJFJngA6rZQkb8UbldAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "['numpy']"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sage.misc.citation import get_systems\n",
    "get_systems(\"point([0,0], size=50).show(figsize=1, axes=False)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize large data sets, we can use matrix_plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAFDCAYAAADVkhLhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgp0lEQVR4nO3df5hkVXng8e/LMAM8MN0qBkcdSDCwoi4GHtB13RXUSPDHxpBsViKrgBhJcIwaNokLZgU0aqJkYpRRokLk0VVJjLqSZFzUABpFDKxR1GhQCD+EAZHH6THB6R7m7B9VXdPT86P7Pbeq61b19/M89fR01X3rnL7Tfd865977niilIEkSwD7D7oAkqT1MCpKkHpOCJKnHpCBJ6jEpSJJ6TAqSpB6TgiSpx6QgSeoxKUiSekwKkqSekU8KEfHKiLgtIn4SETdFxDOG3aelFBEXRkSZ99g07H4NWkScEBFXRcTd3Z/5lHmvR3ff3B0RD0bEtRHxpCF1d2AWsR8+sJvfjy8PqbsaASOdFCLiVOAdwJuBY4EvABsj4rBh9msIvgk8es7j6OF2Z0kcCHwNeNUeXv894Nzu608BNgGfiYjVS9O9JbPQfgD4NDv/fjx/CfqlEbXvsDvQ0LnAZaWU93e/f21EnAycA5w3vG4tuW2llLEfHcxVStkIbASIiJ1ei84TrwXeXEr5ePe5M4B7gdOAP1vKvg7S3vbDHFuX2++H6o3sSCEiVgHHAVfPe+lq4OlL36OhOrI7fXBbRHw0Ih437A4N2eHAGub8bpRStgLXsfx+NwCeGRH3RcQ/R8T7IuKQYXdI7TWySQF4JLCCzqe/ue6lc0BYLm4ATgdOBl5B52f/UkQcPNReDdfs//9y/92AzijivwPPBv4Hnam0v4uI/YbaK7XWqE8fAcxfECJ289zY6k4fzLo5Iq4HvgecAawfTq9aY1n/bgCUUq6c8+03IuJG4HbgBcDHh9MrtdkojxTuBx5i109+h7DrJ8Rlo5Tyr8DNwJHD7ssQzc6f+7sxTynlHjpJYTn/fmgvRjYplFKmgZuAk+a9dBLwpaXvUTt0pwWeANwz7L4M0W10EkPvd6N7DupElvHvBkB3WvFQlvfvh/Zi1KeP1gMf7A6JrwfOBg4DLh1qr5ZQRFwMXAXcQeeT8O8DE8AVw+zXoEXEQcARc546PCKOAR4opdwREe8Azo+IW4BbgPOBfwM+vNR9HaS97Yfu40Lgr+gkgZ8B3kJnlP2JpeynRsdIJ4VSypXdTz5voHP99TeA55dSbh9uz5bUWuAjdE68/wD4MvC0ZbAPjgeumfP97PmTK4AzgbcBBwDvBh5O54T8L5RStixhH5fC3vbDOXTuWTkdeBidxHANcOoY7gf1SZSyrM67SZL2YmTPKUiS+s+kIEnqMSlIknpMCpKkHpOCJKnHpCBJ6hn5pBAR+3UXU1nWBb7cDzu4LzrcD6ox8vcpRMQEsBmYLKVMDbs/w+J+2MF90eF+UI2RHyk0ERHrximmVpt/JvdDfUytNv9MS7kf2iAizomIr0fEVPdxfUQ8b6BtLueRQkR8q5TyxHGIcT/sFFO1L1r+M7kfKmNGWUT8Ip1q0N/tPnUG8LvAsaWUbw6izdbVPuoupfgYYLG1WWbX3F29l+UI92Sf7h/OOMS4H3ao3Rdt/pncD3Uxq4G7y4A+/UbE/sCqitDpUspPFtqolHLVvKdeHxHnAE+jszZ737VupBARjwXuGnY/JI2NtaWU7/f7TSNi/zXs8+AmtteEbwKeDGyd89zW7rKxe2pvBfDf6BQ7PLaU8q2ahhfSupEC3RHCnXfeycRE9kOEJHVMTU1x6KGHwuJnHbJWbWI7d/IYJhKnZ6fYzqHcvQa4b95LF9Epdb6TiDiaztIA+wM/Bn55UAkB2pkUAJiYmDApSGq9CfZPJQV2jCzWsnPC2tMo4TvAMXTKn/9X4IqIOHE5jRQkaYSsBFYktn9o9h9bFnMBQHeVydkTzTdGxFOA1wC/kenlYi3rS1IlaQQFMLAbEkdrpPDp9FU1xP0VJ9K/mA8B4D3zl4te2Kb4TDpmTflBOuap8VPpGKi7rOLvy+XpmL+Ls9Ixzy4vTccAvIEPpWMuekS+nXjgzfmgQ87PxwDlNyr+Nt5U8zvxhIqYMypi4Od5eTrmswzjwplVVI4UFhQRbwE2AnfSuZLq14BnAs9NNJgyWklBklonO32UmqB5FPBBOssNbwa+Djy3lJL/NLlIJgVJamQluUPptkVvWUrJD5caGtg5hYh4ZUTcFhE/iYibIuIZg2pLkoZnVcWjvQaSFCLiVOAdwJuBY4EvABsj4rBBtCdJw2NSWIxzgctKKe8vpfxTKeW1dE6UnDOg9iRpSFaSSwgrh9PNRer7OYWIWAUcB/zhvJeuBp6+m+33Y+fLq1bP30aS2it7oM9fKbaUBnGi+ZF0TsXfO+/5e4E1u9n+POCCAfRDkpbAeCWFQd68Nv+C4djNcwBvBSbnPNYOsE+S1GdOHy3kfjp3Z8wfFRzCrqMHulUBezU/Kso+S9IQtf/kcUbfRwrdOh03AfNv7z0J+FK/25Ok4Rqvq48GdfPaeuCDEXEjnZKvZwOHAZcOqD1JGpLZ6aPFatcaNvMNJCmUUq6MiIOBN9C5PfsbwPNLKbcPoj1JGp6V5M4TVC3Ks2QGVuailPJu4N39fM/9K0pAlYpzFFHyxdkAruCz6Zg119W0dEc64ob7D6lpiNi+NJ9qnlMOqoiq2nl8saIgXjzw6nTMabw+HfO//y0fAxBvek9F1D+kIz5CvtjhQXw+HQPw3aj43Wv3h/CRYO0jSWoke56g3ZnLpCBJjWTPKSzT6SNJWh6yIwWTgiSNMZOCJKknO320+JXXhsGkIEmNZEcKJoVFiYh1wDoGW49JkvrMpDAQpZQNwIaImKCzFqkkjYDs9NHil+MchtYkBUkaTdmRgklBksZYtszF8iudLUnLSHb6aGZQHekLk4IkNZKdPjIp9E1NFfKra8qMPHhZRRC84YB8sbAjT8y3s6Ycn475dwfn2+n4ZD7kQ/mCgjMv/dV0zIpyTDoG4B+rCq3NXx5kYesr2nleOmLWxRUxX0tHvPio/M/06G/XLZx1z7fzMa+taqkpk4IkqSc7fTQ9qI70hfcESJJ6HClIUiPZ6aPluRynJC0TJgVJUk/2nIL3KUjSGHOkIEnqMSkMhFVSJY0mp48GwiqpkkZRYV9K4lCa2XYY2t07SWq57ezL9sShNLPtMLS7d5LUciYFSVKPSUGS1LONfdiWuD4ms+0wjFRS2MJz0jEnv/Mz+YZefUg+Brj9s/kKkj9f8hUkP8e/pmP+iAPTMQCv++Qp+aAn5ENWlC3pmNPi1nxDwEfLQfmgo36cDlnzonwz+165fz4IqKmv+r64MR3zinJdOuYe6v6fPnJURVBNVeSGpsmVuMtsGxHnAb8CHAU8CHwJeF0p5TuJt0lpd8qSpJabYUdiWMwjWTj7RGAD8DTgJDof5K+OiLpPeYswUiMFSVpOSinPnft9RLwMuA84Dvj8INo0KUhSA4OcPtqNye7XB5q9zZ6ZFCSpgQZJYXXETucUt5ZStu4pLjobrwf+vpTyjVwvF89zCpLUQINzCnfRqd4w+zhvgaYuAZ4MvLiP3d+FIwVJaqDBSGEtMPeyu72NEt4FvBA4oZRyV66HOa1JChbEkzSKZshdUTRn2y2llKm9bdudMnoX8MvAM0spt1V0MaU1ScGCeJJG0ez0UWb7hA3AacAvAVsiYk33+c2llAdzb7U4rUkKkjSKpskVw05efXRO9+u1855/GfCB3FstjklBkhoYZFIopaLkQUMmBUlqYMDTR0vOpCBJDUyTO5A2vHlt4EYsKTwrH/LqU9Mh3+UH+XaAI55zaTrmc5+vqOB1wifTIb/Pinw7AKfk2+I7p1TEPDYd8uFyUr4dgHfmi9v9SUX5sT/4dn7kf3rUVXRb/8/5tqLqwsZjKmKuqWmIF5dz0zG/xieq2mrCpCBJ6nH6SJLUMw2pcbgjBUkaY+OWFLx7WJLU40hBkhrYRu48wbZBdaRPTAqS1MA0uSmXtk8ftSYpWBBP0igyKQyIBfEkjaIZcknBS1IlaYxNA5lbBx0pSNIYMylIknpmyCUFp48kaYxlP/k7Uuij03h9OiZf+gyOWFtXlKxmOe33vT1fyOwVvCcdM/O7lVdHv/0z6ZAnPT7/M33z9RX7/Iv5vgGsvjbfv3Mr+le4Lx1zRGX5/BUcWBFVs/9elo44mO9XtAPvrYpaeiYFSVJPdjrI6SNJGmPTQGYc2fak0PcbxSLiwogo8x6b+t2OJLXBTMWjzQY1Uvgm8Jw53z80oHYkaajGbaQwqKSwrZTi6ECSRsygksKREXE3sBW4ATi/lHLr7jaMiP2A/eY8tXpAfZKkvpshN1Joe5XUQRSfuwE4HTgZeAWwBvhSRBy8h+3Po1PraPZRtXKsJA3DdMWjzfo+UiilbJzz7c0RcT3wPeAMYP1uQt467/nVmBgkjYhpYHti+7aPFAZ+SWop5V8j4mbgyD28vpXONBMAEXU370jSMMyQSwptv+pm4Emhe87gCcAXBt2WJC217BrNyy4pRMTFwFXAHcAhwO8DE8AV/W5LkobNpLCwtcBHgEcCPwC+DDytlHL7ANqSpKGaIXegz0w1DUOUUln8bUBmV17bvHkzExMTO7/Gzen3OzuOTse8dzId0vGjUyuCnpCOWMkfpGOeVvn55IsVhda231xRhvDok9Ih09OfzbcDPG9VPuaByP+d/L9fz58f2/dDdX+PL38w39Z7H6xo65fyIXy3IgbY59aD0jEPzSmBOTU1xeTkJMBkKWWqrhd7NnusehSb2YeJBbeftZ0p7mVw/WrK2keS1EB2jea2jxRMCpLUwDZyi+y0a25mVyYFSWoguxynSUGSxtj2FZC5vaoUWn0JkklBkhrYvq9JQZLUVZUUti642dAMoiBelYhYFxHfAr4y7L5I0nLVmpFCKWUDsGH22t9h90eSFqPsC2WMzjS3JilI0khayVjdqNCa6SNJGkmrKh4JEXFCRFwVEXd317w/pW993w2TgiQ1MeCkABwIfA14VV/6uwCnjySpiZUMtExqd+GyjbA0682YFCSpiWxS2DE/s3reQX5rd9GxoRqtpHBxRcXT8vR0zKu4Ph0DcMmpFZcVnJwP+b8vvzAd8+zyuHxDAG+9NR9z3k3pkHJ7vuJp/PSL0jEAL+Qv0jHXlf3TMbE5//vwqfenQwB4ITemY848IP+p86qr0yH8kJfkg4DtF1VU272gqqlmVlE7Upi/7PBFwIXNO9TMaCUFSWqbleSOpDsWaV4LbJnzytBHCWBSkKRmVpE7ku6YPtriegqSNG7qk0IrmRQkqYmV3cdiJU/lRMRBwBFznjo8Io4BHiil3JF7t4WZFCSpiVUMNCkAxwPXzPl+fffrFcCZ6XdbgElBkpoYcFIopVybj6rXmqQQEeuAdbR+xk2S5lhJzV3KrdWaA3ApZUMp5YnAU4fdF0larlozUpCkkVRXz6i1TAqS1IRJQZLUsy+5E80tX0/BpCBJTWRHCq681kfraoLelY64ZKdLghPWVMSc9fJ0yPfOuryiocpKYc/Kh5wdx6Vj4i35v5QXnp8OAeBTJ16ZjinX5f9zv5hao7HjP5eKAoTAz3N8OuZrkd/nPzw8HcJ/uLXuasobLpiuiMr/3zZmUpAk9WQvSXX6SJLGWHakYFKQpDFmUpAk9WSnj5LLcS41k4IkNZEdKZgUJGmMmRQGw4J4kkZSdvpo28KbDFNrkkIpZQOwISImgM3D7o8kLUp2pNDypOCncklST2tGCpI0krLLcWa2HQKTgiQ1MWbTRyYFSWoie6J5ZlAd6Y/RSgoH3FwRdFA64hn8bkU78IVf/52KqJPTEa/4/mXpmN9emw4B4E/KTemYa0q+OBs356uEfeq81+TbAZ7LO9MxwQfTMV/hpekYPltRcQ743Gn5/Vd4bzombj0rHXPDxZUV4P6hImYI9fDSIwWTgiSNMZOCJKknO31UUxF8CZkUJKmJ7Eih5Ut3mhQkqQmTgiSpJzt95H0KkjTGHCkMhgXxJI0k72geDAviSRpJTh9JknqcPpIk9YxZUnD+XpLU40hBkppY8VDnkdm+xUYsKdTcH/6IdMSfR10Br2MrYrbc96J0zE/WRDpm/3JXOgbgSPLF7X5csf/KH+d/pvituv+nT1/7xnzQ5sl0yFMn35aO+V8npUMAeNML8jHx17fkg27OHzLOrKsvyQfKbRVRdQUFm5kmd2xqd52LEUsKktQ2JgVJUs8MuQN9u8ukpk80R8QJEXFVRNwdESUiTpn3ekTEhd3XH4yIayPiSX3rsSS1ynTFo71qrj46EPga8Ko9vP57wLnd158CbAI+ExGrq3ooSa02XkkhPX1UStkIbASI2PnkYHSeeC3w5lLKx7vPnQHcC5wG/Fmz7kpS2yzz6aMFHA6sAa6efaKUshW4Dnj67gIiYr+ImJh9AI4oJI2QmYpHTkS8MiJui4ifRMRNEfGMPnV+F/1OCmu6X++d9/y9c16b7zw6tY5mH3XXTkrSUAx2+igiTgXeAbyZzpXvXwA2RsRhfen+PIO6o3n+BeSxm+dmvRWYnPOoXGJekoZhdvposY/0SOFc4LJSyvtLKf9USnktcCdwTj96P1+/L0nd1P26BrhnzvOHsOvoAehNL22d/X7+eQpJardpcqVPeyOF1fOOd1u7x8OeiFgFHAf84bw3uZo9TMk31e+Rwm10EkPvvszuD3Ui8KU+tyVJLVA9fXQXO0+dn7ebN38ksILclHwj6ZFCRBwEHDHnqcMj4hjggVLKHRHxDuD8iLgFuAU4H/g34MPNuytJbVN99dFaYMucF7buum1PZkq+kZrpo+OBa+Z8v7779QrgTOBtwAHAu4GHAzcAv1BKmfvDS9Jyt6WUMrXANvcDD7HrqGCPU/JN1dyncC2dLLWn1wtwYfchSWNumtyhdPGjilLKdETcRGdK/hNzXjoJ+D+JRhdtxGof1dwJeM3Cm8xzRPmfFe3AsXw/HfNVXpKOeUlNFddSV6rylovybf12RTuPPLci6NyaKppwFY9Lx/zcw/L74bDyrHTMO9IRHU/6m3zMNz/09nTMn7wkfyHIb5dT0jEAj45PLLzRfAOZUFnI4JJC13rggxFxI3A9cDZwGHBp9o0WY8SSgiS1zWDvaC6lXBkRBwNvAB4NfAN4finl9tQbLZJJQZIamaZzgVBm+5xSyrvpnKcdOJOCJDUy+KSwlEwKktTIDLlDabsL4pkUJKmRGXIjBZOCJI2xaXLFIZw+WpSIWAesY3BF+iRpAEwKA1FK2QBs6K6psHnY/ZGkxZkhlxScPpKkMTbNXoo87GH79jIpSFIjJgVJUs8MuaTg9JEkjbHsJ39HCn30vYqYiuJ2N9UtE/3VgyqCHv+ddMjHbsw3849cnA8Cpi/Ixz21pirZ31bEPL9u5dZ/rIj5xVKz8uEr0xFbyokV7cBf8vl0zBM/X1Ps8Oh0DKdWFLYD7inHV0RV/HFoJyOWFCSpbRwpSJJ6tpE7p7BtUB3pC5OCJDUyTW4hB080S9IYMylIknpmyCUFp48kaYxNA9sT25sUFsWCeJJGk0lhICyIJ2k0zZBLCg8NqiN90ZqkIEmjKbscp0lBksaYSUGS1DND7kCfmWpaeiYFSWoku/KaSaFvXhcvScf8UakoiHfc0/MxAPxKRcyqfMiax6dDjrmoouAc8NQL8zGrK9p5RMmUCei4Iep+pjXlB/mgi38qH7MuH8KrrqsIgidedmk+6JqKhk54Vj7myssrGoI/4qaquKWXXXmt3UnByz8lST0jNVKQpPbJLrJTN8JdKiYFSWog2EYkkkKhtDotmBQkqYF9KpJCmy9KNSlIUgMmBUlSzz48NEZnFFqUFCyIJ2kUrSJ/mnnrgPrSD61JChbEkzSKVjJOdym0KClI0ihahUlBktRlUpAk9ezLONVINSlIUiOrGK+kEKW06wKp2RPNmzdvZmJiYufXuL7iHd+WjtgQn6hoBx6oKOp2b0U7l/A36ZifjudXtAS3/2k+5thX5/fDoflm+BQVReoA+E/5kEMqfifuyYfwvIoYgIo6dY8+Px/z4orf8fX/pe4Y86q/zrf1rjkXfE5NTTE5OQkwWUqZqurEXsweq568GVZMLLh5z0NT8PVJBtavphwpSFIDK8kdSAe5QnNEvB54AXAMMF1KeVj2PUwKktTAKnIH0gHfiLUK+EvgeuDlNW9gUpCkMVFKuQAgIs6sfQ+TgiQ10GCksDpip/MmW0spQ7/Z2ZISktTASjqJYbGPlTtC76JTvWH2cd7S9XrPHClIUgPzDvQLmjM2WAtsmfPSbkcJEXEhcMECb/uUUsqNiW7sUWuSggXxJI2iBklhyyIvSb0E+OgC2/xLogt71ZqkYEE8SaOoQVJYlFLK/cD9ybBqrUkKkjSKZs8ptEFEHAY8AjgMWBERx3Rf+m4p5ceLeQ+TgiQ1sJLcSGHANSTeCJwx5/uvdr8+C7h2MW9gUpCkBmavKmqDUsqZwJlN3sOkIEkNZKeP2lVtblcmBUlqIDtSMCn000X/MR2y+oJPpmPWpCM6zq2I2cpL0jH78KJ0zO03Luoc0y5WH5evVPlV3pOPufk30zEcXVXaBV5zWT7mvooyZkdV/Hl9IB8CwNNemg65p+JWqbdU1H1e/6h8DMAlzKRj3lXXVCMmBUlST3b6qO0rr6VvFIuIEyLiqoi4OyJKRJwy7/UPdJ+f+/hy33osSS2SKXHRppPSe1IzUjgQ+Brw58Bf7WGbTwMvm/P9dEU7ktR62QN920cK6aRQStkIbASYV+Fvrq2llE0N+iVJGoJBnVN4ZkTcB/wIuA54fSnlvt1tGBH7AfvNeWr1gPokSX2XPafQ9jWaB5EUNtJZ+ed24HDgTcDfRcRxe6gVfh4LVwCUpFbKTh8tu6RQSrlyzrffiIgb6SSIFwAf303IW4H1c75fTafOuCS1nkkhqZRyT0TcDhy5h9e3MqeO+F7OU0hS6+xLrvZR/u6LpTXwpBARBwOHAvcMui1JWmrZkULFbZBLKp0UIuIg4Ig5Tx3eLc/6QPdxIZ1LVe8BfgZ4C51a4J9o1lVJap9lnxSA44Fr5nw/ez7gCuAc4GjgdOBhdBLDNcCppZS5y85J0ljIXn00dtNHpZRr2fviQSdX90aSRkx2pDB2SWGoLnhLOuRjFc2cnK+713FIRamrz+VDHn50/mT8D487Pt8QsIV81bR9+J10zM8efU465panV5YWe2NFzOUVfyrfvrmioSdUxFR68gfTIZNf/1C+ncvy/7cA+3FpRdTSl5szKUiSerLTR22v+WNSkKQGsiOFcSyIJ0nqMilIknqy00eZG92GIb2egiRpfDlSkKQGVpL79N/2kUJrkkJErAPW4ehF0gjxnMKAlFI2ABsiYgLYPOz+SNJijNs5hdYkBUkaRY4UJEk9+27vPDLbt5lJQZIa2Gdb55HZvs1MCpLUwLglhShl6QtI7c3siebNmzczMTGx82t8J/+Glz4+H/ObdQW8ePA9+ZgDKtq6qaKd496ZjwE28Zp0zBoem47Zhx+lY7bz6nQMwENTb03HHDWZ/zv5cckXLrzn8rq/x/1eno85rKKdW8rl6ZjVcVZFS7DlyoW3ma+8aMe/p6ammJycBJgspUxVdWIvZo9VP7oT5h2q9mpqCh52KAPrV1OOFCSpgdgGkSh9Gi0fKZgUJKmJaXKlT1teJtWkIElNmBQkST0z5FbOafkqOyYFSWpihtynf5OCJI2xMZs+ak3xuYhYFxHfAr4y7L5I0nLVmpGCBfEkjaQxGym0JilI0kgas3MKrZk+kqSRNF3xGICI+JmIuCwibouIByPiexFxUUSkCrM6UpCkJqbJLZIwuOmjo+h80P8N4LvAvwfeBxwI/M5i38SkIElNtGT6qJTyaeDTc566NSIeD5zDuCaFX+WodMzHfrOmGFdlKj/gxHTIFXFdOmZTOgJeV75YEQWP+nBF0A13pUO2/+n5FQ3VlHSDcxLFy2bdkv+RoKa43VnHVzQEW8/6vXTMLfxDvqGb354OWVdRGBBgZVTsvxctvEnfTZM7ku44vKyO2GnfbC2lbO1Tr2ZNAg9kAjynIElN1J9TuIvOlZazj/P62a2I+Fngt4BLM3EjNVKQpNapnz5aC2yZ88puRwkRcSFwwQLv+pRSyo1zYh5DZyrpL0sp70/0zqQgSY3MkDuS7kgKWxa5nsIlwEcX2OZfZv/RTQjXANcDZyd6BpgUJKmZaWBFcvuEUsr9wP2L2TYiHksnIdwEvKyUkl4R2qQgSU3MkEsKA7r6qDtCuBa4g87VRj81eyK7lLLo61NMCpLUxDS5S3YGd5/CLwBHdB/zr5db9CVgrbn6yIJ4kkZSS+5oLqV8oJQSu3tk3qc1SaGUsqGU8kTgqcPuiyQtV04fSVITM+Q+Xre8IJ5JQZKamCYxY4+lsyVprJkUJEk9M+SSgtNHkjTGsp/8HSn0z8eYTMes5vJ0TO3lT58jX5H19IorcKP8bTrmbP4iHQMQJ+UrVZa/qqiK+cN8yPMOzscAUFF981crKn1+7M6KKp+X37jwNrvxjLPy/fsCR+cbOvrn0iGn11Q7BZ5Y1qZj3rjL5flLwJGCJKkne5A3KUjSGJsGMoMhk4IkjTGTgiSpZ4ZcUtg2qI70h0lBkpqYBjIFqk0KixMR64B1tKgekyQtyKQwGKWUDcCGiJigs16pJLXfDLmk8NCgOtIfrUkKkjSSsiuvtTwpOFUjSepxpCBJTYzZSMGkIElNzJA70GfOPwyBSUGSmsgustPypBCl1BWrGpTZq482b97MxMTEsLsjaURNTU0xOTkJMFlKmer3+/eOVYfARCIpTG2HyfsYWL+acqQgSU1kq6S263P4LkwKktREduU1k4IkjTGTwtKYmmrdVJukEbJUx5CpZNXTth/Z2pgUVgMceuihw+6HpPGwmsEci6eBTYfCmorYTbR0Yc42Xn0UwGOALYsMWQ3cBaxNxMz6CvnVN9sa437YoXZftPlncj/UxawG7i4DOtBFxP7AqorQ6VLKT/rdn35o3Uih+5/3/cVu38khAGzJXt4VEdvHJcb9sFPM7D9T+6LlP5P7oS5moLM13QN7Kw/utZZ77aMNYxZTq80/k/uhPqZWm3+mpdwPy1Lrpo+y5pTabuWNIEvF/bCD+6LD/aAa4zBS2Apc1P26nLkfdnBfdLgflDbyIwVJUv+Mw0hBktQnJgVJUo9JQZLUY1KQJPWYFCRJPSYFSVKPSUGS1GNSkCT1/H8UeGlzrsw+iAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.random.normal(0, 1, (20, 20))\n",
    "matrix_plot(A, cmap='hsv', colorbar=True, figsize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already a random coefficient matrix A.\n",
    "Let us generate a random right hand side vector b and solve a linear system A*x = b.\n",
    "The random numbers are uniformly distributed in the interval [-1,+1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00512064],\n",
       "       [-0.70480263],\n",
       "       [ 1.44623146],\n",
       "       [-0.99364059],\n",
       "       [-0.74074156],\n",
       "       [ 2.72227682],\n",
       "       [-1.16717296],\n",
       "       [-0.1015867 ],\n",
       "       [-1.39478794],\n",
       "       [ 1.54941789],\n",
       "       [ 2.09291182],\n",
       "       [ 1.23166181],\n",
       "       [-0.31414951],\n",
       "       [-1.92680986],\n",
       "       [ 0.72882184],\n",
       "       [ 0.62296124],\n",
       "       [-0.80584849],\n",
       "       [ 0.94202707],\n",
       "       [ 1.06721674],\n",
       "       [ 0.39764622]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.random.uniform(-1, 1, (20, 1))\n",
    "x = np.linalg.solve(A,b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27.021346300742454"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = b - A*x\n",
    "r = np.linalg.norm(v)\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residual of a linear system is the norm of b - A*x,\n",
    "but something went wrong, because the residual is too large!\n",
    "What is the problem?  Let us look at the type of A, b, and x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type(A) : <class 'numpy.ndarray'>\n",
      "type(b) : <class 'numpy.ndarray'>\n",
      "type(x) : <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "print('type(A) :', type(A))\n",
    "print('type(b) :', type(b))\n",
    "print('type(x) :', type(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the operator arithmetic to work properly,\n",
    "we must convert to the proper matrix types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<class 'numpy.matrix'>, <class 'numpy.matrix'>, <class 'numpy.matrix'>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mA = np.matrix(A)\n",
    "mb = np.matrix(b)\n",
    "mx = np.matrix(x)\n",
    "(type(mA), type(mb), type(mx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the residual correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.3440770491716685e-15"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = mb - mA*mx\n",
    "r = np.linalg.norm(v)\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Numerical Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most functions do not have symbolic antiderivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\int_{0}^{1} e^{\\sin\\left(x\\right)}\\,{d x}</script></html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\int_{0}^{1} e^{\\sin\\left(x\\right)}\\,{d x}$$"
      ],
      "text/plain": [
       "integrate(e^sin(x), x, 0, 1)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = var('x')\n",
    "f = exp(sin(x))\n",
    "a = integral(f, x, 0, 1, hold=True)\n",
    "a.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute a numerical approximation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6318696084180513"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.n()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see how we got to this ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_168/2732729778.py:2: UserWarning: get_systems() requires Cython profiling to be enabled, otherwise the results will be very unreliable. Rebuild Sage with the environment variable 'SAGE_PROFILE=yes' to enable profiling.\n",
      "  get_systems(\"a.n()\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['ginac']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sage.misc.citation import get_systems\n",
    "get_systems(\"a.n()\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GiNAC is a C++ library is it is not a CAS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes we may want more diagnostic information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(x) = exp(sin(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.6318696084180513,\n",
       " 1.8117392124517587e-14,\n",
       " {'neval': 21,\n",
       "  'last': 1,\n",
       "  'iord': array([          1,       21931,           0,           0,   -17793120,\n",
       "               32687,   -11053888,       32687,  1848425598,   179447140,\n",
       "           -17772848,       32687,   -11053760,       32687, -1094302173,\n",
       "          -925536350,   -17772752,       32687,   -11053696,       32687,\n",
       "          -656120084, -1109646371,   -17772656,       32687,   -11053632,\n",
       "               32687, -1863383552,   851750764,   -17772560,       32687,\n",
       "           -11053568,       32687,  -358732886,  1988466565,   -17772464,\n",
       "               32687,   -11053504,       32687, -1839087321, -1441984478,\n",
       "           -17760208,       32687,   -11053440,       32687, -1447038753,\n",
       "           190633445,   -17772368,       32687,         336,           0],\n",
       "        dtype=int32),\n",
       "  'alist': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       "  'blist': array([1.00000000e+000, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304, 7.74860419e-304, 7.74860419e-304,\n",
       "         7.74860419e-304, 7.74860419e-304]),\n",
       "  'rlist': array([1.63186961e+000, 0.00000000e+000, 6.93649197e-310, 6.93649197e-310,\n",
       "         6.93639347e-310, 6.93649200e-310, 6.93649197e-310, 6.93639349e-310,\n",
       "         6.93645281e-310, 6.93639349e-310, 6.93649191e-310, 6.93639349e-310,\n",
       "         6.93649086e-310, 6.93649200e-310, 6.93649093e-310, 6.93639349e-310,\n",
       "         6.93639349e-310, 6.93649197e-310, 6.93639349e-310, 6.93649194e-310,\n",
       "         6.93649197e-310, 6.93639349e-310, 6.93649197e-310, 6.93649197e-310,\n",
       "         6.93649197e-310, 6.93639349e-310, 6.93639348e-310, 6.93649197e-310,\n",
       "         6.93649200e-310, 6.93649197e-310, 6.93649197e-310, 6.93639349e-310,\n",
       "         6.93639348e-310, 6.93649191e-310, 6.93649201e-310, 6.93639349e-310,\n",
       "         6.93639349e-310, 6.93645281e-310, 6.93639349e-310, 6.93645281e-310,\n",
       "         6.93639349e-310, 6.93639349e-310, 6.93645281e-310, 6.93639349e-310,\n",
       "         6.93649197e-310, 6.93649197e-310, 6.93649197e-310, 6.93639349e-310,\n",
       "         6.93639349e-310, 6.93649190e-310]),\n",
       "  'elist': array([1.81173921e-014, 0.00000000e+000, 6.93646174e-310, 6.93646174e-310,\n",
       "         6.93646174e-310, 6.93646174e-310, 6.93649197e-310, 6.93649200e-310,\n",
       "         6.93649187e-310, 6.93649187e-310, 6.93649200e-310, 6.93649197e-310,\n",
       "         6.93639358e-310, 6.93639349e-310, 6.93649200e-310, 6.93649197e-310,\n",
       "         6.93649197e-310, 6.93649197e-310, 6.93649197e-310, 6.93649200e-310,\n",
       "         6.93649197e-310, 6.93649200e-310, 6.93649201e-310, 6.93649201e-310,\n",
       "         6.93639347e-310, 6.93649191e-310, 6.93649201e-310, 6.93649197e-310,\n",
       "         6.93649200e-310, 6.93649200e-310, 6.93649200e-310, 6.93649190e-310,\n",
       "         6.93639349e-310, 6.93649038e-310, 6.93639349e-310, 6.93639349e-310,\n",
       "         6.93649201e-310, 6.93639347e-310, 6.93645289e-310, 6.93649061e-310,\n",
       "         6.93639349e-310, 6.93649200e-310, 6.93646173e-310, 6.93649200e-310,\n",
       "         6.93649093e-310, 6.93649200e-310, 6.93649200e-310, 6.93649200e-310,\n",
       "         6.93649200e-310, 6.93649096e-310])})"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = quad(func=f,a=0,b=1,full_output=1)\n",
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Rational Approximations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we use a combination of sympy and scipy.\n",
    "For example, the computation of a Pade approximation\n",
    "starts from a Taylor series.  We can compute Taylor series\n",
    "with sympy and Pade approximations with scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "reset()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import series, sin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.abc import x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<generator object Expr.series.<locals>.<genexpr> at 0x7faffef93990>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = series(f,x,0, n=None)\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x, -x**3/6, x**5/120, -x**7/5040, x**9/362880, -x**11/39916800]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "terms = [next(s) for _ in range(6)]\n",
    "terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{x^{11}}{39916800} + \\frac{x^{9}}{362880} - \\frac{x^{7}}{5040} + \\frac{x^{5}}{120} - \\frac{x^{3}}{6} + x$"
      ],
      "text/plain": [
       "-x**11/39916800 + x**9/362880 - x**7/5040 + x**5/120 - x**3/6 + x"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sterms = sum(terms)\n",
    "sterms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 0, -1/6, 0, 1/120]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = [sterms.coeff(x, k) for k in range(6)]\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute a Pade approximation, we must convert the list to an array of floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 0, -1/6, 0, 1/120] has type <class 'list'>\n"
     ]
    }
   ],
   "source": [
    "print(c, 'has type', type(c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.          1.          0.         -0.16666667  0.          0.00833333] has type <class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "nfc = np.array(c, dtype=float)\n",
    "print(nfc, 'has type', type(nfc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we import the pade command from scipy.interpolate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(poly1d([-0.11666667,  0.        ,  1.        ,  0.        ]),\n",
       " poly1d([0.05, 0.  , 1.  ]))"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.interpolate import pade\n",
    "p = pade(nfc, 2)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On return we get numpy polynomials, which have a nice string representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.poly1d'>\n",
      "numerator :\n",
      "          3\n",
      "-0.1167 x + 1 x\n",
      "denominator :\n",
      "       2\n",
      "0.05 x + 1\n"
     ]
    }
   ],
   "source": [
    "print(type(p[0]))\n",
    "print('numerator :\\n', p[0])\n",
    "print('denominator :\\n', p[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numpy polynomials are callable, so plotting is no problem\n",
    "and we compare with the plot of the actual sin(x) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD/CAYAAAD12nFYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7NUlEQVR4nO3dd3yNZxvA8d9JJLEyECL23qOkRmxVMUp5tUqrNi01i1ZRmxr1VouYLTqMDqu1tWLHFiv2SowgRmJEIsnz/nFXXgeJRM45zxnX9/M5nzbn3Oc81yPkyr2u26BpGkIIIcQTTnoHIIQQwrpIYhBCCGFEEoMQQggjkhiEEEIYkcQghBDCiCQGIYQQRiQxCCGEMCKJQQghhBFJDEIABsXDYDAY9I5FCL1lSENb2SIt7FZUVBSenp5ERUXpHYoQ5pSqX3ykxyCszrZt22jevDl58uTBYDCwcuXKl75n69at+Pn5kTFjRooUKcLs2bPNH6gQdkoSg7A6Dx48oGLFisyYMSNV7S9cuEDTpk2pXbs2hw4dYujQofTt25dly5aZOVIh7JMhDUX0ZChJWJzBYGDFihW0bNky2TaDBw/mzz//5MSJE0nP9ejRg8OHDxMcHJyq60RHRycNJXl4eKQ3bCGsVaqGktIyxyCEVQoODiYgIMDouUaNGvHDDz/w+PFjXFxcnntPbGwssbGxSV9HR0ebPU6HER5O7NbdnN55kxsXH3LryiOib8fjrMXj6hSPt9s98vomUryYhlshXyhbFqpWhQIF9I5c/EsSg7B5ERER+Pj4GD3n4+NDfHw8kZGR+Pr6PveeCRMmMHr0aEuFaN9iYri6eAur515h9zF3Dj0swXFa8BjX5N9zDpx3xFPM+QI1ErbTiM94M88JcjStBv/5DzRoAG5ulrsHYUQSg7ALz64yfTJEmtzq0yFDhjBgwICkr6Ojo8mfP7/5ArQ3msaddbtZOPwcS0JKsS+xCc7EUyHXdfz8H9P9zRgq1HLF1xdy5AAPD9A0iI2FyEgIC4PjxzNw9GhxtvxdmAWnumC4moj/L4f48Pvv+dC9O+6d34VPPoGSJfW+W4cjiUHYvNy5cxMREWH03I0bN8iQIQM5cuR44Xvc3Nxwk99I0y4xkSNT/yFwQhS/3GrCY/xoXuYcfTpc5a3ueciePW+Kb8+cWY0YFSgAtWo9eTYDly/Dxo1OLFtWmT4bZjL08Tf0mjeLPtPq4BPwGgwYAAEBINtMLEJWJQmb5+/vz6ZNm4ye27hxI6+//voL5xfEqzkyfSvNvbZTcVBDVt+vxxcdrhF2JQPLjpem/eA8ZM/+6p+dLx906QJr1hi4cMFA556Z+NbpUwq6XKVXSDduNv4Q3nwTDhww3Q2JZEliEFbn/v37hISEEBISAqjlqCEhIYSFhQFqGKhDhw5J7Xv06MGlS5cYMGAAJ06cYP78+fzwww8MGjRIj/DtzrV9l+mQP4jX+tbmRHwxfhlxmov3vBn+YzFy5zH9j5D8+eGbbyAszMCXI5xZHPcuJbNeYXZoHRJerwrt2qmxKGE+mqal9iGERQQFBWmo5dFGj44dO2qapmkdO3bU6tata/SeLVu2aJUqVdJcXV21QoUKabNmzUrTNaOiojRAi4qKMtFd2L7E2Dht3rvrNU/uaDmdbmozux3Q4mITLR7H9eua1rmzpoGm+RW8oe3L0UjTsmbVtJkzNS0hweLx2LhU/byXfQxCIPsYnnV+3Sm6v3eXzfer0an0Hv67tjTZC+n75xIcrOaijx3TGF95OYP2tsapfj34/nsoUkTX2GyIlMQQQqSRprGi9z9Uapqb87H52DDjDAtCq+meFAD8/WHvXhgwwMDgve/QxO8mEWfuQaVKsGKF3uHZFUkMQggA4qMe8Fm5dbQKbEDDQmc5HJ6dgF7F9Q7LiIsLTJoEGzbA4cs5qBy/h/1+H0OrVjBwIDx+rHeIdkESgxCCa7sv8UaeE0wNDeCbD/bz+3k/PHwy6R1WsgIC4NAhyF/AiTq7J/FHx79g2jR44w24eVPv8GyeJAYhHNz+BUepXDMj52LzseXny3y66HWb2C7g6wtbtkCLFgZa/9iMSd3PwunTUK0ahIbqHZ5Nk8QgHFpgYCBlypShSpUqeoeii7+G7KJulyIUzBzJwSMu1PqwkN4hpUmmTLB4MQwfDl/MKsgXrU6jZc6iJiS2btU7PJslq5KEwDFXJS3q8g8dF9Tl7XwHWXSkApmyZdQ7pHSZOlVtkO7ZNY4ZF5vhtHM7/P47NGumd2jWRFYlCSFebF6bv2m/oD7tSx/g9/Ov23xSAPj0U7VydfZ8V3oXXYfWuIkqyLd4sd6h2RyplSSEg/m22d98uuZNelXexbS9/jg528CEQip17ar+262bM1kH/sGkD7pg+PBDiIn5/4vipSQxCOFAJjX8my/+fpPPa+5k4rYaGJzsJyk80bUr3L8P/fs74T56PsN7ZIbu3dVa16dKqYjkSWIQwkHMarWJL/5uyPD6Oxj9Ty2bWHn0qvr1g3v3YPhwJ7y+m0Gfro+hc2d1xkObNnqHZ/UkMQjhAJZ0+4deKxrQr8pOu08KTwwbBnfuQL/+ThRcPoe3Y2NVAb6MGaFFC73Ds2qSGISwc2sGb6PDD3XoUGof3wTXcIikAOrohq+/hosX4f12TmwLWoDfw4fQti388w/UqKF3iFZLlqsKgf0uV90+PYSAviVplC+UP85VIoOr4y1EfPgQ6tdXlbr373hE3i6N4Ngx2LkTSpXSOzxLS9WvBZIYhMA+E8PJdReo/lZ2KnueZ+2lsmT0SOEMZjt3/Tr4+UHBghC0/A6uDWqrSYjgYMiTR+/wLEn2MQjhqG6diqTZ207kdbnJipDCDp0UAHx81F63ffvgs6+ywbp1kJAATZvCgwd6h2d1JDEIh2aPJTHiomJoVSWcqIQsrP47I54FvfQOySr4+6vd0dOmweLt+VVyOHsWOnaExES9w7MqMpQkBPYzlKQlanQtsZ1F56qxed55anYrrXdIVkXT1FaG5cthzx4od3al2h09ejSMGKF3eJYgQ0lCOJqp7+5kwbk6/NDzgCSFFzAYYM4cKFZMHeFwr0FLGDsWRo5U2UIAkhiEsBvbA4/w+YrqfO73Dx/OlKWYycmcGZYtg2vXoH9/1IaH1q2hfXs4fFjv8KyCDCUJge0PJV0/HEGlylDC/Rp/XytHhkwueodk9ebPV+Uzli2DVo0fQs2aaqXSwYNgg38HUkmGkoRwBPEP43i/zhU0DCzdnk+SQip17gwtW8JHH8G1qMzwxx/q9Lfu3dVkhAOTxCCEjRtZJ4it0a+xdHokucvn1Dscm2EwwLx5qrZely6gFSmquhG//QazZukdnq4kMQhhw1YP2clXBxrxVcu91P2krN7h2Bxvb5UL1q+HmTOBd96BPn3U4Q4HD+odnm5kjkEIbHOO4cKOK1Suk4U6uc+wIvx1uzpXwdJ694YffoCjR6FY/lioVUtV4DtwADw99Q7PlKQkhhCpZWuJ4XFMPLV9TnHjoTsHzniQrbCX3iHZtAcPoHx5KFRI1dczXDgPlStDw4ZqaMl+Kg/K5LMQ9mpc4x3sv1eSJYG3JSmYQJYsMHcuBAWpngNFisCCBWpCes4cvcOzOEkMwqHZYkmMXTNDGLetNiPq76Dax6/pHY7dePNN6NQJBg2Cq1dRO6J79ICBA+H0ab3DsygZShIC2xlKig6P4rXCUfhmiWLrjTJkcHPWOyS7cvs2lCmjjmpYvhw1xlSpEnh5qTLdLja/FFiGkoSwN33qHyMywYufV2eXpGAG2bOrInsrVsCaNagxpl9+USuUxo3TOzyLkcQghI1YNvQAP52ryfTOhyhSO6/e4dit1q3VsFK/fvDoEVC1qiqwN26cOr/BAchQkhBY/1BS5Jk7lCmVQM2cp1l+1R+Dk92skrFKJ05AhQowapQqpUR8/P+XsIaEQKZMOkf4ymQoSQh70afhSRI0J2atKyxJwQJKl1Z73MaPV0eCkiGDWqV06ZLKFnZOEoMQVm754D0sveTP9I+Pk7uSr97hOIzhw9Wc88CB/z5RurRKClOmqKPg7JgMJQmB9Q4l3Tp9izKlEvH3OceKK9Wkt2BhixdDu3awaZOadyA+HqpXV5MPBw6Am5veIaaVDCUJYev6NjzBYzIwa30RSQo6eP99qF1blU+Ki0MNKc2fD6dOqXEmOyWJQQgrtWroHhaH1WJaz5P4VsyldzgOyWCAGTPU/rbAwH+frFBBzUhPmKAmou2QDCUJgfUNJUVfuUfpAg+onOMSf0ZUld6Cznr0UCWTzp2DbNlQ3YfXXwdnZ9i715Y2vslQkhAvY60lMb586xBRiVmZuTKPJAUrMGqUygVJo0eurmqV0tGj8PXXeoZmFtJjEALr6jHs+zGUap1KMaXZVgb8VV/XWMT/jRmjEsPJk1C48L9Pfv45TJ8Ox4+rwnvWT8puC5Fa1pIY4h/FUzX7GTQM7LtdjAwZM+gWizD24AEULw716qnVSklPlimjHmvX2kJ5bhlKEsLWzGiznZCYksyZmShJwcpkyQJjx8KSJbB//1NPTp+ujoD74w9d4zMl6TEIgXX0GC4Hh1O6hhcdKxxixuE6usQgUpaQoBYl5csHGzY89cJ//gN79qhxJitYvJAC6TEIYTM0jb7/Ccfd6QHj11TSOxqRDGdnNdewcSNs2/bUC999B9HRaru0HZDEIIQV+HPYHlZcr8G3n4bhmc9d73BEClq1Uqd+DhsGSQMuBQrA6NFq08OBA7rGZwqSGIRVmjlzJoULFyZjxoz4+fmxffv2ZNtu2bIFg8Hw3OPkyZMWjPjVPYiMoc/k/DT23kfryda1bFY8z2BQFbh37HhmOKlvXyhXTm16SEjQLT5TkMQgrM6vv/5K//79GTZsGIcOHaJ27do0adKEsLCwFN936tQprl27lvQoXry4hSJOn4nv7ON6Qg5mLPGWPQs2onFjdcrbl18+1WtwcYHZs9XM9Pz5usaXbpqmpfYhhEVUrVpV69Gjh9FzpUqV0r744osXtg8KCtIA7c6dO698zaioKA3QoqKiXvkzXsX5beGaGzHaMP9/LHpdkX5BQZoGmrZ8+TMvtG+vad7emnb7th5hvUyqft5Lj0FYlbi4OA4cOEBAQIDR8wEBAezatSvF91aqVAlfX18aNGhAUFBQim1jY2OJjo42euhhUNvLeDvdYciKqrpcX7y6evVUxdXhw58ZOZo4UVVfHT1ar9DSTRKDsCqRkZEkJCTg4+Nj9LyPjw8REREvfI+vry9z585l2bJlLF++nJIlS9KgQQO2GS0bMTZhwgQ8PT2THvnz5zfpfaTGP5MPsPxqdSb3OE8Wn6wWv75Iv7Fj1abnFSueejJPHpUtZsxQL9og2ccgrMrVq1fJmzcvu3btwt/fP+n58ePH8/PPP6d6Qrl58+YYDAb+/PPPF74eGxtLbGxs0tfR0dHkz5/fYvsY4h/G8Vq2S3hljGH7nfIyt2DDGjaEmzfh0KGnNj7HxqqJ6IIF1WEO1rMjWvYxCNvj7e2Ns7Pzc72DGzduPNeLSEn16tU5c+ZMsq+7ubnh4eFh9LCkWe12EBpXlOlzM0pSsHFffgmHD8OaNU896eYG334L//wDK1fqFNmrk8QgrIqrqyt+fn5s2rTJ6PlNmzZRo0aNVH/OoUOH8PW1zmMwbx67zoiVlehWZheV2pTQOxyRTnXqQK1aagmr0QDMW29BkyYwYADExOgW36uQxCCszoABA/j++++ZP38+J06c4NNPPyUsLIwePXoAMGTIEDp06JDU/ttvv2XlypWcOXOG48ePM2TIEJYtW0bv3r31uoUUDX8nFM3gxPgVZfUORZiAwaB6DXv2qA6CkalT4coV+O9/dYntVUmVLmF12rRpw61btxgzZgzXrl2jXLlyrF27loIFCwJw7do1oz0NcXFxDBo0iCtXrpApUybKli3LmjVraNq0qV63kKxDi44z93RdprYOJmeJmnqHI0wkIECd2zNu3L9nQz9RsiT07w9ffQUdO4IOixxehUw+C4FliuhpiRp1sh3ldmxmQu4UwiWT/F5mT1atgpYtVQ2l2rWfeiE6GkqUUOtbly7VKbokMvkshDX57fP97IiuwHcjbktSsEPNm0P58k+d8vaEh4c6H/rXX2H3bl1iSyvpMQiB+XsMsffiKJX9OuW8rvDXjWrWtHxRmNCvv0LbtuoYaKPTYhMSwM9Pnd+wY4ee33/pMQhhLaa330N4vC9ff59NkoIde/ddNWo0ceIzLzg7w5QpsGsXLFumS2xpIT0GITBvj+HWubsULW6gXakDBIa+YdLPFtZnzhzo2RPOnIGiRZ958a231GE+oaFqr4PlSY9BiJcJDAykTJkyVKlivnLXY1ofJVEzMPJXWZ7qCDp0gBw51P6253z9NVy8CIGBFo4qbaTHIATm6zGc+SeMMm/6MubN7QzZJL0FRzFqlMoB4eGQPfszL/booSYjzp17wYtmJz0GIfT2Rcdr+DrfoP/S6nqHIizok08gMVEdz/Cc0aMhPl5V4LNSkhiEMJPtM4+y/Eo1vup6nkw5MusdjrCgXLnUkNL06aqenhEfHxgyRA0nnT2rS3wvI0NJQmD6oaTEBI3qnidISDSwL7okThnkdzBHc+oUlCqlDnPr3PmZF2Ni1PKlqlUtvUpJhpKE0MvvA4LZ96AMU8Y9kqTgoEqWVJve/vvfZ4rrAWTKpHbCLV9ulZvepMcgBKbtMTx+EEdpr2uUzHaDNTfMt9pJWL9t26BuXVi3Tp0TbSQhAV57TS1hCgqy1P4W6TEIoYd5XYM5H5+fCbO89A5F6Kx2bVVcb8qUF7zo7KxKZWzdCuvXWzy2lKS7x3DypFqW+1w2FMKGmKrHcD/iPsXyPiSg0Bl+OifVU4Wqm/f++3D0qDrUzYimqQMd7t2DgwfByey/q1umxxAYCN26qaVZQji6bz/cz51ET8b8WFDvUISVaNUKfH2T2dNmMMCkSeoIuCVLLB5bctKdGN57T51DYYXzJ0JY1M2Tt5j8T2U+qbSbQrXy6R2OsBKurvDxx/DTT3D37gsa1KgBb78Nw4dDXJylw3uhdCeGmjVVNvz9d1OEI4RlmbIkxlcfHANg2GIpfSGMffwxPH4MCxcm0+Crr+DSJVVoyQqYZFVSnz6wYgWEhVliiEwI00vvHMPFHZcpWTsnw9/YxZf/1DdDhMLWffAB7Nun9je88Odk586wZo0qleHubq4wLLcqSYaThKMb0fES2Zyi+HSxLE8VL9a7t9rovHFjMg1Gj1anvX3zjUXjehGTJAYZThKO7MiyM/xy3p+R74aSxSer3uEIK+XvD5UqqTIZL1SgAPTqpda23rhh0dieZZLE4OQE77yjEoOsThKOZugndymaIYxu82voHYqwYgaDGnZfty6FEklDh6ofqM+dD2pZJpsRkOEk4Yh2zT7CmhtVGNfjMi5ZXPUOR1i5tm0hWzaYNSuZBjlywOefqwYXLlg0tqeZrCRGYiLkywdt2sDUqekPTAhLeqXJZ02jQbaDRMa6c+heMamJJFJl8GCYOxcuX1ZHQD/nwQMoVgzefBN+/tnUl7dsSQwZThKOZvPEvWyO8mPs5/ckKYhU69lTzTEvWpRMgyxZYMQI1eDIEYvG9oRJi+ht3652d+/cqfZsCGEr0tpj0BI1ankc4bHBhT1RpTE4WaQAmrATLVqobQuHDiVTO+/xYyhdGsqWhVWrTHlpyxfRe7I66bffTPmpQlif9SOD2fWgIuNGPJakINLs449VFYz9+5Np4OKizgf980/Yu9eSoQFmKLstm92ELUpLj0FLSOT1rCfJkjGerbcrWKhasrAnCQlQuDA0agTz5qXQqHx5yJ8fNmww1aX1Kbstq5OEvVv5+U4OPirDuPFOkhTEK3F2hq5dVd28e/dSaDRmjNoRt22bReMzeY/hyeqk996Db7995biEsIjAwEACAwNJSEjg9OnTL+0xJMTGU9HjPL7uD9gUWcmCkQp7Ex4OhQqplakffZRMo8REdaBD1qzq3Ib0/yaSqg8wywluMpwkbE1qh5KW9NzGB7PrELzgJNU7lbJghMIeNWsG16+rGkrJWrNGNdywAQIC0ntJ/RKDrE4StiY1iSH+YRxlPK9Q0vsWf1173cIRCnv0559qhdKBA1C5cjKNNE39IE1IgD170ttr0O9oT1mdJOzRzz13cSa+MGOmeekdirATTZtCnjwpTECDSgTjxqluxV9/WSQus/QYQIaThG15WY8hLvoRJbPf4PU81/g9rJoOEQp7NXw4fPcdXLuWzE7oJ954A27dUpsfXv2Hqn49BlClMa5cgR07zHUFISxnfrddXErIx+jZPnqHIuxM165w/z78+utLGo4dq3ZCW6CMtdl6DImJap1u06YpFIwSwkqk1GOIiXxAMZ9o3ih8gZ/PyqSZML3GjdWxny9d5t+0KZw/D8eOQYYMr3IpfXsMTk7w/vtqnsFKjjEV4pXM6byb64k5GTkvv96hCDv10UdqXvmlpZHGjlVHwCVbaMk0zDr6//77cPs2bNpkzqsIYT73r0YzYU15OpfaTbH6khiEeTRvDj4+L5mEBvDzg1atVLkMM/7GbdbEUKEClCkDixeb8ypCmM+Mjvu4q3kyfGFRvUMRdszFBTp1Uh2B2NiXNB49WlXgmz/fbPGYNTEYDOoA7FWrVIlxIWxJ1MU7TP67Mt0r7qVANV+9wxF2rmNHuHMnFStSy5VTwzFjx0JMjFliMftC0vffV0nBQstvhUiTwMBAypQpQ5UqVZ57bWr7g8SQkaE/yg5nYX6lS0O1arBwYSoajxqltkzPmWOWWMy2Kulp1atDrlxql58Q1ujZVUm3TkVSuJQrH1UJYcreOnqHJxzE7NnQu7eqo+T7sk5qt27qh+r586qWUurouyrpaR98AOvXq4loIWzB1+2PkIgTg38up3cowoG0aaNWoaZq0dHw4WqN6/TpJo/DIonhvfdUmY9lyyxxNSHSJyIkgmn7qtO/9kFylsyudzjCgWTLBi1bquGklw7mFCyo1rlOnqwShAlZJDHkzg0NGsjqJGEbJnY6gavhMQN/rqh3KMIBdeoEx4+rwnovNWwYPHoEU6eaNAaLVTH64ANVTvzyZUtdUYi0u7LvCrMP+zOoQQjZCnrqHY5wQA0bqsJ6qZqE9vVVkxLffAORkSaLwWKJoVUrcHOTXoOwblN6nier4QH9fvLTOxThoJydoX179bPypXsaAAYPVv/9+muTxWCxxODhocbOfvwxFWNnQujkxzPV+aLZcdx9U73KQwiTS/WeBgBvb+jbF2bMUEtYTcCiBbE7dIDQUFU1VghrlMNwh08WPL+nQQhLStOeBoCBA9VypkmTTHL99CcGTYMbN1LVtGFDVQ/kp5/SfVUhTOrMhnMADGxxlsw5MukcjRBqEnr9enVOw0tlzw4DBqhS1levpvvaqdrgZjAYDFFRUYkvfPHTT2H/fti2LVVHzg0dCkuXqgKBLi5pjlcIk4iNjSX2qQHcnq/vYfX1dzkXeg7vvN46RiaEcucOlCgBI0aog89eKipKFahr3RqmTHlhE09PT0/gnvaSH/ypTQweQFQqQhNCCGHdPDVNi06pwSv3GKpUqcK+ffvUFy1bqv7Orl1qSv0p0dHR5M+fn/Dw8KQDUGrUgOLF1UT0s4w+NxXS0t5cbV90j3rEYa7PTuv9mSuOtLZNqf3TPYaPKp/g+C13whL9CQ0NJW/evBaLw5Jt5e+pZeJIa9uU2q9eDe3aqbMaSpVK//cwtT2GVB0B9KIPcXZ2/n9gEyaogkjr1qkNCy/g4eGR1L5zZ7UvIyFB7fRL9nNTIS3tzdX2iafvUa84zPnZqb0/c8Zh6j+P/T8eZ93Nhkz7cA19fwF3d3eTfw/T2l7+nqbvs+3p7+k770CvXqpCddWq/3/+Vb+HL+spPPHKk8+9evX6/xfVqqmTJkaOhPj4l773gw/g8eMXH11q9LlpjUOntmllzjjs/R5N/ecxfOBDSrme4+0xaauJZC3fF1v8Hpr7s60hDlP9ebi5qZJCixap45LT6pX/7DRNS+0jZYcOaRpo2g8/GD0dFRWlAVpUVJTR840ba1rNmi/9VJuQ3D3aC3u9v+0zQjTQtN8G7NLCw8M1QAsPD9c7LLOx1+/jE/Z6f9u2qR+t27aZ5B5T9fPedPsYXnsN3n0XxowxOnLOzc2NkSNH4ubmZtS8QwfYuRPOnTNZBLpJ7h7thV3en6YxfFgiFTOe4p1J1ZLuza7u8Rl2+X18ir3eX82aUKAA/PKL5e7RtOcxhIaq04UCA6FnzxSbPnyoiusNGKDOnBDCkjZ/fYAGn/vx5/C9NB9T9bnzGISwJkOHqrMarl1Tw0vpkKrzGEx/UE/79rB5M5w9C5lS3ijUtatqeu4cOFl0D7ZwZFqiRk3PoyRozuyOLoPBySCJQVi10FAoWxZWrFCLQNNBp4N6Ro5U9Tpmz35p086d4eJFCAoyeRRCJGvdmH0E36/A2GGPMDil6t+JELoqU0aN1v/yi2WuZ56jPbt3V+urXnLknKapmiCVK0vVVWEZWkIifu6nyOoSx9Y7FZISg/QYhLX773/VkNL16+Dl9cofo+PRnl9+qU4UmjEjxWYGA3TpAsuXq+3fQpjbii/2cCimNOPGIb0FYVPatlXL/FetMv+1zJMYnj5yLur/lTTefvttChQoQMaMGfH19aV9+/YEBEQQH5/KM06t3MWLF+natSuFCxcmU6ZMFC1alJEjRxL31CotezB+/Hhq1KhB5syZ8UrHry6WlhCXwIjp3jTMfoA6fdTpbIGBgZQpU4YqVey3ouq2bdto3rw5efLkwWAwsHLlSr1DMqkJEyZQpUoV3N3dyZUrFy1btuTUqVN6h2Uys2bNokKFCpQu7YGT0w769Qtm3bp1Zr2m+aZ8hw6FmBijI+fq16/Pb7/9xqlTp1i2bBnnzp3jk09a0awZ/PCD2SKxmJMnT5KYmMicOXM4fvw4U6dOZfbs2QwdOlTv0EwqLi6O1q1b0/MlK8+szdJ+wRyPLc7YrzMmPderVy9CQ0PTVL7A1jx48ICKFSsy4yU9eFu1detWevXqxe7du9m0aRPx8fEEBATw4MEDvUMziXz58jFx4kT279/PsGHFiY6uyttvd+H48ePmu2hqNzy80laKgQM1zd1d0yIjX/jyqlWrNIPBoC1f/lgDTTtw4JWuYtUmT56sFS5cWO8wzGLBggWap6en3mGkStyDOK1ohova27l3v/B1e90c9SxAW7Fihd5hmNWNGzc0QNu6daveoZjc9eua5uSkaZkz99e+//77V/kIC29we5HBg9U+7hccOXf79m0WLVpEjRo1aN48A76+9tFreFZUVBTZs2fXOwyHt+CjYM7H52fsDPle2Luof4ev7fHfXY4cCZQtG0FMTAv8/f3Ndh3zJoacOaFfP5g+PenIucGDB5MlSxZy5MhBWFgYq1atIkMGdZTdokVq9MlenDt3junTp9OjRw+9Q3Foj+4+YsyS4rQtGEyFd4rrHY4wI03TGDBgALVq1aJcubTVv7JmR48eJWvWrLi5uXH27DigLl5eZcx2vXQnhlGjRmEwGJJ9ZPvqK+KdnGDiRAA+++wzDh06xMaNG3F2dqZDhw5omkaXLmqeevnydN+Tyb3sHg0GA/v37zd6z9WrV2ncuDGtW7emW7duOkWeeq9yj7ZiVqfdRCTmZPS8l5fTFratd+/eHDlyhCVLlugdikmVLFmSkJAQdu/ezccf50TTHjNjRoTZrpfufQyRkZFERkam+MaiS5bgMmmS2g2dL1/S85cvXyZ//vzs2rULf39/6tZVxzls3pzakCwjNfdYqFAhMmZUk5pXr16lfv36VKtWjYULF+JkA9u603qPAAsXLqR///7cvXvXzNG9unvX7lM0bwwtip9g3qk6ybZzlH0MBoOBFStW0DKd22etUZ8+fVi5ciXbtm2jcOHCeodjVjlz7sLZ2YeIiKJpfWuq1min6jyGlHh7e+Pt/ZKjEAcOVPWTxo9XZ5L+60lSenJgSteuakjp3Dkomub7NZ9U3eO/rly5Qv369fHz82PBggU2kRQgbfdoS75rv58ozZ/hC63oL5QwKU3T6NOnDytWrGDLli12nxQAcuUKIjR0GBcvQqFCpv98y/zU8vAg7P33SZg7l9A1a7h06RJBQUF88MEHFC1aNGkSpXVrtaNvzhyLRGVyV69epV69euTPn58pU6Zw8+ZNIiIiiIgwX5dPD2FhYYSEhBAWFkZCQgIhISGEhIRw//59vUMzcvv8Xab8U4merwVTwN9xh5Hu37+f9D0CuHDhQtL3zx706tWLX375hcWLF+Pu7p70by7GTiYshw4dyvbt27l48SJHjx5l2LBhnDgxCTe3BH77zUwXTe3ypVddXvXE0T17tFuurtoiV1fNzc1NK1SokNajRw/t8uXLRu3699e0HDk0LSYmvVe0vAULFmioIbfnHvakY8eOL7zHoKAgvUMz8oV/kJaZ+1rE0RsvbWvPy1WDgoJe+P3q2LGj3qGZRHL/5hYsWKB3aCbRpUsXrWDBgpqrq6uWM2dOrUGDBtrGjRu11q01rVKlNH9cqn7em6dWUnKmTYNPP4UTJ6BEiRc2OXVKnW3688/w4YfpvqJwUBFHb1K0Qmb6+e/lq131X9reUeYYhP1YtkwdgXPyJJQsmeq36VgrKTkffQR58qR4AEPJkvDGG0ZTEUKk2YQOobjwmM9+eS3Fdo5QEkPYp6ZNVY3SX381/WdbtscAMHcu9OgBISFQocILm/zxh5pvOHw42SZCJCts91WK++dgRINdDPv75b0FkB6DsE3t28PBg5CG6hhW2GMAdQhD0aIwbFiyTVq0AF9f6TWIVzOm0zk8DdH0+/l1vUMRwqxat1aH+ISGmvZzLZ8YXFxg7FhYvRp27Ei2Sbdu6lCKe/csHJ+waac3XWLhKX+GNj9GVl93vcMRwqwCAsDdXY2ymJLlh5JA1U/y81MDZNu2qYMZnhEertbnzpjx0uOjhUjyfqFgdoQX5ExkNjJmS/lo2afJUJKwVe3awdGjcORIqppb6VASqAOeJ0xQPYZk6ornzw/Nm6vhpNTnLuHIDi87y9JL/oxoezpNSUEIW/buuyoxmPIICv225TZqBHXqwJAhqgfxAj17qhvetcvCsQmbNLzXLYpluEinuTX0DkUIi2ncGLJkMe1wkn6JwWBQvYYjR2Dp0hc2adhQzVPLJLR4meBZIfx1vRqje1zDJYur3uEIYTGZMkGzZqZNDPrMMTytRQs4dkxtenN9/h/011+rI6QvX1ZVvIV4lpaoUcfrMPfiM3EwujhOGdL++47MMQhb9mSJ/5kzUKxYik2teI7haePHw4ULyZ7S07mzmpL4/nsLxyVsxurhe9hx7zUmDX/wSklBCFvXtClkzmy6XoP+PQaADh1g0yZVljtLlude7tYN1q9X+cPFxWxRCBsU/yieip4XyZ3lHn/fqvSiBW6pIj0GYetat1Y/I19ybIqN9BgARo+GW7dULaUX6NcPrlxRtUGEeNpPPXcRGleMSd9lfKWkICUxhL1o3RoOHIDz59P/WdbRYwDo2xd++kkdxpAjx3MvN2gADx9CcLBZoxA2JOZ2DMVz3qVW3gssDUvfSiTpMQhbd/++mocdMwY++yzZZjbUYwA1w5yYCOPGvfDlfv1g927Ys8fCcQmrNe3DPVxP9GbcAsc9a0GIJ7JmhU8+gWzZ0v9Z1tNjADURPXq0WqH0zBFuCQmqUne1arB4sdkjEVbu9vm7FClqoH35EKYfqZvuz5Meg3AQNtZjAHVWQ65cMHTocy85O0OfPvD772q+QTi2r9oeIQEnhi8urXcoQtgd60oMmTOrAnu//fbCMaMuXdRmDtnw5tguBV9l+r5qfFZ3H7nK5dI7HCHsjnUNJYEaM6pUSR3+vHXrcwX2+vdXp7uFhb1wZatwAB2L72TDueKcvZyRrHlMM+wjQ0nCQdjgUBKoMaPJk2H7dvjzz+de7t8foqJg/nzLhyb0d+SP0/x81p+R7x43WVIQQhizvh4DqHKqAQGq9vbRo8/tavvwQ1WY9cwZ2fDmUDSNpjn3cjY6F8fv5sMls+m++dJjEA7CRnsMoIaPJk+G06dfWAvjs8/g0iU1FSEcx8axe1h3qxoTBt4yaVIQQhizzh7DEx06wIYNqlSGu/FpXE2awNWr6ujoVy2DIGxHfMxjXvO6QI7MMWy5VQGDk2m/6dJjEA7ChnsMT4wbpyYUvv76uZcGD1YVuzds0CEuYXHzO2/neFwJvpmZyaRJQUpiCPE86+4xAHzxBUyfriYU8uT5fzCa2uyWNSts3qxLZMJCoi/doXjhxzQudpYfT5vnEB7pMQgHYQc9BlAnvGXKBMOGGT1tMKheQ1AQ7NunU2zCIia8e4B7WlbGL0250LwQwjSsPzF4eqpNbwsXPpcBWraE4sVh0iRdIhMWcHHzeabur8XnbxwgX2XZzCaEJVj/UBJAfDxUrqzGjXbuNJptnjsXevRQB2EXL65bhMJM2ubbzraIEpyJ8CCLdyazXUeGkoSDsJOhJIAMGeDbb1XN7SVLjF7q0EGVV5oyRZ/QhPkEzzjAr1dq81W3C2ZNCkIIY7bRY3iiVSs1nHTypFE9jIkTYeRIdZRDvnw6xidMJvFxAjW8jvPY4Mq+qJI4OZt3TbL0GISDsKMewxNTpsCNG2rz21M++USNMslcg+27c+cO7du3p2GWj9nzsAIVyvxA9L2oFN/TqVMnDAaD0aN69eoWilgI+2NbiaFIERg4UCWGS5eSnvbwgAED1HyDlOS2bR988AGhe09zLP4rWvlsZn/Metq3b//S9zVu3Jhr164lPdauXWuBaIWwT7aVGEAtX82WTa1VfUqfPmp0SXoNtuvEiROsX7+eSoYxPNQyMW1VCebNm8fq1as5depUiu91c3Mjd+7cSY/s2bNbKGoh7I/tJQZ3d5gwAX79VVVg/ZeHh+pMSK/BdgUHB1M2UzUWnHqTEY32kLdaPqpXr46npye7du1K8b1btmwhV65clChRgu7du3Pjxo0U28fGxhIdHW30EEIotpcYANq3hypVoG9ftZT1X336qLN+pNdgm65dvUaG2CmUcLlAv99qJT2fK1cuIiIikn1fkyZNWLRoEZs3b+a///0v+/bt44033iA2NjbZ90yYMAFPT8+kR/78+U16L0LYMttMDE5OMGMGHD4MM2cmPS1zDdZp1KhRz00OP/vYv38/N9d4cjixFtNH3cbVI2PS+zVNw5BCpcQ2bdrw1ltvUa5cOZo3b866des4ffo0a9asSfY9Q4YMISoqKukRHh5u0nsWwpbZ1nLVZ/XsCYsWqd1tvr6AqrlXuLA6s2HaNJ3jEwBERkYSGRmZYpvsrt5UKBZPMXawI/Fdo9e8vLyYOnUqnTt3TvU1ixcvTrdu3Rj8zFxUcmS5qnAQqVqumsHcUZjVV1/B8uWqm/DvxjdPT/XluHEwaBAUKKBzjAJvb2+8vb1TbDOoylbuaa9zkQHs3VuAqlWrArBnzx6ioqKoUSP1xfNu3bpFeHg4vv/+siCESBvbHEp6Ils2tbdh6VLYtCnp6X791LDSqFH6hSZSL/Svc3y3vwbDGuyhfOOydO/end27d7N79266d+9Os2bNKFmyZFL7UqVKsWLFCgDu37/PoEGDCA4O5uLFi2zZsoXmzZvj7e3Nf/7zH71uSQjbpmlaah/WKTFR0+rW1bTixTXt0aOkp6dN0zQnJ007dky/0MTLJSYkavW9DmrFXC5oj+7GaLdu3dLatWunubu7a+7u7lq7du20O3fuGL0H0BYsWKBpmqY9fPhQCwgI0HLmzKm5uLhoBQoU0Dp27KiFhYWlKY6oqCgN0KKiokx0Z0JYpVT9vLftOYYnQkOhYkUYMQKGDwcgLg5KlYLy5WHVKp3jE8n6tX8wbb/zZ+3ofTQZod9hOTLHIBxEquYY7CMxgNr4NnUqHD8ORYsCsHgxtGsHO3ZAzZo6xyeec/fiXUoXjaVarousvFZN11gkMQgH4WCJ4cEDKFsWSpeGtWvBYCAxEfz81I7o7dvlbGhr81HpbSw9+Rqhex+Qr4q+E8WSGISDsMMieinJkkWtT12/Xq1UQm13mDhRHeGwerXO8QkjW789xLyTdZjUNkT3pCCEMGY/PYYnWrRQpblDQ8HLC02DBg1UUdbDh8HZWe8AxaM7MVTwiSBXxntsu10Opwz6/34iPQbhIBysx/BEYKAaVho0CFDDR5Mnq6mH77/XOTYBwNhmu7n0OA/zFmfRPSkEBgZSpkwZqlTRb+JbCGtjfz0GUDUxPv5Y7W14800AOnWCNWvgzBnw8tI1Ood25PdT+L1XhOH1dzJicz29w0kiPQbhIBxs8vlpT8aPLlyAo0cha1auXoUSJeCjj+Cbb/QO0DElxCXgn+0kDxNcORhZENesrnqHlEQSg3AQDjqUBGr8aN48uH4dhg0DIE8eGDoUpk9XpZWE5U1/bzv7H5Zm3rRHVpUUhBDG7LPH8MTUqeqQhu3boWZNHj1Sq1nLlpVVSpZ2ZtNFKgbkoluFfUw7XFfvcJ4jPQbhIBx4KOmJhASoVQtu34ZDhyBzZpYtg3ffhb/+gmbN9A7QMSTEJVDH+zjXYzw5fMWbLLmy6B3ScyQxCAfhwENJTzg7w8KFEB6edBRoq1YQEAC9e6vFS8L8vvnPdoLvlePHaVFWmRSEEMbsOzEAlCwJX3+tDvbZsAGDQZ3tExEBY8fqHZz9O77qLF+u9Wfg69uo2bOC3uEIIVLBvoeSntA0aNIEjhxRq5Ry5GDcOBg9Wo0wlSund4D2KTY6Fv/c53mU4MrB63nJ6JXx5W/SiQwlCQchQ0lJDAaYPx8ePVKnvmkan30GxYpBjx6QmKh3gPbpywbBHIspyqL5sVadFIQQxhwjMYBarzpnDvz+OyxahJsbzJql6ijNn693cPZn08QDTNlfj4nNdlKpXRm9wxFCpIFjDCU9rX17+PNPVTipUCE6dVJfHj+edGy0SKebJyKpUC6BCtnCWRdRWfeyFykJDAwkMDCQhIQETp8+LUNJwt7JctUXunsXKlUCHx/Yto1b91wpW1aV5169Wkpzp5eWqPG27z723CzMkUOJ5K7oo3dIqSJzDMJByBzDC3l5wa+/wsGDMGQIOXKo4npr18IPP+gdnO2b2XYbq29UZcHwCzaTFIQQxhyvx/DEd99B//6wciW0aEHXrvDbb2rhUuHCegdnm/YvOkmtDwvRvfweph+xvt3NKZEeg3AQMpSUIk1Tu922bIFDh4jOXogKFaBQIdi8WR3yI1Iv8lwUfqXuk9v1DtuuFsPN07ZWIUliEA5ChpJS9GQJq5cXtG2LR8Y4FiyArVvVQXAi9RLiNd73v8DDBDf+2Ohhc0lBCGHMcRMDQLZs/59vGDCA+vWhb18YMgROnNA7ONsxImA3m2+W59cxp8lfs4De4Qgh0smxEwNA1aqqFndgIMybx4QJULAgvP8+xMToHZz1WzXyIF8F+TOh/kbe+LKG3uEIIUzAcecYnvXJJ2p50ubNHPGoRbVqasvD3Ll6B2a9Tm+4QJXG2Xkz9zH+CK+OIYPtHqgtcwzCQcjkc5o8fgwNG0JoKOzfz/y/C9C1K/z0k0oQwlhUeDT+xW6QiBN7z+bAI7+n3iGliyQG4SBk8jlNXFxUuYzMmaFFCzq/94COHVUtpePH9Q7OusTHPOa9Sme4FufNqpXYfFIQQhiTxPC0nDlVfYzTpzF0aE/gtAQKF1YH+0RH6x2cddASNfpW2s7mWxVYNvkcJZsU0TukdAkMDKRMmTJUqVJF71CEsBoylPQif/0FLVtC796c7PEt1aobqFNH7YVztt1hdJOY9MYGvghqxLyuu+n2fXW9wzEZGUoSDkKGkl5Z8+bqYJ9p0yi19huWLlUlM4YO1TswfS3osJkvghox8s2ddpUUhBDGMugdgNXq2VMdCTpoEE2+92LKlK4MGAClS0OnTnoHZ3l/DQ2m+891+LjcTkZukGWpQtgzSQwpGT9eVWPt3p3+i7NyonsbuneH3LmhcWO9g7OcrdOP8N6E12iR9wCBB/0xOEkJWiHsmSSGlBgMakjp/n0M7T9k5h9ZiYh4i3fegaAgtTfO3m2ZdYK3+halplcoi45VxNlFRh+FsHfyr/xlnJxUTaXmzcnQ5h2WdlrPa6/BW2/ByZN6B2deW+ac4q1PClDD8zh/nSohx3MK4SAkMaRGhgywZAk0akTmNs35q9tKfHygfn37TQ5bAo/TtEd+anoc489TpciUy13vkIQQFiKJIbXc3OCPP6BNG7J3e4fN3Rbj7Q316tnfBrh1E0Jo2rswtbyOs+p0aTL5yPJNIRyJJIa0cHFRNTI++ohcn7YjqO0ccudWPYcjR/QOzjTm99xH86HleDPnYVadLStJQQgHJIkhrZycYOZMGDwY7y978E+d0eTLp/HGG7B3r97BvTotUWNMw+10nV2FbsW2sPxCZTLlyKx3WEIIHUhieBUGA0ycCN98Q47AMfyT6wNKFkugXj1VUcPW3L/zmHYl9jLy79qMq/c3s041IEMWN73DsggpiSHE86QkRnqtXQtt2xJToCQf5tvCio1ZGD0ahg2zjeNBT267wTtNHnDpYU7m99zPezPr6R2SLqQkhnAQUhLDIpo2heBgMsXc5vddeRnxznFGjoS334Y7d/QOLmW/jzxGlXqZ0eIes++nkw6bFIQQxiQxmELZsnDwIE5NGjHqj3KsaTSN4GANPz91aqi1uXf7MT1f38d7Y8rxVs597D3lRen2r+sdlhDCSkhiMBVPT1i6FObOpcmWwRzIWo/sLveoVg2+/BIePdI7QGXj7POU843k5wOlCXxrLUuu1CFrkVx6hyWEsCKSGEzJYIDu3eHwYQoV1Nh1OgfDK6xk8mSNihVh+3b9QgsLvc/7pQ7SqGcRimW4xNFVF/hkdVObPo5TCGEekhjMoUQJ2LIF15nfMeJ8Z0Jcq+Ede5k6daBNGzh71nKhXL8Yw+AG+ylV1omgU3mZ/+5a/r5dmcJvl7dcEEIImyKJwVycnFTp7rNnKdOxCtsvF2F+5t7sWh9FqVIaHTrAsWPmu/ylU4/oXecIhQrDrM0l+LTCZk4fjaXz700xuLma78JCCJsnicHccuSAwECcLpyjc/cMnI4rzBTnwQStuEP58lCrFixYAFFR6b9UdDQs/OoqDQufpUgpF5Zuz8Ow19ZyaX8k4w83w6NcgfRfxMzGjx9PjRo1yJw5M15eXql6j6ZpjBo1ijx58pApUybq1avHcXurUyKEBUlisJT8+eHbb8kUdor+n7txzq0sS2lD5kM76dIFcnonEtAwkTFjYNOm1J0xffcu/LMhnom9wmle+iy5vGLpMiw38eHXmFN3CZcO3ubLQ++Qzc92zmWOi4ujdevW9OzZM9XvmTx5Mt988w0zZsxg37595M6dm4YNG3Lv3j0zRiqE/ZINbnqJi4N16+Cvvwj7M4RVN/1Z59yM3QZ/7sSrDVZ5sz8kb87HeHpoeGZNgMePiXmocT3SiWu33LjywAsAd6Kp4nSQpuXDadM7J/na11dF/2zYwoUL6d+/P3fv3k2xnaZp5MmTh/79+zN48GAAYmNj8fHxYdKkSXz88cepup5scBMOIlUb3OSgHr24ukKLFtCiBQUSE+lz6BB9/v6bxP3zOX0sjr3nvTl9uyARt3MThSdReGJAw41YKhlu0DRHHEVLPqZK/ayUfKccTq/XVEX+HMyFCxeIiIggICAg6Tk3Nzfq1q3Lrl27kk0MsbGxxMbGJn0dnZoumhAOQhKDNXByAj8/8PPDCSgFlNI0tXX6/n24dw8SE9VeCQ8PcHcHZ1lmChAREQGAj4+P0fM+Pj5cunQp2fdNmDCB0aNHmzU2IWyVzDFYK4MBsmeHAgXUzury5dX/e3nZXFIYNWoUBoMhxcf+/fvTdQ2DwbiHrGnac889bciQIURFRSU9wsPD03V9IeyJ9BiE2fXu3Zu2bdum2KZQoUKv9Nm5c+cGVM/B19c36fkbN24814t4mpubG242Pg8jhLlIYhBm5+3tjbe3t1k+u3DhwuTOnZtNmzZRqVIlQK1s2rp1K5MmTTLLNYWwd2lZlSSE2RkMhgJAduBt4DOg9r8vndU07f6/bU4CQzRNW/Hv14OBIUBn4AwwFKgHlNQ0LVVrVg0GgwcQBXhqmiYz0cKhSY9BWJsxQMenvj7073/rA1v+/f+SgOdTbSYDmYCZQDZgDxCQ2qTwr3v/fqZsfhAOT3oMQgghjMiqJCGEEEYkMQghhDAiiUEIIYQRSQxCCCGMSGIQQghhRBKDEEIII5IYhBBCGJHEIIQQwogkBiGEEEYkMQghhDAiiUEIIYSR/wHmPKDoFmyKlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "Graphics object consisting of 2 graphics primitives"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sp = plot(sin(x), (x, -pi, pi), color='red')\n",
    "pp = plot(p[0](x)/p[1](x), (x, -pi, pi))\n",
    "(sp+pp).show(figsize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the function to bring the polynomials from scipy\n",
    "into SageMath objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x*(1.0 - 0.116666666666667*x**2) has type <class 'sympy.core.mul.Mul'>\n"
     ]
    }
   ],
   "source": [
    "fp0 = p[0](x)\n",
    "print(fp0, 'has type', type(fp0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.116666666666667*x^2 + 1.00000000000000)*x has type <class 'sage.symbolic.expression.Expression'>\n"
     ]
    }
   ],
   "source": [
    "sfp0 = SR(fp0)\n",
    "print(sfp0, 'has type', type(sfp0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Numerical Solving of Ordinary Differential Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We end with the numerical solving of an ordinary differential equation,\n",
    "taking the pendulum problem: ``diff(theta,t,2) = -diff(theta,1) - 9.8*sin(theta).``\n",
    "We need to define a system of first-order differential equations,\n",
    "introducing an extra variable velocity which is the derivative of theta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[velocity, -velocity - 9.80000000000000*sin(theta)]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t, theta, velocity = var('t, theta, velocity')\n",
    "oderhs = [velocity, -velocity-9.8*sin(theta)]\n",
    "oderhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right hand side of the system has to be defined as a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pend(y, t):\n",
    "    return [y[1], -y[1]-9.8*sin(y[0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate the function symbolically, as a verification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[velocity, -velocity - 9.80000000000000*sin(theta)]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pend([theta, velocity], t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "Then we import odeint of the scipy.integrate package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.31415927,  0.        ],\n",
       "       [ 0.31400806, -0.03015808],\n",
       "       [ 0.3135566 , -0.05998775],\n",
       "       ...,\n",
       "       [ 0.00157   ,  0.00380359],\n",
       "       [ 0.00160711,  0.00361063],\n",
       "       [ 0.00164228,  0.00341607]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.integrate import odeint\n",
    "from numpy import linspace\n",
    "npts = 1000 # number of points\n",
    "tspan = linspace(0, 10, npts)\n",
    "sol = odeint(pend, [pi/10, 0], tspan)\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAApiElEQVR4nO3dd3xUZfY/8M+TQmhJXAuG0GGRYkWCSrEhoK8AuuKiC3b5iQUWWOsiriKKARVFTRYX8Mta0MV1AQuioLIIS2+CFGkRolKkmNAygeT8/jhJDClk7sy9c+/MfN6v17w0M/c+9wyBOXOfch4jIiAiIgKAGLcDICIi72BSICKiUkwKRERUikmBiIhKMSkQEVEpJgUiIirFpEBERKUiLikYlWSMMW7HQkQUbuLcDiBIFVbe5ebmIjk5Gbm5uW7EQ0TkNZa+IEfcnQIREQUuYpJCVlYW2rZtiw4dOrgdChFR2DJhXvuoQvB5eXml3UdJSUluxERE5CXR3X20a5f+d9gwV8MgIgpLEZcU6tbV/06ZAnz6qbuxEBGFm4jtPmrdOhfx8UlYu9aNsIiIPMNb3UfGmAeNMdnGmHxjzEpjzOWnOLaLMeZ/xpj9xphjxphNxpi/BHLdZ58Ftm8Htm0LPHYiomjjaFIwxtwCYDyA0QDaAVgAYLYxpnEVpxwBkAngCgBtADwH4DljzECr1+7RA2jRApg5M4DAiYiilKPdR8aYpQBWicgDZZ7bCGCmiAz3s43pAI6IyO2VvHzK2UcvvZSE778Hpk0L9B0QEYU9b3QfGWNqAGgPYE65l+YA6ORnG+2Kj51f1TE+nw95eXknPUpccw2wbBlQVGQ5fCKiqORk99GZAGIB7Cn3/B4AKac60RjzozHGB2AFgCwRmVzVsRkZGUhOTi59NGrUqPS1zp2B3buBr74K+D0QEUWVUExJLd/FYyp5rrzLAaQBuB/AMGNMv6oOHD58OHJzc0sfOTk5pa/FxQEDBgCrVwcYORFRlHGyIN4+AIWoeFdQDxXvHk4iItnF/7vOGHM2gJEA3q/s2ISEBCQkJFTZ1hVXADt2+BkxEVGUc+xOQUQKAKwE0L3cS90BLLLQlAFQ9ad+NWrWBN58M9CziYiii9Ols18G8I4xZgWAxQAGAmgM4A0AMMZkAGggIncU/zwIwE4Am4rP7wLgEQCvBxpA9+5ASgqweTNwzjkBvw8ioqjgaFIQkWnGmDMAPAWgPoDvAKSLSEmHTn1okigRAyADQDMAJwBsA/BXAP8INIYaNbQLqaT8BRERVS1iy1yUrZKang6cfz4wdmzI4yMicpuldQrhvvOaX8aOBX7+2e0oiIi8L+KqpFZm5Urg4YfdjoKIyPsi5k4hKysLWVlZKCwsrPDaHXfoqub8fJ2NRERElYuKMYX9+4G0NGDuXOD3vw95jEREbvJG7SMvOeMMYNIk1kAiIqpOVCQFAJgwQe8UiIioahEzplCdN95guQsioupEzZ3CtGnAoEFuR0FE5G1RMdBcYvduLXlBRBRFONBcmexsIDUV+OEHtyMhIvKuqEkKzZoB334LNG3qdiRERN4VMUkhKysLbdu2RYcOHao8ZuxY4MMPQxgUEVGYiaoxhbffBi6+GDjvvJDFR0TkNo4pVCUtDSizWycREZUTVUlh/ny9WyAiospFVfcREVEUYvdRVTZt0jpI+/a5HQkRkTdFVVJo0UIL451xhtuREBF5U1Qlhfh4IDaWu7AREVUlqpICALz4IrB4sdtREBF5U8RUST3VzmtlLVwYooCIiMJQ1M0+evNNYPt2YPTokMRHROQ2S7OPIuZOwV8tWgC1a7sdBRGRN0VdUrjqKuDIEbejICLypqgbaF66VKekMjEQEVUUdUkhLU3LXdSp43YkRETeE3VJITZWH0ePuh0JEZH3RF1SAIAbbgBmzXI7CiIi74m6gWYA2LoVqFXL7SiIiLwnKu8UNm4EFi1yOwoiIu+JmKTgz3acJWbOBN5/3/mYiIjCTdStaCYiijLcT6E6e/cCr77qdhRERN4TlUkhJweYNg2opnYeEVHUYfcREVFk81b3kTHmQWNMtjEm3xiz0hhz+SmO7WOMmWuM+cUYk2eMWWyMudaJuMaPB3btcqJlIqLw5WhSMMbcAmA8gNEA2gFYAGC2MaZxFadcAWAugHQA7QHMA/CJMaadnXGJADNmaDcSERH9xtHuI2PMUgCrROSBMs9tBDBTRIb72cZ6ANNEZFQlL7P7iIjo1Lyxn4Ixpgb02/6Yci/NAdDJzzZiACQCOFDVMT6fDz6fr/TnvLw8v+L75hsgLg7o5FckRETRwcnuozMBxALYU+75PQBS/GzjYQB1AHxQ1QEZGRlITk4ufTRq1Mivhj/7DPjiCz+jICKKEo51HxljUgH8BKCTiCwu8/wIALeLSOtqzu8HYDKAG0TkyyoOk8ruFBo1asTuIyIi5Y3uIwD7ABSi4l1BPVS8ezhJ8QD1mwD6niIhAAASEhKQkJBgObi8PGDJEqBHD8unEhFFLMe6j0SkAMBKAN3LvdQdQJXl6IrvEP4JoL+IOFbgeuFCYMAAp1onIgpPTpfOfhnAO8aYFQAWAxgIoDGANwDAGJMBoIGI3FH8cz8AbwMYCmCJMabkLuOYiOTaGVh6OqekEhGV5+g6BRGZBmAYgKcArIGuQ0gXkR3Fh9SHJokS90ETVRaAXWUejlQq2rJF1ywQEZGK2jIX+flAcjIwZw5w5ZWOxkhE5CbPDDR7Ws2awPbtQIMGbkdCROQdUVkltQQTAhHRySImKVjZea1E//7AE084GBQRUZiJ2jEFAFi9GqhbF2jZ0rH4iIjcxjEFf7WztfYqEVH4i5juo0AsWgQMHOh2FERE3hHVSaFuXQ42ExGVFdVjCkREUcBb23F6mQgwZQpw9KjbkRAReUNUJ4W8POC554AdO6o/logoGkT17KPkZGDbNrejICLyjoi5Uwhk8RqgXUcnTjgUFBFRmIn6gebu3YGOHYFRoxyJj4jIbVy8ZsU//qHdSERExKSA5s3djoCIyDsiZkwhUF9/DfzlL25HQUTkDVGfFOrUAc46y+0oiIi8IeoHmomIIhxXNFv15ZdAUZHbURARuS/qk8LPPwO9e+vWnERE0S5iZh9lZWUhKysLhYWFls5LTdVyF/HxDgVGRBRGOKZARBTZOKZg1S23AP/8p9tREBG5L2K6j4LRrx/QqpXbURARuY9JAcAf/uB2BERE3sDuI2j57Dlz3I6CiMh9TAoAvvoKyMx0OwoiIvdx9hERUWTj7CMiIgpMxCSFQHdeA4DsbKBBA2DPHgcCIyIKI+w+AnD8OPDBB8Cf/gTExtoeIxGRm9h9ZFV8PHDrrZGREF5/HTj9dCAxEbj/frejIaJww6RQbMYM7UYKZwMHAkOGALm5QEGBbjV64YVuR0VE4YRJodgbbwBLl7odReA++giYNAno2BE4cQLw+TRBrF0L3Hab29ERUbhwfEzBGPMggEcB1AewHsAwEVlQxbH1AYwD0B5ASwCviciwUzTPKanFUlOBmjV1IZ4p04PYrx8wbRqwfDnQvr178RGRa7wzpmCMuQXAeACjAbQDsADAbGNM4ypOSQDwS/Hx3zoZWyQZNw7YvVvvFky5X//UqZowBgxwJzYiCi9Odx89BOBNEZksIhuLv/XnAHigsoNF5AcRGSoibwPIdTi2k4wbp7OPwtHEiUCvXsD551d8LSYGeO01YNMm3WGOiOhUHCuIZ4ypAe0GGlPupTkAOtl1HZ/PB5/PV/pzXl5eQO306gWkpdkVVeh8/jmwdSvw4YdVH9OnD9C8OTBqFNCtW+hiI6Lw4+SdwpkAYgGUXxK2B0CKXRfJyMhAcnJy6aNRo0YBtdOqFXDllXZFFTqTJgHt2lV+l1DWs88C69dz21EiOrVQzD4qPxhsKnkuYMOHD0dubm7pIycnJ6B2jh0D3nkHCKe1fIcOAatWAUOHVn/sDTcAycnaTUZEVBUnk8I+AIWoeFdQDxXvHgKWkJCApKSkkx6ByM4GnnxS92sOF3PnAnFxQO/e1R8bFwf07Qt8/LFOWSUiqoxjSUFECgCsBNC93EvdASxy6rqBatsW2LFDv02Hi88+Ay65BDjtNP+OHzRIF7ZxwJmIquJ099HLAP6fMeYeY0wbY8wrABoDeAMAjDEZxpi3y55gjLnIGHMRgLoAzir+ua3DcYadY8eA776ztmtc48ZAp066boGIqDKOJgURmQZgGICnAKwBcAWAdBHZUXxIfWiSKGt18aM9gP7F//+Zk3GWuOoqYNasUFwpeOvXa1fX5ZdbO69fP2DFCuDwYWfiIqLw5vhAs4j8XUSaikiCiLQXkW/KvHaXiFxV7nhTyaOp03ECWjvoggtCcaXgff01kJKiDyvS04HCQmDhQmfiIqLwxtpHZfTvDwQ4ozXkZs0Crr7a+nlnnQWce64OOBMRlcekUMbWrcC3YVBcY98+YNEi4PrrAzu/d2/tQiqz5o+ICEAEJYVgdl4r8eabwNixNgblkBUrdMC4ugVrVenWTUtrr1tnb1xEFP6481oYyswEdu4EXngh8DZ699Yy2088YV9cRORJ3qmSSs5YuxZo2jS4Nq66Cli82I5oiCiSMCmUsW6dTvEsKHA7kqodPQpMmRL8jmpXX62VU/fvtycuIooMTAplpKTo4G18vNuRVG31av1A79w5uHYuvhg4ckRLZRARlWBSKOOss4BHH624UY2X5OZqVVQ79OgBfPKJPW0RUWRgUihnyRLtovGq//wHOH7cnrZ69dL6SUREJZgUyund27sF44qKgJUrgZtvtqe9664DkpKAZcvsaY+Iwp9jO6+Fqx9/BBIS3I6ickeP6kZALVva017durreYflyrbZKRBQxdwp2LF4DvJsQAP3wnjYNOP10+9q8+GKd4kpEBERQUhg0aBA2bNiA5cuXB9XO+PHAiy/aE5PdVq0C/vxnewfC09J0Ki433iEigN1HFfz+9979gDzjDP831PFX+/bAmjWacNiFREQRc6dgl169rG1cE0oPP2x/wkpK0gV7H31kb7tEFJ6YFMopKNAZPl6zZw9w993OJKwbb9Q7BSIiJoVyvvwS6NoV8FqdwAMHtKx3UZH9bXfurO/34EH72yai8MKkUE56OvDTT95b1TxrFtC6NVCjhv1tN2mi012zs+1vm4jCC5NCJerWdTuCis47TxOWE5KSgLg44IsvnGnfDV6dLEDkdUwK5RQW6mCz1zaguf9+YPdu59rv2ROYMcO59kPhueeAxES9y4uP1//WqqV/dk50uxFFoohJCnYtXouJ0Wmav/udTYHZID8feOgh3UPaKb16aU2lX35x7hpO2b0bqF8f+NvfdIvRCy4AbrhBCweKAP/4B1CnDus8EfmDO6+FgY0bge7dge++s3+dQomjR4F+/YARI8JrvcL33wMXXaSJc/Bg4PXXKx4zcaK+dvw48Ne/AhkZIQ+TyE3ceS1Y27YBOTluR/Gbn3/WD2unEgIA1K6t36a9OB23Kj//rGU6CguBBQsqTwgAMHCglhxv1gwYMwa4997QxkkUTpgUKvG3vwW3/7HdDh50djyhROvWwL//7fx17FBYCFx6qQ4oL10KdOly6uNr1QK2b9dzJk/WfTOIqCImhUq8/bbWQPKKWbO0RpHTbrxRi+4dOOD8tYJ1yy06dXjqVGubDi1Zon+WL72kYw1EdDImhUrExQGxsW5H8ZtLL9W1BE5r00ZLXmzY4Py1gvHFF7rZ0N13A3/8o/Xzly4FmjcHBg0CVqywPz6icMakUIl584Brr3U7CvXDD8Czz+pWoU6Li9Pks36989cKlAhw110a5+TJgbURE6NlPerU0QH8w4dtDZEorDEpVKJ5c+8UxUtNBUaPBurVC831mjUD5s8PzbUC8dhjWgdq+vTgVp0nJ2tJk7w8TQxEpJgUKtGkCfDAA25HoTIzte87VN1ZnToBH34IHDsWmutZceCAzjC66SaddRSsDh30LmzJEk28RBRBScGuxWslpk3Tuftu69ULyMoK3fW6dAFatQIWLgzdNf11331AzZrApEn2tfnEE0DHjsDIkd4fSyEKBS5eq8SRI/rB+OmnujDKTR07Ar1764dXqAwdCrRoAQwZErprVmfnTh0IHzTI/unCeXl6d5iUpEUBY8L4q9KxY8CyZcCOHfo+WrTQuyovbzNLjrPU0cqd1ypRpw7w449uR6Gefz70K4zPOcd7+zYPHaqFCkeNsr/tpCSdhnz99brQLdABbDeI6Paxr70G7NpVdY0nY3Syws036wK+OnVCGyeFD94peNiuXUCDBtqt0bp16K47f76uA/j5Z298a965Ezj/fOCRR3RhoVNuuUUHsOfNq34xnNvy83VsZfZsTQzGAGeeqWswrrhC7xAALQOycKFuubp372/7hNSrB4wdqzO5KOJZm5IhIuH8qCA3N1cASG5ubmUv++2mm0SeeiqoJoJ26JDIggWhv25hoUjt2iKzZ4f+2pW57TaRM88UOXrU2esUFIikpIicdZaIz+fstYIxZIiIMSKAyOmni2Rl+X/u/PkiF1742/kJCSKPPupYqOQNlj5XPfA90JuGDwfuvNPdGN59VwdXQy0mRovjLV4c+muXd/Ag8Pnn+udQq5az14qP1zuFX3/VuwavWb8eOOMM7SpKTgbmzAH27wcefND/Nq64Qu8aTpzQPb+N0e6nGjWAxx93LHQKI0wKVWjfXtcruOmaa3Sw2w2dOwNbtrhz7bJGj9bqpo89FprrdewI3HMPMHOmt/aXGDFCu9AOHtRKrwcPBre+IiZGS30cOwa88oouXHzhBZ3d5aW6X1U5fhx4/32gTx/9d5qYqIktJkYTXflHTIy+npQENG0KdOsGPPOMViCmk3FMoQpr1mh/7fDhwYQXnGuv1bn0zz0X+muvXq17OHzyiXs70fl8OujdvXtoB39F9INm/36dcODm0FRRkSaqZct0HGDFCqBRI2euNWaMTs31+XQg+tVXgQEDnLmWVYcPa3xTp+rvpPzOevHxWum3bl1NEDVr6tqeEyd0avmhQzqr8Nixynfli48HTj9dZxv26wf07avteVVhoVZy3rpVZ8ydcw5w5ZVVHm5pTCFsk4IxxuTm5hb5fD74fL7S5w8dOoS2bdsiJycnqKSwcKEu4nKzMN7ixXrH4sS+zNU5dkyvPXGie4Ou//d/wFNP6Qdiampor/3DD1por1074OuvQ3vtEvv365eC/fv1C8IHH4TmuiNGABMm6AdPYiIwbpw73WkrV/72+y8o+O35xES9a/rjH3U2VWKi9bYLC7X444wZwKJFWkG3snInycm6yv+aa/R6bdsG/n78tXOnxrR0qXYZ/vyz3hkeO6ZxVyY+Hti3r/LXkpOTkwEcEj8/7MM5KSQByHU7DiKiMJAsInn+HBjOSaHSO4Vdu3bhkksuwYYNG9CgQYOA2u7QoQOWL1+OFSv027LVGjsl5weqQ4cOWLBgOTZutFYWGtDus0aNGgV1p1QS/1tvaVkJq5VEg3n/Zf/sb7pJq6H6Wzbczvdeols3LZ73xRf6rd3q+YFc/5FHluO++7Sf/9NPgcsu8+/cYN9/ZbEXFekakalT9VtqTIwOVk+YUPHuLZD3fuKErtifNAnIyekAQM9PStJ9w59/Xrt1quPE776EiN61TJ+u3+B37NDuqOPHK7RQGr8VcXFAYWEHNGiwHI0aAeeeq7/zTp10Snp1qnvvVu8UwnbxWnVvMDExMeC/HLGxsfD5knDttTq2cO651s8PpusqNjYWc+YkYfjwwHeAS0pKCur9JyUloXdv7cvfv19voa2eH8y1X3xRVxl37Wq9DTvee4l587QPv08f/TCobve7YH/3Bw7EYuDAJCQn6/arDRtabyPQ919V7G+9pY+MDO3X/+9/dXV5rVr6wfXYY0CPHv6992PHdKHglCnAunUnl5KJjY3FY48lYeTIwLtM7fzdl9W1a9V/Fw8c0N/Vn/4Ui/vuS8KBA7+NXcTE6J/TaacBKSn6d/r3v9cxq7LvsW3bWGzYENzgVVXv3d87hLInhPOjgpycHAEgOTk5fk/iLS8zM1NERPLygjs/mOvv2SOSnW39XDvWaZTEX1Agkp4uMmNGYOcHeu1ffhGpWdPa/HsRe997WatWicTFibRsKVJUZP18f918swiQKY0bixw5Yv38YN+/v7Fv3izStatIjRq61qHkEROTKUlJIo0aibRqJdK6tUjTprrGpFatk48tWSPRrp3Ie+9Zu35lnPrdh+p8h9+7pc9Vtz/UPZkUvODGG0X+/Gfr59m1eK/EkCEiDz9sS1N+e/ppkYYNNSlZYfd7L+v11/Vfy9VX2960FBaKXHyxtn/ppdUnnqo4+f6r8uuvIiNH6oK45GRNnuU//GNjNSmkpopcd53Iu++KHD9ubxxuvHevYFKoJins3btXAMjevXur+7OsVu/eIuPHB91MQLZs0RXNVuXn58vTTz8t+fn5tsTx/vsil12mH1yhUFQkcvbZgSVEu997efffr/9i+va1r829e3UFNaArt4Ph9Pv3Mr73U753S5+rYTvQXMzR2kfz52t/cqgXsfl8OhXu88+Bq64K7bXL+/FHoGVLHWy2OrYSiPffB269VacINm3q/PWs6tVL98zu2zf4KaKzZulYRUGBLiR7+GF7YiQqJzrWKRSL2IJ4a9YAF1zgjYJ07dvrh9eIEc5fq1s37XT46ivnrxWobt00vo4dgf/9L7Ad4G6+Gfj3v3WR1ddfa1tEDrH0N9QDHzneNW8e8PTTob/up5/qKl4vJAQASE/Xb7VO27RJ70geesj5awXjyy/1Q33xYp0u+d13/p87Z47eBf7737pnx759TAjkLR752Ame3TuvAb/VTAm1rl2Bjz4K/XWrcsMNeueyZ4+z13nxRf3AvO46Z69jh2nTtMsnN1dX13bsCPzyS9XHf/aZzuu/9lpdOTt6tCZB7mtAXsPuIw8aMkTrtrzyituRqMJC7coaMQLo39+Zaxw+rHtGDBighcrCxU8/6bjP1q36c+3aupdBw4b6nrZt02RaUp6ga1etJ+XlujoUcdh9ZJfjx7X2T9m6K6Fwzz3e2vwkNlbrH338sXPXmDpVV88OHuzcNZzQoIFWk12/XnfIO35cF2XNng0sWKAbJSUlaXnr/Hwdi2BCIC9jUjiFX37RLo1TdQs44cEHtaKjvzIyMtChQwckJiaiXr16+MMf/oDvv//e1pjS07UwWWVFw+wwfjxw6aW6ZWQwMjIyYIzBsGHD7AjLb23bagGzgoKTZ+kXFemK16ws5/ZJ/umnn3DbbbfhjDPOQO3atXHRRRdh5cqVzlzMQ06cOIEnn3wSzZo1Q61atdC8eXOMGjUKRVXtSRrmvvnmG/Tu3RupqakwxmDmzJknvS4iGDlyJFJTU1GrVi1cddVVWL9+veXrMCmcQmqqfgsMsIRSwG6/3do02Pnz52PQoEFYsmQJ5s6dixMnTqBHjx44YiWzVOPqq7VGy4IFtjVZauFCnYL6yCPBtbN8+XJMnDgRF1xwgT2BhYGDBw+ic+fOiI+Px+zZs7FhwwaMGzcOp1VXjyMCjB07Fm+88QYyMzOxceNGvPDCC3jxxRfx+uuvux2aI44cOYILL7wQmZmZlb7+wgsv4OWXX0ZmZiaWL1+OlJQUdO/eHcYYa3VkrS5s8NijArtXNm7apKs2Q2XDBl3NfOJE4G2ULOCbP3++fYGJblF61122NikiIr16iVxySXBtHDp0SFq2bClz586VK6+8UoYOHWpLbF73+OOPS5cuXdwOwxU9e/aUe+6556Tn+vTpI7cFuwowDACQGWXqzxQVFUlKSoqMGTOm9Ln8/HxJTk4WAPeJhc9V3ilUo39/Ld4VKnXr6oYZwcx6ys3ViuKn+1Ne0oKePbVaqJ1jLD/+qDNz7rknuHYGDRqEnj17olu3bvYEFiY+/vhjpKWloW/fvqhXrx7atWuHSZMmuR1WSHTp0gVfffUVNm/eDAD49ttvsXDhQqSnp7scWehlZ2dj9+7d6NGjR+lzCQkJuFJ33ulkpa2wrZIaKt984/zewGUtX66zcAJZEAXond9DDz2ELl264LzzzrM1tptuAoYN05XW119vT5svvKD7Dg8cGHgb//rXv7Bq1aqgSlaHq+3bt2PChAl46KGH8MQTT2DZsmUYMmQIEhIScMcdd7gdnqMef/xx5ObmonXr1oiNjUVhYSFGjx6Nfv36uR1ayO3evRsAcPbZZ5/0fPHPKVbaYlKoRo0aOrgaqtmteXm6cXygBg8ejLVr12LhwoW2xVQiKUkHgydPticp+Hy6u9rQoYEnwZycHAwdOhRz5sxBzZo1gw8qzBQVFSEtLQ3PP/88AKBdu3ZYv349JkyYEPFJYdq0aXj33Xfx3nvv4dxzz8WaNWswbNgwpKam4s4773Q7PFeYcv+QtKep4tT9U7LS1+TBR6nMzExp06aNnHPOObaOKYwaJXL55bY05ZcJEwIrmS0iMnjwYGnYsKFs377d1pjKmjJFK14GUtq5vGee0fLLwYzZzJgxQwBIbGxs6QOAGGMkNjZWTgQzOBMGGjduLAMGDDjpub///e+SmprqUkSh07Bhwwolp5999llp1aqVSxGFDsqNKWzbtk0AyKpVq0467vrrrxcAb0k0jikMGjQIGzZssL0LYcgQ4J13bG3ylGbO1Jk4VogIBg8ejOnTp+Prr79GMys74ljUv7/OQpo8Obh2CguBzExtLzk58HauueYarFu3DmvWrCl9pKWl4dZbb8WaNWsQGxsbXKAe17lz5wrTjzdv3owmTZq4FFHoHD16FDHlBt9iY2MjdkrqqTRr1gwpKSmYO3du6XMFBQWYP38+ACyy1JiVDOLBRwVO1FUP1ZfNo0dFPvjA+nkPPPCAJCcny3//+1/ZtWtX6ePo0aP2Bykid98d/GyhceNEYmICvys6lWiafbRs2TKJi4uT0aNHy5YtW2Tq1KlSu3Zteffdd90OzXF33nmnNGjQQD799FPJzs6W6dOny5lnnimPPfaY26E54tChQ7J69WpZvXq1AJCXX35ZVq9eLTt27BARkTFjxkhycrJMnz5d1q1bJ/369ZP69esLgESx8Lnq9oe655PC+vXaXWLD9gzVWrNGpH596/soQPsMKzymTJniSJzLlunOaN98E9j5BQW6G9eNN9obV4loSgoiIp988omcd955kpCQIK1bt5aJEye6HVJI5OXlydChQ6Vx48ZSs2ZNad68uYwYMUJ8Pp/boTli3rx5lf47v/POO0VEp6U+/fTTkpKSIgkJCXLFFVfIunXrRCx+rrL2UTUKC7U0QZmZXo7Zu1c3BG/RwvlrBatNG91v9vPPrZ87cqRuyL5unVYKJSJHsfaRnWJjQ5MQAC31MGRIaK4VrEcfBZYsAaxW08jL0/fZqxcTApEXMSn44cILtWCb0/78Zy3JHA5uv13Xb1jdLWzgQODECSBCKxEQhT0mBT+MHx+au4UrrwTee8/569ghPl5rFc2Zo4Xy/LFokW4uM2BA6OtJEZF/OKbgIVu3ah3+cFmDVVioe1gnJuqGMadagFZQoMfGxAA5OTqtlYhCIjrHFJzYea3EU0/pHsVO2rRJN7LJy3P2OnaKjQUmTNBkVt1YSPfuwP79utk9EwKRd/FOwQ/Z2Vp6ol07W5qrVGGhbuDevbtz13BKr166h/PEicC991Z8vW9f4MMPde/lceNCHx9RlLN0p8Ck4KcjR5zdT/fvf9ek8OGHzl3DKYWFQNOmWvH03ns1OQB6Z9Cli94Fpadr4iCikGNSsDspLF2qg8D79mlpayd8/z2wY0fopr/azefTtQvZ2fpzfLxuTQloKYtQzN4iokpF55iCk9LSgHnznEsIADB2rC5cC1cJCVqz6fXXdbC8Vi3g4ot1gRoTAlH44J2Cn7Kz9cMuPt62Jk8yZYp+iF54oTPtE1HU4p2CEzp2BP7zH2faPnhQ+9sbNnSmfSIif3FyoJ/WrAHq1XOm7bg44OyznbsLISLyF+8U/LR/P7B6tTNtr12rdwlhttaOiCJQxCQFJxevAbrRzhtvONI0jh3TLiQiIrdxoNkDpkwBzjkH6NzZ7UiIKAJxoNkJ27drFVMncuh33wE7d9rfLhGRVUwKfjp+XEtd2J0Ujh8HzjoL6NnT3naJiALBpOCnVq10XCHG5j+xo0eBjz4CDhywt10iokAwKVjQpInWJ7LT/v3AAw9o7SAiIrc5mhSMMb8zxrxjjMktfrxjjDmtmnP6GGO+MMbsM8aIMeYiJ2O04rXXgPbt7W1z/frw2W2NiCKf03cK7wG4CMB1xY+LALxTzTl1APwPwF+dDCwQLVsCu3fb22Z8PPDMM/a2SUQUKMdWNBtj2kATwWUisrT4uXsBLDbGtBKRSrd8F5F3io9t6s91fD4ffD5f6c95Du5S88orQI0aQFaWfW1+8okmm7Q0+9okIgqUk2UuOgLILUkIACAiS4wxuQA6Aag0KViVkZGBZ0L0VXvSJHvbE9EB7DvusLddIqJAOdl9lAJgbyXP7y1+zRbDhw9Hbm5u6SMnJ8eupitYskQ3i7FLfj7w1lvATz/Z1yYRUTAs3ykYY0YCeLqaw0pqTVQ2q99U8XxAEhISkJCQYFdzp1S/PtC1q33t7dkD3HUXcO659rVJRBSMQO4UMgG0qebxHYDdAM6u5PyzAOwJJFi3NWkCXHSR7jJmh507gc8+s6ctIiI7WL5TEJF9APZVd5wxZjGAZGPMJSKyrPi5SwEkA1hk9bpe4PMBN90EzJ4NdOoUfHu7dwMPPhh8O0REdnFsoFlENhpjPgcwyRhzX/HTEwF8WnbmkTFmE4DhIjKj+OfTATQGkFp8SCtjDADsFhGbJ4Rak5Cg1UztWtW8Y4ezW3wSEVnl9DqFWwGsAzCn+LEWwO3ljmkFvXsocT2A1QBmFf/8r+Kf73c0Uj/ddx8walTw7YgAmzYB11wTfFtERHZxdOc1ETkA4LZqjjHlfv4ngH86F1Vw7rrLvs1watSwv5YSEVEw+JFkUXIy8PHHwbezciWwZQvQokXwbRER2YVJwaJffwW2bQu+nfr1gWuvBYyl7S+IiJzlaPdRKGVlZSErKwuFhYWOXqdLF00MBQXa/ROol17SbTiJiLyE23Fa9Ouvul5hwQLgggsCb+fLL4Ezz9R1D0REDuJ2nE467TRg69bgxwJuuUXbISLyEiaFAPTrBzz/fODnHz4MjBwJ3HijbSEREdmC3UcB2Ftc5q9evcDOX7EC6NMH2LwZqFnTvriIiCrB7iOnTZ8e3KKz7duBRx9lQiAi74mY2Ueh1L8/0KFD9cdVJSaGi9aIyJv40RSAw4eBjh118VkgnnxSZx4REXkNk0IAUlOBjz7SqalWFRQAt98OdO9uf1xERMGKmKSQlZWFtm3bokMw/ToWPPMM8PLL1s/bvFkXrnElMxF5EWcfBejbb4Hjx4G0NGvnvfcesHYtMGaMM3EREZXD2UehMGcOMGSI9fOMAZo2tT0cIiJbMCkE6IEHgF69dIzAiuefBxwuz0REFDAmhQD9+quuV8jL8/+cvXuBRo2AP/3JsbCIiILCpBCghg2BVq2ARRZ2mz58GPjd77hojYi8i0khCJddBtSp4//xM2YA2dnWziEiCiUmhSDs3QtkZPh//MKFukaBiMirmBSCcPfd2h3kz2Y5x44BcXHApZc6HxcRUaAiJimEevEaANSurf/1Jyns3g389JNuw0lE5FURkxQGDRqEDRs2YPny5SG7ZkqKlsGeNav6Y2fN0oFmJgUi8rKISQpuue024Jtvqj/u4491XQMRkZcxKQSpUyftGjrVgrTcXGD1aqBnz9DFRUQUCCaFILVpA+zfD/z4Y9XHLFoEnH020L596OIiIgoEk0KQmjYF8vOBL76o+pgZM4Dzz+eiNSLyPiYFG3Tpoh/8lSkoAD7/nPsnEFF4YFKwwR//CKxcqWML5X34odZJ6tMn5GEREVnGpGCDzp2BpCRg4sSKr736qu7nfNppIQ+LiMiyiEkKbixeKxEbC6SnA6+9phvvlFi4EFi/Hnj00ZCHREQUEO68ZpOcHB10HjBA7xiKioAWLYCYGGDrVm6/SUSusfTpE+dUFNGmZJ+EyZM1QWzaBPzwg85KYkIgonDBOwUbFRUBrVsDW7bozw8/DLz0UkhDICIqj3cKbomJATZvBnbsAE4/HUhMdDsiIiJrmBQc0KSJ2xEQEQUmYmYfERFR8JgUiIioVLgPNFdgjEkCkAsgWUTy3I6HiCicRGJSMAASARySSHtzREQOi7ikQEREgeOYAhERlWJSICKiUkwKRERUikmBiIhKMSkQEVEpJgUiIirFpEBERKX+P4SCvgqoBv32AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "points = [(tspan[k], sol[k][0]) for k in range(npts)]\n",
    "list_plot(points,size=1).show(figsize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A phase portrait shows the displacements and the velocities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4LklEQVR4nO3deVxU9fc/8NcIAqIyLgimEmImirsoi0vuqGkalUsoWm5Z+jWX+qmZhdonzT6VmqiZS4umflxy+egHpdyVJRHcQFASQQRR1AFEWe/vjxPLsA5w5947M+f5eMyD4TJ37mHEOfPezlslCAIYY4yxArXkDoAxxpiycGJgjDGmhRMDY4wxLZwYGGOMaeHEwBhjTAsnBsYYY1o4MTDGGNPCiYGJTkVsVCqVSu5YGGNVZ16Fx/JKOKYTjUYDtVoNjUYjdyiMMW06fVjjFgNjjDEtnBgYY4xp4cRgJIqXvMrLky8Oxpjh48SgME+fAjk5dP/WLeDgQbr/8CHw1lvAkyf0/dChwK5ddP+LL4ARI+j+5cuAjQ3w+DF936ULsHcv3f/+e2DRIrp/9y6wYAGQn0/fnz8PPH+uz9+MMWYoODHI5N49+pqbC0ybBty5Q997ewPLltH9oCBg61a6b24ONG1KXwFgxgzA3Z3uT5wILF9O99u3Bw4cABo2pO9XrgT69qX7nTsDHh50PyMDiImhxCAIwPDhwJkz9LP33gOWLqX7t24VJSDGmGlQVaHsNs9KqoGQEHoj9vUF/vyT3ojT0wGVCpg+HViyBHByAuLi6BN/o0bSxpeVBVha0v0LF4C6dSmR7NsH/PQTcPgwJYlXXqFWSZMmQGAgJZr69bWfKy0trXBWko2NjbS/CGOsIjrNSuLEILKMDHpTVamADz+kN/s5c4Dt2yk5fP89fUqPjQVeflnuaKsmO5sSxdixQK1a1IL5+WdgyBBg1SqgXTvgtdc4MTCmYJwYpHLtGuDiQm+WrVoBfn7UvXP0KNC4cVGXjzH797+B2Fh/nD7tj/v38/DoUQw0Gg2srW2Qnl7UtcUYkxUnBn06dw7o3RvIzAQaNABOngR69QKuXqXkULeu3BHKJygoDT17Uoth2zYbbNsGRERQd9XNm0CHDnJHyJjJ4sQgpvx84NIloHt3Ggdo04be5BwdacaQra3cESpH8a4kMzMbpKRQl9revcAHHwApKdQtFRpKyZUxJhlODGJITqa+9PBwwM2N3tQaNqSWgrW13NEpU0VjDE+eUAvr9GmaYvvoEY3HhIRQi4sxpldcEqOmgoMBBwcaUO7alVoGBX3lnBSqp0ED+tq3L72etWsDFy8CgwdTss3PB6KiZA2RMZPHiaGE69eBFi0AjYYGja9dA+rVo5+p1fLGVl3r16+Hk5MTrKys4OrqirNnz+p03vnz52Fubo4uXbroJa6C6bEeHkBqKiXbs2dpUV7BYruMDL1cmjFWAU4MoO6NKVMoGbz0Ek29tLGhLg5nZ7mjq5ndu3djzpw5WLx4McLDw9GnTx8MGzYM8fHxFZ6n0WgwceJEDBw4UJI469Shr337AvHxgJUVcOQIJWnGmLRMOjGcO0flJwSBZswIAr0h+fhQUjAG3377LaZMmYKpU6eiXbt2WL16NRwcHLBhw4YKz3vvvffg4+MDT09PiSItYm9PXwcNAk6coPt79wLDhkkeCmMmyWQTw/PntDArJITGDbZvL+r/NhbZ2dkICwuDl5eX1nEvLy9cuHCh3PO2bduG2NhYfP755zpdJysrC2lpaVo3MVhaAt260f1u3YB336X7hw4BmzeLcgnGWBlMLjGMGUOrd62saEaMMU+XfPjwIfLy8mBf8BH8H/b29khOTi7znJs3b2LhwoXYsWMHzM1128dpxYoVUKvVhTcHB4cax15Sq1b0bwcAaWlFxQTDw4GkJNEvx5hJM4nEEBcHREfT/VGjqNAcUDT4aexK7rApCEKpYwCQl5cHHx8fLF26FG3atNH5+RctWgSNRlN4S0hIqHHMFZkwAfjoI7o/fz7wyy96vRxjJqcqW3sarNWraf3Bb78B48fLHY10bG1tYWZmVqp1kJKSUqoVAQDp6em4ePEiwsPDMWvWLABAfn4+BEGAubk5jh8/jgEDBpQ6z9LSEpYyZdnAQMDMjMqVt25NA9YF3U+Mseox2hbDjRtUr0gQgO++ozEEU2NhYQFXV1cEBgZqHQ8MDETPnj1LPd7GxgZXr15FRERE4W3GjBlwdnZGREQE3BVY9MnMjL7WrUsFCtu1A27fpn9zxlj1GF2LQRBoRlHt2kWDySqV8cwyqqp58+bB19cX3bt3h6enJzZt2oT4+HjMmDEDAHUDJSYm4pdffkGtWrXQoUQhIzs7O1hZWZU6rkRvvUVfExOpdtXcuTTJwNy8aB8LxljljO6/S8uWwJo1wOuvA2vXyh2N/MaOHYvU1FQsW7YMSUlJ6NChA44ePQpHR0cAQFJSUqVrGgxN795FkwpGj6YxpZUr5Y2JMUNiFLWS0tKovEKrVsCpU9THzNsAyEdJ+zHExtJGQtbWlBwWLy5aTMeYCTKdWkkbNtCiNADo14+TAivy0kuAnR2V1ggKonpMjLGKGXSL4dAhSgTW1vQfnhOCMiipxVCWadNoDGr9erkjYUxyxt1iEARg5kyqzGluzklBCfz9/eHi4oIePXrIHUqF5s2jbVefPwd27aK/JcZYEYNrMSQmUo2jsWPljoSVR+kthgKRkVSPKTLS+MqhMFYO42wxhIcD69bJHQUzBi4uwL17lBRmzwb+9z+5I2JMGQxmuuqpU1TsbsQIujEmJicnoHFj2igoNxewsJA7IsbkYzAthh07gIAAuaNgxmruXNq69YsvgOHD5Y6GMXkpfozh2DGgZ0+ai84Mg6GMMZQlJYV2k3NyorGsQYPkjogxURnHGMOkSbT3MmNSsLOjekuhoVRiIydH7ogYk55iE0NMDH1NTqaN4hmT0iuv0H4deXm0WlqkvYcYMwiKTAzx8bTXcmys3JGwqjCUdQy6qlWL1jqEhtLKacZMhWLHGBITgebNpbwiE4shjzGUJzeXCjOuWkXTXBkzUIY3xhAZSRvBP37MSYEpi5kZ0LUrL4RjpkFRieGll2jxGv/nY0qjUgHLlwNNm9Lf6ZkzckfEmP4oIjGkpABTp9JA3+jRprupDlO+WrWofHenTjzuwIyXIhJDfj4N8nExM2YIRo8GLC0BW1vg/Hm5o2FMfLInhvPnae749u20by9jhqBOHepOcncHbt6UOxrGxCVrYnj4EBgwALh+Xc4oGKseNzeaWt2mDX1lzFjImhhsbWkGUseOckbBWPW1agUkJABNmgBhYXJHw5g4ZEsMr78OHD1Ku68x42BsC9x01aIFlex+9VW5I2FMHLItcNu6lcoOtG4t5rMyJTDGBW66yMwE7t6l2XXt2skdDWNl0mnOp+T7McTF0TS/yZOlvjJj+mVtDaxdCzx7BmzZInc0jFWf5C2GhQupBtKePWI8G1MiU20xFHfwINC/P+9FzhRHpxaDLF1J+fm0UIgZJ1NPDIIAODjQFOx+/eSOhjEtyqqVlJxMO2M9ecJJQWrr16+Hk5MTrKys4OrqirNnz5b72P3792Pw4MFo0qQJbGxs4OnpiWPHjkkYreFTqWiswdmZ95Fmhkmyt+g6dagqJe+lK63du3djzpw5WLx4McLDw9GnTx8MGzYM8eVMvD9z5gwGDx6Mo0ePIiwsDP3798drr72G8PBwiSM3fGfOAMuWyR0FY1UnSVdSZiZ9iqpTp7rPwKrL3d0d3bp1w4YNGwqPtWvXDq+//jpWrFih03O0b98eY8eOxWeffabT4029K6mkAweAXr1orQNjMlNOV9J77wH/939SXIkVl52djbCwMHh5eWkd9/LywoULF3R6jvz8fKSnp6NRo0blPiYrKwtpaWlaN1ZkyRLgr7/kjoIx3UmSGL7+GvjySymuxIp7+PAh8vLyYG9vr3Xc3t4eycnJOj3HN998g6dPn2LMmDHlPmbFihVQq9WFNwcHhxrFbWyuXqWyGRUM7TCmKHpPDNu3A/XqUaE8Jg9ViTrmgiCUOlaWnTt3ws/PD7t374ZdBf+AixYtgkajKbwlJCTUOGZjc+gQ8O23ckfBmG70mhgyM6kZzdUn5WFrawszM7NSrYOUlJRSrYiSdu/ejSlTpuA///kPBg0aVOFjLS0tYWNjo3Vj2ubNA/7zH2DfPi4vz5RPr4nB2hq4fZu2RGTSs7CwgKurKwIDA7WOBwYGomfPnuWet3PnTrzzzjv47bffMHz4cH2HaTJiY4Hp04HUVLkjYaxieiuJERMDrFkD+Pvr6wpMF/PmzYOvry+6d+8OT09PbNq0CfHx8ZgxYwYA6gZKTEzEL7/8AoCSwsSJE7Fs2Ro0b+6BoKBk5OQAtWrVgZWVGvXqAY0b04pec3PaC5nppm1bSgonTwI9elAXK2NKpLfEkJdHu7IxeY0dOxapqalYtmwZkpKS0KFDBxw9ehSOjo549gyIiEjC33/HY/Fi2lPgwIEfkJubi08+mYlPPplZ+Dz16k3C8+c/oVYt4MUXqWS6hQWtYm/YEGjfnjat6daNK+ZWxscH2LyZFnwypkSyVVdl0svNBf7+G/jvf+l2/jzQoAF9em3Tht7wmzWjN/3mzemrlRW1CvLz6fxnzyjhZ2QADx5QV+GVKzTzJiYGUKuBhg3TcPasGqtWaTBkiA06dZL7N1cWQQAiI2lCBq9tYBKTr1bSrl30xrN9u65nMH2KjqZZMadOUWKoXx/w8gIGDQJ69hR3NXp8PBAamobRo9UYOFCD8HAb2NgAY8YAkybR6ncG9O0LDBtGRSUZk5B8ieHaNfoU6eOj6xlMH/74A9i0iT7VN21KG8kMGUK7julT8ZXPT5/aYO9eqqZ7/jwlokmTKFGY8uSl7GyqH1anDrcamKSUW12V6dfJk8Dy5TTQ2acPMGUK0KULlSXRJ39/f/j7+yMvLw8xMTGlSmLExQHbtlH/eoMGgLc38P771G1liry9AScnXt/AJCVPYvj3v6nL4r//1fVpmVju3AEWLAAuXqSBzTlz6I1HapXVShIEYPdu4JtvAI0GGDkS+OgjatWYksxM2vO8fn3Tbj0xScmTGOLigHv3qMuASWfTJkrKnTsDfn40S0guVSmid/gwLYJMSwOmTaNkZkrFFkeNoq69776TOxJmIuRJDNnZXFpbSvfv05vpnTvA1KmAry9Qu7a8MVWnuurWrbQ62MmJWhIDBug5SIV4/Bh4+hSwt5f/342ZBOmrq27YQHPZmTQuX6YB/qQk4IcfaB9tQ31zmTyZklyPHsDgwcD8+TRF1tg1bAiMHk2tPcaUQtQWw6NH9MmVS2DoX2goDSp37UrdSFZWckdUpKb7MRw/DnzyCa2f2LpV3m4xKdy5A1hamt4YC5OF9C0GS0tOClI4d45mGw0eDPzyi7KSghi8vGjdxQsvAOPGAUeOyB2Rfjk6Au3a8YQNphyiJYYLF2il7LNnYj0jK8vFi8CECbTxkTF3PzRrRjufDRsGjBgBrFsnd0T6FRgIDBwodxSMEdFqJbm50TRVU5pRIrXbt4G5c4GhQ2nzI32vS1CCVauo5fDNN1SK46OP5I5IPywsaKFbQgKNOzAmJ9ESQ14eDzzrU34+MGMGFajbuFHuaMpWfIGbmObOpeQwaRLN4Pn8c1GfXhFcXIC9ezkpMGUQrSvJwQE4eFCsZ2MlrVwJpKfT7COlmjlzJiIjI/GXHjY4HjcOWL0aWL8e+O030Z9edubmwKVLVCqEMbmJ1mI4dAhcRVNPrl2jPvYNG4CWLeWORj7vv09VXSdPpvEsLy+5IxLX6NFA//5yR8EY10oyCP3702Ds9u2GMa5Q0+mqFREEaj2EhQGnTxtfnaUOHehDQL9+ckfCjJR001X//W+aPcLEd/QoEB5OXUmGkBT0TaWiBGltTSu9jc0nn1ByYExOonQljR4N9O4txjOxkpYuBcaPpzEcpXn2jPZ6uH+fBsft7Wk+vr7Vrg389BPg4QF8/z1N3TUWKhVNMti7V+5ImCnjriQFCw6mue2RkbQISm7Z2cCXX9KiusRE+r4sZmZpyMtTY/x4Db76ykZv3T3vvw8cO0bdSsYymyc6mvYyGT1a7kiYkZKuiJ6HB82r79NH16diunj3XSqyduCAvHE8fw689Rbwv/9Ry8DKiloGffoA3btT0srPB+7eBf76Czh/Pg1hYWoAGgA2aNaM+s29vcWNKy2N9pnw9qZ1Dsbg+XMaX9i6lXe7Y3oh3RjD1KnAyy+L8UyswPPn9EY8YYK8cezeTXsFHDkCODvTLmzPntHUyjVrqJrrK6/Qm5lG44/AQBekp/cAANy4Abz5Ju0N/cYbQIsWwK1b4sVmY0PdbBs30lReY2BhQa0FrpvE5MRdSQp14gR1lVy6BNStK08MCxcCX30F1KtH+yboOlOm5Kyk/Hxq/fz6K80qWrmSNhQSQ2YmzdhasABYtEic55TbJ59QYpg9W+5ImBGSpsUQE0PTB0Ve7GryTp8GOnaULymsW0dJoUUL4OHDmk2frFUL+Plnai3Ur08JZ8oUceK0tqbNbrZtE+f5lMDdndcEMXnVODFYWdGnG55KKa7r12lvAjncv0+fVhs0oPpMlpbiPG+rVjRm8sIL1If+ySfiPK+vL3DzJg3SGwOVCtiyRe4omCmrcWJ48UUqVVBL1ALepi03F0hOpm065TBkCHX5BAdTqQZd5OdTV063blTkrzxmZlQorm5dYMUKcd7MBw0CGjUynpIsrVrRuA1jcqnx2/nlyzSnnIknLQ148gR46SV5rn35MnVnODvrdo6tLb3hr1xJi/GCgui4Wg20bl368WZm9DhAvDdAFxfgjz/EeS65FSziK286MGP6VuPEcPMm9Ycz8dy7R4nB1lb6axdM+9ywQbfHm5kBqal0v18/2legePXX2NiyuxlffplW+KamijNTqXVrGu8yBnZ2VEnAULdpZYavxonhrbeMa+BPCe7dAzQa+sQttaNHqftIl534Zs6kLiQzM+p6OnmSdpWbMaPoMQVdUWUNohes7p01q+Zxd+5sPJtEZWcD+/YVJVzGpKZTD7JKpVJpNJoyf5acTAOKUpRCMBX37tGAb0aG9NdOTKSZPmlplT+2oGUQG5uFu3ezipXtKFhUkIY33gD+8x+aVlryOV94gb6Ghup2vYrUqUOtrNRUw/+kbWYGvPYa3a/p68JMS2gobfjk5FT2z9VqtQ2AdKGSdQo6rWNQqVQ2oGWsjDHGDJtaEIQKP3LomhhUGo0mv6yfCQKtYShr9kqPHj2qtWmLoZyXlpYGBwcHJCQkVLm8dEXXPH0aWLYM+PNP8WLV9byuXYGUFGo5VHZey5bUWjx2LAsuLlno0qWg+yMJgBu2b4/Eu+82R04OPb5kozM/n2ocvfACrZKuaqzFbdsGzJnTA48e/QUzsyqdKunfjS5/M1lZtO7jnXdoJXRN46zJuUp7bcS8Xk3Ok/qaYr02arVaDR1aDDp1JVX0JIIA5ORo/wEXMDMzq1Y9fkM5r4CNjU2Vz6/oms2aUdeItXXphKvv37FDB+C//6WFaCpVxecdOAD07QsMH05/Aw8fUrwFix0nTKgPgM7dupVKWBT39tv0dcYM7Z9V53dMSABq1TJDw4aG8XdT0d/MjRs0BXzKlJq/LjU9V2mvjdjXM6TXFKj5a1NZS6FAjQef160DevUq+2czZ86s1nMaynk1UdE11Wqq/VPWYKq+f8fx4+mT/LFjlZ/3yis0cyo3F2jblo7l5gIXL9J9Cwt6jCBQSYzi4uOBXbvo/mefVS/W4q5dA+ztjePvxsmJ1ng0aCDe9Yz9/5QccRr1ayMIgq63Mj18KAgREeX91LhpNBoBgKDRaER93gcPBMHGRhCuXRP1aXWSlycItWsLQq9eup9jaSkIgCBYWwtCXJwgJCQkCACEhISEMh8fEECPBwRhyxZx4m7aVBDGjBHnufRJl7+ZQ4cEwd5ewqAUQl//n4yBrq/N8+eVPpVO7/c1bjE0bizfCl25WVpa4vPPP4elWDUj/mFrS4XrChaBSalWLVqPEBwMPHqk2znPn9Nq3cxMGndo355ej9zcotclJ4dqLzVoULQyeulS2r+5pq5cobENX9+aP5e+6fI34+lpmlPA9fX/yRjo+tp06yZOOZUaV1eNjgbmzwd+/93wpwkqycSJtJ/xihXSXzsxkUqdDBpU1KWki1OnaG+EJ0/SABTtx1BSkyZASEj5U+qqyseHxkWePDGO0iwLFlCy+9//5I6EGZrwcPqQVsEaKGmqq9rYUCXIqs4EYRXr2xeIi5Pn2s2bU2I6fhwICND9vH79aJZSQUvDzY0SjIMDbejzxRc04yYlRbyk8PgxfSiZMME4kgJAyXX1armjYIaoa1dxFsbyfgwKdfEi8OmnwJ49NENIanl59IaemUklK6pSnqPkfgz65O1NK67//psK6RmDzp2BadPEWRHOWAnS7eC2YwetgGbiadsWuHOHuhTkYGZG6ymeP6dP+0+fVn6Ov78/XFxc0EOieuEHD1IJj48/Np6k8Pw5lRopXlaEMamJ0mJwdweWLwe8vMQJipFx42gBmK4F7fThxAlap9CiBfVf1qtX+TlStBhiY2m/ipdeon2mjUV0NJU9Dw+nf3vGdHX3LjBpErB/vwLGGAAaSOSkIL5Zs4CdO3X7tK4vAwbQH9rdu1SGOzpavlgK3L8P9OlDiwCrMjhuCK5eBb7+mpMCqzobG1pTZmVV8+cSJTHk5pYud2CMHj9+DF9fX6jVaqjVavj6+uLJkycVnrN//34MGTIEtra2UKlUiIiI0Pl6PXpQwj18uGZx19SwYfSpPDsbcHUFfvhBvliiomjL05wc4Px5ZXYhrV+/Hk5OTrCysoKrqyvOnj1b7mOTkpLg4+MDZ2dn1KpVC/7+c4y6W7Yqr83+/fsxePBgNGnSBDY2NvD09MQxY/skUExVXptz586hV69eaNy4MerUqYO2bdtiy5bvsGyZODsuipIYvv8e6N9fjGdSNh8fH0RERCAgIAABAQGIiIiAbyWT558+fYpevXph5cqVVb6epSV9Ajh5Uv49tTt0oAHe7t2pD9zTk8ZApOTvTzPg6tensZeWLaW9vi52796NOXPmYPHixQgPD0efPn0wbNgwxMfHl/n4rKwsNGnSBIsXL0b79p1x4wYwdqzEQUukqq/NmTNnMHjwYBw9ehRhYWHo378/XnvtNYTLscBHz6r62tStWxezZs3CmTNnEBUVhU8//RSffvopNm3aJE5Auq6Eq2gpXVqaINy6VemKO4MWGRkpABCCg4MLjwUFBQkAhBs3blR6/u3btwUAQnh4eJWuGx0tCF27CkJUVFUj1p8tWwShcWNBMDMTBB8fQbh7V/vnYq9gDQ8XBDc3Wik9apQgZGaK8rR64ebmJsyYMUPrWNu2bYWFCxdWem6vXn0FT88PhWfP9BWdvGry2hRwcXERli5dKnZoshPjtfH29hYmTJhQ2cOkWfkM0Cc4ObahlFJQUBDUajXc3d0Lj3l4eECtVuPChQt6u26bNjR9cdUqvV2iyiZPpj0jPvwQOHSIZlANGFB2NdiauHCBWqK9etGiu0OHqHBfnTriXkcs2dnZCAsLg1eJATcvLy+d/kaioqi7Tow+YqWp6WsDAPn5+UhPT0cjJfYf1oAYr014eDhOnLgAW9u+osQk2pKg6dONZ8/dsiQnJ8POzq7UcTs7OyTruVP4+++B3buBc+f0epkqsbCgbUCTk2nl+99/02yaRo2oXDRAi8+qIjubus1mzaJFdgMGALdvA//6Fy32K9i8RqkePnyIvLw82Nvbax23t7ev9G/k8WPa0lOXnfMMUU1emwLffPMNnj59ijFjxugjRNnU5LVp0aIFLC0t0b17d3TrNhOdOk0VJSadym7rok0befYorik/Pz8sXbq0wscU1DJXlbF5sSAIZR4XU716NFNl2TJajawkdesCfn50Cw6mOi1Hj9LP2ral2F9+GbC3pwVzTZrQCmVLS1oB/fAhtQZiYqikhUZD9bf69AE++AAYOLDsPaOVrOTfgy5/I3/+SclvyBA9BqYA1XltAGDnzp3w8/PDwYMHy/yAZgyq89qcPXsWGRkZCA4OxsKFCzFtWmsAb9c4FtESw0cfifVM0po1axbGjRtX4WNatmyJK1eu4P79+6V+9uDBg1KZXh9GjAAuXaJtMpX6gSkszB/nz/ujXj0aKV+/nrYavH2bymKHhlJJjMxMSij5+bTnRPPm9Em5Tx/aM7pDB5l/kWqytbWFmZlZqU95KSkpFf6NZGZS6ZFu3fQdoXyq+9oANDA7ZcoU7NmzB4MGDdJnmLKoyWvj9E9tmY4dO+L+/fvw8/PD228rKDHk59MG5t7eZe/mplS2traw1aGp4+npCY1Gg9DQULi5uQEAQkJCoNFo0LNnT32HiRdfpKmiJ07QG2jBfslKMnPmTMycObNwgdubb5beh8GYWVhYwNXVFYGBgfD29i48HhgYiFGjRpV7Xn4+JUzd15oanuq+Njt37sTkyZOxc+dODB8+XIpQJVfd16YkQRCQlZUlTlC6jlJXNtSdmCgItraCcPu2zoPoBmfo0KFCp06dhKCgICEoKEjo2LGjMGLECK3HODs7C/v37y/8PjU1VQgPDxeOHDkiABB27dolhIeHC0lJSdWKoVcvQVi5ska/ht6Zcl39Xbt2CbVr1xa2bNkiREZGCnPmzBHq1q0rxMXFCYIgCAsXLhR8fX21znFyChe++CJccHV1FXx8fITw8HDh+vXrcoSvV1V9bX777TfB3Nxc8Pf3F5KSkgpvT548ketX0Juqvjbr1q0TDh06JMTExAgxMTHCypVbBWtrG2Hx4sWVXUqn93vREoMpSE1NFcaPHy/Ur19fqF+/vjB+/Hjh8ePHWo8BIGzbtq3w+23btgmgciJat88//7xaMeTkCIKVlSAcP17930PfTDkxCIIg+Pv7C46OjoKFhYXQrVs34fTp04U/mzRpktC3b9/C74OChDL/PhwdHaUPXAJVeW369u1b5mszadIk6QOXQFVem7Vr1wrt27cXrK2tBRsbG8HRsavQocN6IS8vr7LL6PR+z9VVDVBQEPDGG/RViYu8pKyuaujefpvGVBYvljsSZiKkq5VUIDSUpiuWtVcxE4+bG/DJJzRrJz1d7mhYda1YQWtUOCkwpRE1MXTqBPz6q3IXIBkLMzNg/Hja7OZf/6JaVcyw3L1Ls7G6dJE7EsZKEzUxWFlRiWamf40a0cKysDBg+3bjntFijJYvBwIDi/a/Zqy6UlLErzIs+maIycnAwoXyF30zBa1aAWvXUrfS6dNyR8N09fXXtPhPjE3bGTt9mvYJF5PoiSEjA7h+nbs3pNKuHQ1CDxwIXL4sbyxS7+BmiK5epdZe69ZUVoSxmho9GqhCNX+d8KwkI5GURGMOgYG0AE5OPCupfEOGUFHAzz6TOxJmoqSflVTc33/r65lZWV54gVoOH35IYw5MeXr1AmbPBpYskTsSxiqml8SQnExluCMj9fHsrDxdu1L/dWgoDW4yZcjKAlauBKZOpX8jQysKyJTr5k3ayVBsekkMTZtSpUgXF308O6vIwIFUr6pevaLy10w+OTn0AenAAepGatZM7oiYMdm/H1i9Wvzn1VtXkqOjvp6ZVaZ/f2DQIMDR8TFatvSFjY1ue1Tn5ORgwYIF6NixI+rWrYtmzZph4sSJuHfvnnTBG5nt2ylRBwdzUmDiW7AA2LNH/OfVW2IAaG9guTeyN1UdOwJBQT5IT4/AW28FYPnyyveozszMxKVLl7BkyRJcunQJ+/fvR0xMDEaOHClh5MZjxgzqVg0KkjsSxqpGr7OS/P2pfAPPXpReVFQUXFxcEBwcjOPH3aFWA1euBGPLFk/cuHEDzs7OOj3PX3/9BTc3N9y5cwcvvviiTueY+qykZ89okHnoUFprYqy7sjGDJO+sJIBaDJwU5FF8j+olS6j0Qt26HrCwUOPQId33qNZoNFCpVGjQoEG5j8nKykJaWprWzVRFRgJHjlBy6NmTkwLTn+3baRMsfdBrYgCA1FRelSuHkntUv/IK8O23gKWlHc6dS8aaNZU/x/Pnz7Fw4UL4+PhU+Ml/xYoVUKvVhTcHBwcxfgWDc+4cbaS0fz/9p1XiZkrMeGzfTouJ9UHviSEggFoOTBx+fn5QqVQV3i5evAig9B6yZmaAvb0AGxsVMjOBceOA+Piyr5OTk4Nx48YhPz8f69evrzCmRYsWQaPRFN4SEhJE+V0NxdOn9DqOGkUrmn/7Te6ImCkICKAtf/VB75twjh9PNyYOMfaoHjDAHhMmAHPnAgcPUnP022+LHpOTk4MxY8bg9u3bOHHiRKXjBJaWlrC0tKzW72PoEhKAnTuBvXtpoLl2bbkjYqzmDGh3ZgaIt0d17drAunXAhQtAgwbUF/7dd4CbWw7Gjx+Dmzdv4uTJk2jcuLGefyPDlJ5OLYWXXqLKllOnclJgxkPvXUkF+vXjUg1SateuHYYOHYpp06YhODgYwcHBmDZtGkaMGKE1I2ny5LaoV+93zJsHxMbmolGjtxAaehEbNuxAXl4ekpOTkZycjOzsbBl/G+UQBCAmhlpb8+cDUVH0t92okdyRMVMREgJ8/rl+ryFZi2HBAtrCkElnx44dmD17Nry8vAAAI0eOxLp167QeEx0dDY1Gg3feAa5fv4usrEO4dw945ZUuWo87efIk+vXrJ03gChUaCjx+DIwcSSv7LSwAblAxqWVnU2tVn7i6KitFo6FSDt260XL7oCCqwaQrY1vHEBJCLQVvb1qb4+HBq5iZwdJpHYOkYwyZmdSPPX8+7fbGlEmtpq/x8VTn/ehRqgiank7VW5s2NY3tWy9fpnIDOTmUGOLiABMdY2cmRtLEoNHQ4p/p0zkxGIouXYDNm4GLF6lK6PDhwPvv0/qU2bO1+9b9/f3h7++PPAPevi8+nqqfDh1KraSHD4GNG+WOijFpcVcSq5LUVMDcHHjjDWo9fPghcOYMfZou2CDIELuSdu4E2rSh6rSbNtE+utOn8y5rTDnu3qWJDmfP1mjxpPK6kgoIAr3B6DDrkilMwWDrn3/SNq5ffAGcPAl88w0tZIyNpYVzShcXR7ve/f477U/+9CltR3v1KtC8OVBLsvl6jOnGzo4m8UjxvilLi+HgQWDSJKCCCtDMAJ09S4u89uxJw549akycqEH//jZo3Ji6pOzs5Omjz82lbszAQBob2bABGDyYylcsXQrk51OxR8ZMgPxF9MozfDjVp2fGpU8f2pj8hx/o+wEDgM6dqTVx/jxVGt27F5gwgWq8/PorkJYm3jaw9+7RQPGuXdTsnjiRWjMNG9IYyfLlQKdOdHz+fBrv6t6dkwJjJcmSGMzNgbZt5bgyk4KZGX319qYV1adOUffSvn3AsGFAkyY0wLtsGc16cnamT/MNG9LUWHd36tLx9qatCydOpOTxzjvUVfXWW7SwrHdv4NIlwN6eWivNm9P5c+bQegMLC8DVFTh0iHZPu34dcHICfHzke20Yq46wMOqCl4qsJTFCQoD164Gff5YzCiYVDw/6+t139PXmTfoaGwu8+CL9HXTuTG/89vZAy5Y0e83cnPr8VSq637QplfEYOpRKUqxZA3h6AjduUJJJTqbn3byZvvbvL+VvyZi4Hj+m/zuhodKVcZd1VlJUFDX7ly4V+5mZnAxxVhJjSqbRFK0vqiGdxhh4uioTHScGxhRLuYPPJT17xntDGwN/f3+4uLigB2/bx5hBU0SLISSEZirdvcsroo0BtxgYq7kHD6iK74YNQP36oj2tche4leTuTqUHGGOMkVq1aN2NHIstFdFi0LqIQLNPmOHiFgNjimU4YwwFnj8vmsvOGGNMHopKDFZWwO7dtA8AY4yZoowM8aoBVJeiEgNAK1S55j1jzFT9+COVlpGT4sYYCmzdSqUVJk2S8qpMDDzGwFj15eZSq6FBA708veGNMZTEg9CMMVNjbq63pKB7DPJevnyTJ8sdAasqY9jBjTGm4K6kAnl5tH/D66/z5imGgruSGKu6fv2Azz6jcvV6ZPhdSQCthp4yhWrtM8aYsfL1pe1llUDxLQZmeLjFwJhiGUeLobi8PCAzU+4oGGPMuBlUYvjyS+DVV+WOwnA8fvwYvr6+UKvVUKvV8PX1xZMqbLT93nvvQaVSYfXq1XqLkTFTduwY8NNPckdRmkElhpkzafEH042Pjw8iIiIQEBCAgIAAREREwNfXV6dzDxw4gJCQEDRr1kzPUTJmuu7fBxIS5I6iNMVOVy1Lo0Z0Y5WLiopCQEAAgoOD4e7uDgD48ccf4enpiejoaDg7O5d7bmJiImbNmoVjx45h+PDhUoXMmMmZOFHuCMpmUC2G4mbNAjZulDsK5QoKCoJarS5MCgDg4eEBtVqNCxculHtefn4+fH198fHHH6N9+/Y6XSsrKwtpaWlaN8aY4TLYxDBwIBfbq0hycjLs7OxKHbezs0NycnK553311VcwNzfH7Nmzdb7WihUrCscx1Go1HBwcqhUzY6YgKYn2ulcyg00M3t6Am5vcUUjPz88PKpWqwtvFixcBAKoyaooIglDmcQAICwvDmjVr8NNPP5X7mLIsWrQIGo2m8JagxE5TxhQiIgJQ+nwOgxpjKM+//gUMGkQ7wRm7WbNmYdy4cRU+pmXLlrhy5Qru379f6mcPHjyAvb19meedPXsWKSkpePHFFwuP5eXlYf78+Vi9ejXi4uLKPM/S0hKWXBKXMZ0MG0Y3JTOKxJCaCjx9KncU0rC1tYWtrW2lj/P09IRGo0FoaCjc/mlahYSEQKPRoGfPnmWe4+vri0GDBmkdGzJkCHx9ffHuu+/WPHjGmEEwisTw7bdyR6A87dq1w9ChQzFt2jT88MMPAIDp06djxIgRWjOS2rZtixUrVsDb2xuNGzdG48aNtZ6ndu3aaNq0aYWzmBhjFXv0CDh7Fhg1Su5IdGOwYwzl+f13ICZG7iiUYceOHejYsSO8vLzg5eWFTp064ddff9V6THR0NDQajUwRMmYa/voL+PhjuaPQndHVSvL2Bt58E5gwQe5ITBfXSmJMsXSaVWIUXUnF/f673BEwxphhM7qupOKiooBbt+SOgjFmimJigAMH5I6ieow6MXz1FbB2rdxRmA5/f3+4uLigR48ecofCmOyCg4ESQ3oGw+jGGIoTBN43Wg48xsCYYhnffgxVVTwpZGYC16/LFwtjjBkKo04Mxe3ZA7z2mtxRMMaM2dq1tMeCoTOZxDBxInD5stxRMMaMWXq6cewyadRjDBWJjAScnQEzM7kjMT48xsCYYvEYQ3ny8oDu3YETJ+SOhDHGlMckE4OZGRAfDwweLHckjDFDFhMDODhQIU9jYpKJAQCKFyjVaICbN+WLhTFmmBwdgS+/BBo2lDsScZlsYijuxx+5tpIYeIEbMzWWloCvL1DLyN5JTXbwubi8PODZM6BePbkjMQ48+MyMWWgo0KwZ0KKF3JFUCw8+68rMTDsphIQADx7IFw9jTLk+/RTYuVPuKPSLWwxl6NEDmDIFmDFD7kgME7cYGFMs0yy7LYaQEOPrM2SMMV3x218ZiieFvDxg3TogI0O+eBhj8tm0CRg/Xu4opMWJoRIZGfSHkZwsdySMMTn07g28/bbcUUiLxxiY6HiMgTHF4llJ+nDvHjBpEvD8udyRKA+vY2DGIDYW+O47uaOQFyeGKlKpAHMesi/TzJkzERkZib/++kvuUBirtsRE4MwZuaOQF3clMdFxVxJjisVdSVKIiKBVkDxriTHDlJYGnD4tdxTKwomhhpydgTVruJwGY4bqxAnayIsV4cRQQ3XqAKNHF30vCMDJk/LFU9zjx4/h6+sLtVoNtVoNX19fPHnypNLzoqKiMHLkSKjVatSvXx8eHh6Ij4/Xf8CMyeD114G4OLmjUBZODCKLiwOGDqX9HuTm4+ODiIgIBAQEICAgABEREfD19a3wnNjYWPTu3Rtt27bFqVOncPnyZSxZsgRWVlYSRc2YfpVVZl+lU8+76eDBZz3IzgYsLOSNISoqCi4uLggODoa7uzsAIDg4GJ6enrhx4wacnZ3LPG/cuHGoXbs2fv3112pfmwefmZJ9/TWwbx8QHCx3JLLgwWe5lEwKvr7A5cvSxhAUFAS1Wl2YFADAw8MDarUaFy5cKPOc/Px8HDlyBG3atMGQIUNgZ2cHd3d3HDhwoMJrZWVlIS0tTevGmFLNmwf88YfcUSgbJwYJvPACYG0t7TWTk5NhZ2dX6ridnR2Sy6nvkZKSgoyMDKxcuRJDhw7F8ePH4e3tjTfeeAOnK5i2sWLFisJxDLVaDQcHB9F+D8ZqKjqaFqYWKFlmn5XGiUECq1YBL79c9H1ycvX3e/Dz84NKparwdvHiRQCAqoyOU0EQyjwOUIsBAEaNGoW5c+eiS5cuWLhwIUaMGIGNGzeWG9OiRYug0WgKbwkJCdX75RjTgwULgLVr5Y7CsPAaXhnMm0efWDZtqvq5s2bNwrhx4yp8TMuWLXHlyhXcv3+/1M8ePHgAe3v7Ms+ztbWFubk5XFxctI63a9cO586dK/d6lpaWsLS01CF6xqS3Zw+1EpjuODHI4Mcfqz8LwtbWFra2tpU+ztPTExqNBqGhoXBzcwMAhISEQKPRoGfPnmWeY2FhgR49eiA6OlrreExMDBwdHasXMGMSyssDNm4E3nkHqFuXjtWuLWtIBom7kmRQt672mMO1a8CgQUBurnjXaNeuHYYOHYpp06YhODgYwcHBmDZtGkaMGKE1I6lt27b4/fffC7//+OOPsXv3bvz444+4desW1q1bh8OHD+ODDz4QLzjG9CQjgxJD8TEFVnWcGBSgSROgf3/xm7s7duxAx44d4eXlBS8vL3Tq1KnUNNTo6GhoNJrC7729vbFx40asWrUKHTt2xObNm7Fv3z707t1b3OAY0wO1Grh6VXtMj1Udr2NQqMOHgVatgPbt5Y6k6ngdA5PKzz/TrKMvv5Q7EoPB6xgM2a5dQFCQ3FEwpmyOjkCJuRJMBNxiMCDXrxtGC4JbDExfrl0DOnSQOwqDxi0GYxIRAXTqBDx6JHck5eMd3Jg+JSbS/4GYGLkjMX7cYjAgKSlA8cXM+fk0k0nuukwlcYuB6UvJ/wOsyrjFYGxK/ofw8wNGjpQlFMb0buZMYPdu7WOcFKTBC9wM2MyZwPjxckfBmH706EEz85j0uCvJyAwaBAwbBsyfL18M3JXEqiooiLqJRo2SOxKjp1NXErcYjMyXX9Ie1MXl5gLm/C/NFOzSJSA2lhODUnCLwQTY2wNbtgAjRkhzPW4xsIqkpACRkUC/fnJHYpJ48JmRffuo5EZxSth6lJmmY8eowjBTLm4xmKCnT4EGDYBTp4BevcR/fm4xsAJPngBhYcDAgXJHwv7BLQZWtrp1qb5MyaSwaxfw7Fn1n5cXuLGSTp8Gpk6VOwpWVdxiYACA9HSgXTvg6FFaXVoT3GIwTfHxwIkTtBcCUyxuMTDd1a8P3L2rnRTy8oA33wSiouSLixmOqChg+3a5o2Bi4MTAylWrFtC6demN0+PjgexseWJiynD4MDBnjvaxIUOAP/6QJRwmMk4MrFwqFfDVV4CDg/bxvn35k6EpSUsDiu3lBIBKU7z0kjzxMP3jMQZWZcnJgK2t9qK5L76ghXWTJ/MYg7GZMgXIyuIPA0aCVz4z/WjatPQxR0faorS4v/+mJDJ0qDRxsZqbPx8YPhwYMKDo2DffUCVfZjq4K4mJwte3dAI4exbw9y/92OvXpYmJlS8vD9ixA8jJ0T7eqBFgba19rEEDOs5MB3clMdFV1JV07RrQuTO1JIq3MB49oo3czcwkDtYEZGTQ6926ddGx5GTA1RU4dw5wcpIvNiY5nq7KpKXLArcOHYCkpNLdTj17Aps3ax/LygIePtRDoEYsIYFWthe3bRvwxhvax5o2pR3ROCmwsnCLgYmuOoPPsbE0eF2nTtGxL7+kaZFBQdqPPXkSeOUV025d5ObSquKSpSa6dwfGjQM++qjoWF4etRrUamljZIqkU4uBEwMTnVizkjIzqdZO8TLiqalAixZAaCjQsWPR8QMH6Ovrr1f7cop14ADNACteHTcigkqaJCbSGECB1FRarKi07V6ZYnBXEjNs1tal95Zo3Jjm1BdPCgBw+zYQF6d97OZNetO8f1/7+JEjwJkzpa93+3ZNIy7b3bulj507V3oxWFISlUgv+XvExAC3bmkf69KFXofiSQGg14eTAqspTgzM4JT1xjd3bumVuE5OwE8/ld4n+Px5qvhZ3LVrtGArJUX7+NSp2t0yAL3Rt2lDA7glY5g1S/tYYiItEIyO1j5+5gzw55/ax+ztgbVrqUVU3P/7f6V/N4A3X2L6w11JRuzx48eYPXs2Dh06BAAYOXIkvv/+ezQo+TGzmIyMDCxcuBAHDhxAamoqWrZsidmzZ+P999/X+bqGusDtzh1aj1FcRAQlIheXomPZ2TSgO3kyULt20fErV2i+f5cu2s9x8ybw8sv6ipqxKuExBlM3bNgw3L17F5s2bQIATJ8+HS1btsThw4fLPWfatGk4efIkNm/ejJYtW+L48eP44IMPsG/fPozScd9FQ00MjJkATgymLCoqCi4uLggODoa7uzsAIDg4GJ6enrhx4wacnZ3LPK9Dhw4YO3YslixZUnjM1dUVr776KpYvX67TtTkxMKZYPPhsyoKCgqBWqwuTAgB4eHhArVbjwoUL5Z7Xu3dvHDp0CImJiRAEASdPnkRMTAyGDBlS7jlZWVlIS0vTujHGDFdVWgzMgKhUqk8AvCMIQpsSx2MAbBMEYUU551kA+BHARAC5APIBTBUE4dcKruUH4PMyfqQWBIGzBGMGhlsMBkalUvmpVCqhklv3fx5eVtZXlXO8wGwAHgBGAnAFMB/AepVKNaiCc1YAUJe4NQGQXrXfjjGmBNxiMDAqlcoWgG0lD4sD4APgW0EQGpQ4/wmAuYIgbCvjuesA0ADwFgThSLHjmwG0EASB66QyZgJ4JrSBEQThIYBKKwipVKogAGqVSuUmCELoP8fcQZ/myxtkqP3PrWSR5Txw65Ixk8H/2Y2UIAhRAAIA/KhSqTxUKpUHaOzgv4IgFC63UqlUN1Qqlfc/56QBOA3ga5VK1U+lUjmpVKp3QOMNv0v+SzDGZMEtBuM2HsBaAMf/+f4QgBJrc+EMakUUGAcaM9gBoBGAOwAWA9io10gZY4rBYwyMMca0cFcSY4wxLZwYGGOMaeHEwBhjTAsnBsYYY1o4MTDGGNPCiYExxpgWTgyMMca0cGJgjDGmhRMDY4wxLZwYGGOMaeHEwBhjTMv/B5SVs/IBzLo/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords = [(sol[k][0], sol[k][1]) for k in range(npts)]\n",
    "list_plot(coords, size=1, figsize=4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.3",
   "language": "sage",
   "name": "sagemath"
  },
  "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.10.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}