{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taking Derivatives with Vandermonde Matrices\n",
    "\n",
    "Copyright (C) 2020 Andreas Kloeckner\n",
    "\n",
    "<details>\n",
    "<summary>MIT License</summary>\n",
    "Permission is hereby granted, free of charge, to any person obtaining a copy\n",
    "of this software and associated documentation files (the \"Software\"), to deal\n",
    "in the Software without restriction, including without limitation the rights\n",
    "to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n",
    "copies of the Software, and to permit persons to whom the Software is\n",
    "furnished to do so, subject to the following conditions:\n",
    "\n",
    "The above copyright notice and this permission notice shall be included in\n",
    "all copies or substantial portions of the Software.\n",
    "\n",
    "THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n",
    "IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n",
    "FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n",
    "AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n",
    "LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n",
    "OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n",
    "THE SOFTWARE.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as pt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are a few functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "if 1:\n",
    "    def f(x):\n",
    "        return np.sin(5*x)\n",
    "    def df(x):\n",
    "        return 5*np.cos(5*x)\n",
    "elif 0:\n",
    "    gamma = 0.15\n",
    "    def f(x):\n",
    "        return np.sin(1/(gamma+x))\n",
    "    def df(x):\n",
    "        return -np.cos(1/(gamma+x))/(gamma+x)**2\n",
    "else:\n",
    "    def f(x):\n",
    "        return np.abs(x-0.5)\n",
    "    def df(x):\n",
    "        # Well...\n",
    "        return -1 + 2*(x<=0.5).astype(np.float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5f98d5dd30>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VHXe/vH3J50EAoQk1EBCJ7QAQxEVRUVBlAA2wIIVH11XH3VVLKuuBXFdn8UCKpYVy1IEUVQUAWEVBST0EEpCDy2Nml6+vz9m2F8SAwnMJGfK53Vdc2XmzDkz96HMnTntK8YYlFJKqdP8rA6glFLKvWgxKKWUqkCLQSmlVAVaDEoppSrQYlBKKVWBFoNSSqkKtBiUUkpVoMWglFKqAi0GpZRSFQRYHeB8REZGmtjYWKtjKKWUR1m7dm2WMSaquvk8shhiY2NJSkqyOoZSSnkUEdlbk/l0U5JSSqkKtBiUUkpVoMWglFKqAi0GpZRSFWgxKKWUqsAlxSAiH4lIhogkn+F5EZE3RSRNRDaJSO9yz40XkVTHbbwr8iillDp/rvrG8DEw9CzPDwM6OG4TgHcARCQCeA7oD/QDnhORxi7KpJRS6jy45DwGY8zPIhJ7llkSgU+MfRzRVSLSSESaA5cCi40xOQAishh7wcx0RS5ljZLSMg4dL2Bvdh6Zpwo4WVDCyYISikvLCPT3I9BfaFgvkOjwEKIbBNM6IpQGIYFWx1ZKOdTVCW4tgf3lHqc7pp1p+h+IyATs3zZo3bp17aRU58wYw66sXNbszmH9vmOs33+UXZm5lJSd21jiMRH16NwsnISYRgxs14TuLRsS4K+7wJSyQl0Vg1QxzZxl+h8nGjMdmA5gs9nO7VNHuVRZmWHNnhwWbTnCT9uOsCc7D4BGoYEkxDTi8i5NiW0SSuuIMJo1DKFBSAD1gwMI8vejuKyM4lLD0dwiMk4WcuREAbuzctl66ARbD51gccoRABoEB3Bxx0iGd2/BZZ2jqRfkb+UqK+VT6qoY0oGYco9bAQcd0y+tNH15HWVS5+jAsXzmJqUzd91+9ufkExTgx8B2TbjrojgubB9JXGQYIlV1/f8X7OdPcADUDw4gJiL0D89nnypk5a5sfk3LYnFKBgs3HyY0yJ+h3Zpx64A2JMQ0qvY9lFLOEftmfxe8kH0fw7fGmG5VPDcceAC4GvuO5jeNMf0cO5/XAqePUloH9Dm9z+FMbDab0Wsl1Z3kA8d57+ddLNx8iNIyw8B2TbjB1oqrujYjNKj2frcoLTOs3p3NNxsP8c3Gg5wqLKF7y4bccWEsI3q20E1NSp0jEVlrjLFVO58rikFEZmL/zT8SOIL9SKNAAGPMu2L/Fe9t7DuW84A7jDFJjmXvBJ5yvNTLxph/Vfd+Wgx1Y8P+Y/xj0XZWpGVRPziAcf1bc+uANlX+pl/bThWWMH/9AT75bQ+pGaeIbRLKnwa3Z1SvlloQStVQnRZDXdNiqF1pGaf4x6Lt/LDlME3CgpgwqC1j+7cm3A2OHDLGsDjlCG8sTWXLwRPERYbx9NVduLxLtG5iUqoaWgzqnOUWlvDm0lQ+XLGbkEB/7rm4LXddHEf9YPe7Ovvpgnj1h23szMzl4g6R/PWaeDo2bWB1NKXclhaDOieLthzmbwu2cPB4ATfZYnh8aCea1A+2Ola1ikvL+GzVXqYsSSW3sIT7B7fnT4PbERygRzEpVZkWg6qRY3lFPPNVMt9uOkTnZg14aWQ3bLERVsc6Zzm5Rbz0bQpfrj9Ax6b1+fv1PUmIaWR1LKXcihaDqtbPOzJ5bO5Gsk8V8fCQjkwY1JZAD9+Ru2xbBk/N38yREwX8+bIO/Pmy9rpzWimHmhaD/o/xQUUlZbzwTQq3ffQ74SGBfPWnC/nT4PYeXwoAgztH8+PDgxjVqxVvLE1l3AerOXQ83+pYSnkUz/8kUOfkyIkCxr6/io9+3c34C9rwzZ8volvLhlbHcqkGIYG8fmNP/nlTT7YcOM6wN37hp21HrI6llMfQYvAhq3ZlM/zNFWw9dII3x/bib4ndCAn03p20o3q14tsHL6Zlo3rcNSOJqcvS8MRNp0rVNS0GH/Hpyj3c/MFqwusF8NWfLmREzxZWR6oTcZFhzLtvICN6tuC1Rdt5YOZ68opKrI6llFtzvwPUlUuVlhkmLdzKhyt2c3nnaKaMSfC5S1yHBPoz5aYE4puHM/mHbezOzOVfd/SlaXiI1dGUckv6jcGL5RWVcN9na/lwxW5uHxjL9NtsPlcKp4kI917Sjo9u78ve7FxGT/uNtIyTVsdSyi1pMXipnNwixk5fxeKtR3ju2nieH9EVfz+9ZMTgTtHMvvcCCkvKuO6dlSTtOev1GpXySVoMXujIiQJuem8lWw+f5L1b+nDHhXFWR3Ir3Vo2ZP79A2kSFsTNH6zmxy2HrY6klFvRYvAy+3PyuOHdlRw8ls/Hd/Tlyq7NrI7klmIiQpl730C6NA/nvs/X8c3Gg1ZHUsptaDF4kbSMU9zw7kqO5xfz+T0DGNgu0upIbi0iLIjP7u5Pn9aNeWjWer5cl251JKXcghaDl9iZeYox01dRUmaYfe8AvU5QDdUPDuDjO/syoG0THv1iI3PW7K9+IaW8nEuKQUSGish2EUkTkYlVPP9PEdnguO0QkWPlnist99wCV+TxNXuychn3/irAMGtCfzo3C7c6kkcJDQrgo9v7MqhDFI/P26TloHye0+cxiIg/MBUYgn0M5zUissAYk3J6HmPMw+Xm/zPQq9xL5BtjEpzN4av25+Qx7v1VFJWUMWvCBbSP1vEIzkdIoD/Tb+vDPZ+s5YkvNxES5O8zJwEqVZkrvjH0A9KMMbuMMUXALCDxLPOPBWa64H193oFj+Yx9fxW5RaV8dnd/OjXTUnBGcIA/793Sh76xETwyewOLU/T6Sso3uaIYWgLlv3unO6b9gYi0AeKAn8pNDhGRJBFZJSIjXZDHJ+TkFnHrB6s5nl/Mp3f1o2sL77oQnlXqBfnz4XgbXVuE86fP17EiNcvqSErVOVcUQ1VnTZ3pSmVjgLnGmNJy01o7rg8+DpgiIu2qfBORCY4CScrMzHQusYfLKyrhjo/XcOBYPh/d3pcerXRHsys1CAlkxp39aBsVxj2fJLFu31GrIylVp1xRDOlATLnHrYAzHRQ+hkqbkYwxBx0/dwHLqbj/ofx8040xNmOMLSoqytnMHqu4tIz7P1/H5vRjvDW2F309cLQ1T9AoNIhP7+pPdHgwd89IYndWrtWRlKozriiGNUAHEYkTkSDsH/5/OLpIRDoBjYGV5aY1FpFgx/1I4EIgpfKyys4YwxPzNrF8eyaTRnXXk9dqWVSDYGbc0Q+A8R/9TtapQosTKVU3nC4GY0wJ8ACwCNgKzDHGbBGRF0RkRLlZxwKzTMUL4ncBkkRkI7AMmFz+aCZV0as/bOfLdQd4dEhHxvRrbXUcnxAbGcaH421knCzgro/X6CW7lU/QMZ89xOw1+3hi3mZuGdCaFxO7IaIXxKtLi1OOcO+nSQzuFM17t/bRcaSVR9Ixn73Iyp3ZPD0/mUEdo3j+2q5aChYYEt+UvyV2Y+m2DP72jX6pVd5NB+pxc7uzcvmfz9YSFxnG2+N66W+qFrp1QBv25+Qx/edddGrWgFsGtLE6klK1Qj9l3NjxvGLu+ngN/n7Ch+P7Eu6jg+y4kyeGdmZwpyieX7CFlTuzrY6jVK3QYnBTxaVl3P/vtew/mse7t/ShdZNQqyMpwN9PeGNsL2Ijw7jv87Xsy86zOpJSLqfF4KZe/m4rv6Zl88roHvSL03MV3El4SCAf3GbDGLj7kzWcLCi2OpJSLqXF4Ia+XJfOx7/t4e6L4ri+Tyur46gqxEaG8c7NvdmZmcvDszdQVuZ5R/cpdSZaDG4m+cBxnvxyMwPaRjBxWGer46izGNg+kr8O78KSrRlMW55mdRylXEaLwY0cyyvivs/X0jg0iLfH9dYjkDzA+IGxJCa04PXFO/gl1bev4aW8h37yuInSMsODszZw5Hgh79zSm8j6wVZHUjUgIrwyujsdoxvw4Mz1HDiWb3UkpZymxeAmpizZwc87MnluRDy9Wje2Oo46B6FBAbxzS29KSg33f7aWwpLS6hdSyo1pMbiB5dszeOunNG7o04pxeg0kj9Q2qj6v3dCTjenH9cxo5fG0GCx25EQBj8zZSOdmDXhxpF4DyZMN7daMey9py79X72P++nSr4yh13rQYLFRaZnhw5nryi0p5e1xvQgL9rY6knPTYlZ3oFxvB0/OT2ZV5yuo4Sp0XLQYLvbE0ldW7c3hpZDfaR9e3Oo5ygQB/P94Ym0BQgB8P/Hu97m9QHkmLwSK/pWXx1k+pXNe7FdfpSWxepXnDevzj+p6kHDrBKwu3WR1HqXOmxWCBzJOFPDR7A20jw3ghsavVcVQtuCK+KXdcGMvHv+1hccoRq+ModU5cUgwiMlREtotImohMrOL520UkU0Q2OG53l3tuvIikOm7jXZHHnZWVGR6Zs4ET+cVMvbk3YcF65XNvNXFYZ7q1DOexuRs5qOc3KA/idDGIiD8wFRgGxANjRSS+illnG2MSHLcPHMtGAM8B/YF+wHMi4tUH8X+wYhe/pGbx3LVd6dws3Oo4qhYFB/jz1tjeFJeU8dCs9ZSUllkdSakaccU3hn5AmjFmlzGmCJgFJNZw2auAxcaYHGPMUWAxMNQFmdxSysETvLZoO1d1bcrYfjFWx1F1IC4yjEmju7Nmz1GmLd9pdRylasQVxdAS2F/ucbpjWmXXicgmEZkrIqc/FWu6rMcrKC7loVnraRwaxCuje+j5Cj4kMaEliQkteGNpKhv3H7M6jlLVckUxVPUJV/kaxN8AscaYHsASYMY5LGufUWSCiCSJSFJmpuddrGzy99tIzTjFazf0JCIsyOo4qo69kNiNpg2CeXj2BvKKSqyOo9RZuaIY0oHy20VaAQfLz2CMyTbGFDoevg/0qemy5V5jujHGZoyxRUVFuSB23fnPjkw+/m0Ptw+M5ZKOnpVduUbDeoG8fmMCu7Nzefm7rVbHUeqsXFEMa4AOIhInIkHAGGBB+RlEpHm5hyOA0/8zFgFXikhjx07nKx3TvEZObhF/+WIjHaLr6/gKPu6Cdk245+K2fL56Hz9t00NYlftyuhiMMSXAA9g/0LcCc4wxW0TkBREZ4ZjtQRHZIiIbgQeB2x3L5gAvYi+XNcALjmlewRjDxHmbOJZXxJQxCXrJC8WjV3akc7MGPD53M9mnCqtfQCkLiDGeNyShzWYzSUlJVseo1pw1+3l83iaeurozEwa1szqOchPbDp9gxFu/ckmnKKbf2kcPRFB1RkTWGmNs1c2nZz7XkvSjebzwbQoD2kZw90VtrY6j3EjnZuE8PrQTi1OOMCdpf/ULKFXHtBhqgTGGJ+ZtwhjDa9f3xM9PfyNUFd15YRwXtG3CC9+kkH40z+o4SlWgxVALPl+9j1/Tsnny6i7ERIRaHUe5IT8/4e/X9wBg4rzNeOImXeW9tBhcbH9OHpMWbuWi9pHc3F9HY1NnFhMRylPDu7AiLYt//77P6jhK/ZcWgwuVlRken7sJPxEmX9dddyqqao3r15qL2kcy6but7M/RTUrKPWgxuNBnq/eyclc2zwzvQqvGuglJVU/K/RLxxLxNlJXpJiVlPS0GF9mbncsrC7cxqGMUN/XVC+SpmmvVOJSnru7Cbzuz+Vw3KSk3oMXgAmVlhsfmbiLAX3hVNyGp8zC2XwwXd4jklYW6SUlZT4vBBWas3MPvu3N49pp4mjesZ3Uc5YHsm5R64CfC43N1k5KylhaDk/Zm5/LqD9u4rHM01+vYzcoJLRvV45nhXVi5K5vPV++1Oo7yYVoMTjDG8OSXmwn082PSKN2EpJx3U98YBnWM4pXvt3FAhwNVFtFicMIXSen8tjObiVd3plnDEKvjKC8gIkwa1Q2Ap+friW/KGloM5ynjRAEvfZdCv7gIxvbVE9mU67RqHMpjV3Vi+fZMFmyscngSpWqVFsN5em7BFgpKypg8urteC0m53G0XxJIQ04i/fZNCTm6R1XGUj9FiOA8/JB/m++TDPHR5B9pG1bc6jvJC/o5rKZ0sKOaFb7ZYHUf5GC2Gc3Q8v5hnv04mvnk4Ewbp5bRV7enYtAH3X9qerzYcZNn2DKvjKB/ikmIQkaEisl1E0kRkYhXPPyIiKSKySUSWikibcs+VisgGx21B5WXdzeTvt5J1qpBXr+tBoL/2qqpd9w9uR/vo+jwzP5lThSVWx1E+wulPNhHxB6YCw4B4YKyIxFeabT1gM8b0AOYCfy/3XL4xJsFxG4EbW7kzm5m/7+eei9vSvVVDq+MoHxAc4M+r1/Xg4PF8/rFou9VxlI9wxa+8/YA0Y8wuY0wRMAtILD+DMWaZMeb0ef6rAI87E6yguJQnv9xEmyah/O8VHa2Oo3xInzaNGX9BLDNW7mHt3qNWx1E+wBXF0BIoPz5humPamdwFfF/ucYiIJInIKhEZeaaFRGSCY76kzMxM5xKfhylLUtmTnccro7pTL8i/zt9f+ba/XNWJ5uEhTJy3icKSUqvjKC/nimKo6ljNKs/KEZFbABvwWrnJrR2DU48DpohIu6qWNcZMN8bYjDG2qKgoZzOfk+QDx3n/l13cZIthYPvIOn1vpQDqBwfw8ujupGacYtqynVbHUV7OFcWQDpS/znQr4A9n5YjIFcDTwAhjTOHp6caYg46fu4DlQC8XZHKZktIyJn65iYiwIJ66uovVcZQPG9wpmsSEFkxbnkZaximr4ygv5opiWAN0EJE4EQkCxgAVji4SkV7Ae9hLIaPc9MYiEuy4HwlcCKS4IJPLfLJyL8kHTvDctfE0DA20Oo7ycX+9Jp7QoAC9XIaqVU4XgzGmBHgAWARsBeYYY7aIyAsicvooo9eA+sAXlQ5L7QIkichGYBkw2RjjNsVw6Hg+r/+4nUs6RjG8e3Or4yhFZP1gnhzWmdW7c5i7Nt3qOMpLBbjiRYwxC4GFlaY9W+7+FWdY7jeguysy1Ia/LUihpMzwYmI3vXKqchs32mKYty6dSQu3cnmXpkSEBVkdSXkZPUPrDJZuPcIPWw7z0BUdaN1Ex29W7sPPT3h5VHdOFpTw8ndbrY6jvJAWQxXyikp49ustdGxan3su1steKPfTsWkD7r2kLfPWpfPbziyr4ygvo8VQhTeWpHLgWD6TRnXXy14ot/XnyzrQOiKUZ+Yn67kNyqX0U6+SrYdO8MGK3YzpG4MtNsLqOEqdUUigPy+N7MaurFzeWa7nNijX0WIop6zM8NT8zTSqF8jEYZ2tjqNUtQZ1jGJEzxZMW7aTnZl6boNyDS2Gcmau2cf6fcd4engXGoXqkR7KMzxzTRdCAv303AblMloMDhknC3j1+20MbNeEUb3OdqknpdxLdIMQJg7rwqpdOcxbd8DqOMoLaDE4vPTtVgqKy3hxpJ6zoDzPmL4x9GnTmJe/06FAlfO0GICfd9gHXb/v0na006E6lQfy8xMmOc5teGWhntugnOPzxVBQXMpfv04mLjKM+y6t8sKuSnmETs0acM+gtnyxNp2VO7OtjqM8mM8Xw9RlaezNzuPlkd0ICdRxFpRne/CyDsRE1OPprzbruQ3qvPl0MaRlnOTd/+xkVK+WOs6C8gr1gvx5MbEbuzJzeXf5LqvjKA/ls8VgjOHp+cn2SxgP13EWlPe4tFM01/RoztTlaezOyrU6jvJAPlsMc9ems3p3DhOHdSayfrDVcZRyqWeviSfY349nvtJzG9S588liyMktYtLCrdjaNOYmW0z1CyjlYaLDQ3h8aCd+Tcvm6w1/GFBRqbNySTGIyFAR2S4iaSIysYrng0VktuP51SISW+65Jx3Tt4vIVa7IU51XFm61X7J4VHf8/PScBeWdxvVvQ8+YRrz0XQrH84qtjqM8iNPFICL+wFRgGBAPjBWR+Eqz3QUcNca0B/4JvOpYNh77UKBdgaHANMfr1ZrVu7L5Ym06d1/clk7NGtTmWyllKX8/YdKobhzNK2byD9usjqM8iCu+MfQD0owxu4wxRcAsILHSPInADMf9ucDlYj+9OBGYZYwpNMbsBtIcr1crikrKePqrZFo1rsdDl3eorbdRym10bdGQOwbGMvP3fazdm2N1HOUhXFEMLYH95R6nO6ZVOY9jjOjjQJMaLusy03/eSVrGKV5M7Ea9ID1nQfmGh4d0pEXDEJ76Mpni0jKr4ygP4IpiqGojfeXDIM40T02Wtb+AyAQRSRKRpMzMzHOMaHfoeAHDuzdncOfo81peKU8UFhzA3xK7sf3IST74ZbfVcZQHcEUxpAPlD+1pBVQ+DOK/84hIANAQyKnhsgAYY6YbY2zGGFtUVNR5BX15VHfeGJNwXssq5cmGxDflyvimvLF0B/tz8qyOo9ycK4phDdBBROJEJAj7zuQFleZZAIx33L8e+MnYD65eAIxxHLUUB3QAfndBpjMK0KE6lY96fkRX/EV49utkPbdBnZXTn5KOfQYPAIuArcAcY8wWEXlBREY4ZvsQaCIiacAjwETHsluAOUAK8APwJ2OMXuBFqVrQolE9Hh7SkWXbM/k++bDVcZQbE0/8zcFms5mkpCSrYyjlcUpKyxjx9q9k5xay5JFLaBASaHUkVYdEZK0xxlbdfLpdRSkfEuDvx6TR3ck4WcjrP+6wOo5yU1oMSvmYhJhG3DqgDTNW7mFT+jGr4yg3pMWglA/6y1WdiKofzFPzN1Oi5zaoSrQYlPJB4SGBPHdtV5IPnGDGyr1Wx1FuRotBKR91dfdmXNopiv/7cTuHjudbHUe5ES0GpXyUiPBiYjdKjeH5BVusjqPciBaDUj4sJiKUBy/vwKItR1iScsTqOMpNaDEo5ePuubgtHZvW57kFW8grKrE6jnIDWgxK+bhAfz8mjerOgWP5TFmSanUc5Qa0GJRS2GIjGNM3hg9X7Cbl4Amr4yiLaTEopQCYOKwzjeoF8vRXmykr87xL5SjX0WJQSgHQKDSIp4d3Yf2+Y/z7931Wx1EW0mJQSv3XqF4tGdiuCa/+sI2MkwVWx1EW0WJQSv2XiPDSyG4UFpfx4rdbrY6jLKLFoJSqoG1Ufe4f3I5vNh7k5x3nN4yu8mxaDEqpP7jv0na0jQzjma+SKSjWsbN8jVPFICIRIrJYRFIdPxtXMU+CiKwUkS0isklEbir33McisltENjhuOiCzUm4gOMCfl0Z1Y19OHm//lGZ1HFXHnP3GMBFYaozpACx1PK4sD7jNGNMVGApMEZFG5Z5/zBiT4LhtcDKPUspFBraLZHSvlrz3807SMk5aHUfVIWeLIRGY4bg/AxhZeQZjzA5jTKrj/kEgA4hy8n2VUnXgqeFdCA0K4Kn5yXjiMMDq/DhbDE2NMYcAHD+jzzaziPQDgoCd5Sa/7NjE9E8RCXYyj1LKhSLrB/PksM78vjuHL9amWx1H1ZFqi0FElohIchW3xHN5IxFpDnwK3GGMOT1k1JNAZ6AvEAE8cZblJ4hIkogkZWbqkRJK1ZUbbTH0jW3MpIVbyTpVaHUcVQeqLQZjzBXGmG5V3L4Gjjg+8E9/8GdU9RoiEg58BzxjjFlV7rUPGbtC4F9Av7PkmG6MsRljbFFRuiVKqbri5ye8Mro7eYWlOm6Dj3B2U9ICYLzj/njg68oziEgQMB/4xBjzRaXnTpeKYN8/kexkHqVULWgf3YA/X9aebzcd4scth62Oo2qZs8UwGRgiIqnAEMdjRMQmIh845rkRGATcXsVhqZ+LyGZgMxAJvORkHqVULfmfS9vRuVkDnvkqmeP5xVbHUbVIPPFIA5vNZpKSkqyOoZTP2Zx+nMSpK7jRFsPk63pYHUedIxFZa4yxVTefnvmslKqx7q0acs+gtsxas59f07KsjqNqiRaDUuqcPHxFR+Iiw5j45SYdCtRLaTEopc5JSKA/k0d3Z39OPq//uMPqOKoWaDEopc5Z/7ZNuGVAaz76dTfr9h21Oo5yMS0GpdR5eWJoZ5qHh/DE3E0UlugVWL2JFoNS6rw0CAnk5VHdSc04xdRlO6tfQHkMLQal1Hkb3DmaUb1aMm1ZGikHT1gdR7mIFoNSyinPXhNPo9AgHv1iI0UlZdUvoNyeFoNSyimNw4J4ZXR3th46wds/pVodR7mAFoNSymlD4ptyXe9WTF2+k437j1kdRzlJi0Ep5RLPXhtPVP1gHv1io44T7eG0GJRSLtGwXiB/v74HaRmn+L/FeuKbJ9NiUEq5zKCOUdzcvzXv/7KLNXtyrI6jzpMWg1LKpZ66ugutGtfjL19s1GspeSgtBqWUS4UFB/Da9T3Zl5PH5O+3WR1HnQctBqWUyw1o24Q7L4zjk5V7+SVVx2j3NE4Vg4hEiMhiEUl1/Gx8hvlKy43etqDc9DgRWe1YfrZjGFCllBd47KpOtI+uz6NzNpKTW2R1HHUOnP3GMBFYaozpACx1PK5KvjEmwXEbUW76q8A/HcsfBe5yMo9Syk2EBPrzxpgEjuUV88S8TXjiaJG+ytliSARmOO7PAEbWdEEREeAyYO75LK+Ucn9dWzTk8aGdWJxyhJm/77c6jqohZ4uhqTHmEIDjZ/QZ5gsRkSQRWSUipz/8mwDHjDGnD1tIB1qe6Y1EZILjNZIyM3WbpVKe4s4L47i4QyQvfLuFtIxTVsdRNVBtMYjIEhFJruKWeA7v09oxAPU4YIqItAOkivnO+F3TGDPdGGMzxtiioqLO4a2VUlby8xP+cUNP6gX687+z1+uF9jxAtcVgjLnCGNOtitvXwBERaQ7g+Jlxhtc46Pi5C1gO9AKygEYiEuCYrRVw0Ok1Ukq5nabhIbx6XQ+SD5zg9cXbrY6jquHspqQFwHjH/fHA15VnEJHGIhLsuB8JXAikGPueqGXA9WdbXinlHa7s2oxx/Vsz/edd/JaWZXUcdRbOFsNkYIiIpAJDHI8REZuIfOCYpwuQJCIbsRfBZGNMiuO5J4BHRCQN+z6HD53Mo5RyY38dHk9cZBgPz9lA1qlCq+OoMxBPPITMZrOZpKQkq2Mopc5DysETjJz2K/3jIphxRz/8/Kra3ahqg4isdezvPSs981kpVafiW4Tz/LVd+SU1i2nL06yOo6oFGzQMAAAMgUlEQVSgxaCUqnNj+8WQmNCC/1u8g5U7s62OoyrRYlBK1TkRYdKo7sRGhvHgrPVkntT9De5Ei0EpZYmw4ACmjuvNifxiHp69gdIyz9vf6a20GJRSlunSPJwXEruyIi2Lt35KtTqOctBiUEpZ6kZbDKN7t+SNpan8tO2I1XEUWgxKKYud3t8Q3zych2ZuYFemXk/JaloMSinLhQT6896tfQgM8OPeT9dyqlCHBLWSFoNSyi20ahzK2+N6sSsrl0fnbKBMd0ZbRotBKeU2BraL5MlhnVm05Yie/GYhLQallFu566I4EhNa8PriHSxJ0Z3RVtBiUEq5FRFh8ugedG/ZkAdnrSf5wHGrI/kcLQallNupF+TPB7fZaFgvkLtnJHH4eIHVkXyKFoNSyi1Fh4fw4fi+nCwo5q4Za8jVI5XqjBaDUsptxbcI561xvdh66AQPzdLLZtQVLQallFu7rHNTnr0mniVbj/Dityl44hgynsapYhCRCBFZLCKpjp+Nq5hnsIhsKHcrEJGRjuc+FpHd5Z5LcCaPUso73X5hHHdeGMfHv+1h2vKdVsfxes5+Y5gILDXGdACWOh5XYIxZZoxJMMYkAJcBecCP5WZ57PTzxpgNTuZRSnmpZ4Z3YWRCC15btJ2Zv++zOo5Xc7YYEoEZjvszgJHVzH898L0xJs/J91VK+Rg/P+G1G3pyaaconp6/mR+SD1kdyWs5WwxNjTGHABw/o6uZfwwws9K0l0Vkk4j8U0SCz7SgiEwQkSQRScrMzHQutVLKIwX6+zHt5t4kxDTiwZkb+G1nltWRvFK1xSAiS0QkuYpb4rm8kYg0B7oDi8pNfhLoDPQFIoAnzrS8MWa6McZmjLFFRUWdy1srpbxIaFAAH93elzZNQrlnRhJJe3KsjuR1qi0GY8wVxphuVdy+Bo44PvBPf/BnnOWlbgTmG2OKy732IWNXCPwL6Ofc6iilfEGj0CA+v7s/TcNDGP/R76zdq+XgSs5uSloAjHfcHw98fZZ5x1JpM1K5UhHs+yeSncyjlPIR0eEhzJwwgOjwEMZ/tIZ1+45aHclrOFsMk4EhIpIKDHE8RkRsIvLB6ZlEJBaIAf5TafnPRWQzsBmIBF5yMo9Syoc0DQ9h5j0DiKwfxPgPf2e9loNLiCeeLGKz2UxSUpLVMZRSbuLQ8XzGTF9F9qki3r/NxgXtmlgdyS2JyFpjjK26+fTMZ6WUx2vesB6zJ1xA84YhjP/X7yzWy3U7RYtBKeUVmjUMYc69F9CleTj/89la5q1NtzqSx9JiUEp5jcZh9qOVBrSN4NEvNvLuf3bqtZXOgxaDUsqr1A+2n+dwTY/mTP5+G0/M20RRSZnVsTxKgNUBlFLK1YID/HlzTC/aRobx5k9p7MvJ491b+tAoNMjqaB5BvzEopbySn5/wyJWdmHJTAuv2HmPUtN/Yfvik1bE8ghaDUsqrjezVkn/f059ThSUkTl3Bl+t0p3R1tBiUUl7PFhvBdw9eREJMIx6Zs5Env9xEQXGp1bHclhaDUsonRDcI4bO7+nPfpe2Y+ft+rn1rBZvTj1sdyy1pMSilfEaAvx9PDO3MjDv7caKgmFHTfmXKkh0Ul+pRS+VpMSilfM4lHaP48X8v4ZoezZmyJJWRU3/Vi/CVo8WglPJJDUMDmTKmF+/c3JusU4WMnvYbE+dtIie3yOpoltNiUEr5tGHdm7P00Uu55+I4vlibzuB/LOe9/+wkv8h3d05rMSilfF794ACeHh7PwgcvpmdMI175fhuXvLaMT1fucaujl+rqDG4tBqWUcujUrAGf3NmP2RMG0KZJKH/9egsXTv6JKUt2kHWq0LJcOzNP8dzXyfSftISMkwW1/n5OXRJDRG4Ange6AP2MMVUOkiAiQ4E3AH/gA2PM6QF94oBZ2Md7XgfcaozRDXxKKUv1b9uEOfdewMpd2Xz4y26mLEll2vKdDOnSlNG9WzKoYxSB/rX7e/WpwhJ+3HKYL9cdYEVaFkH+flzTo3mdfGtwaqAeEekClAHvAX+pqhhExB/YgX2Et3RgDTDWGJMiInOAL40xs0TkXWCjMead6t5XB+pRStWltIxTfLZqLws2HiQnt4gmYUEM7hzN4E7RXNQhkob1Al3yPgeO5bMiNZPl2zP5aVsGhSVltGxUjzF9YxjbvzWR9YOdev2aDtTjkhHcRGQ5Zy6GC4DnjTFXOR4/6XhqMpAJNDPGlFSe72y0GJRSVigqKeM/OzJZsPEgP+/I5Hh+MX4CHZs2oGerRvSIaUhcZBgxjUNp1jDkjN8q8opKOHy8gPSj+Ww9dIItB0+w+cBxdmflAhDdIJhh3ZoxIqElvVs3QkRckr+mxVAXV1dtCewv9zgd6A80AY4ZY0rKTW9ZB3mUUuq8BAX4MSS+KUPim1JSWsbG9GP8Z0cWG/cfY1HKYWYn/f+POj+BsOAAQoP8CQ0KoLTMUFhSSl5RKScLSiq8bouGIcS3aMjN/VtzcYcoOjat77IyOB/VFoOILAGaVfHU08aYr2vwHlWtnTnL9DPlmABMAGjdunUN3lYppWpPgL8ffdpE0KdNBADGGNKP5rM/J4/0o/mkH83jREEJ+UWl5BaVEOAnBAf4Uy/In+jwYJo3DKFZeD06NWtARJh7XQ682mIwxlzh5HukAzHlHrcCDgJZQCMRCXB8azg9/Uw5pgPTwb4pyclMSinlUiJCTEQoMRGhVkdxWl0crroG6CAicSISBIwBFhj7zo1lwPWO+cYDNfkGopRSqhY5VQwiMkpE0oELgO9EZJFjegsRWQjg+DbwALAI2ArMMcZscbzEE8AjIpKGfZ/Dh87kUUop5TyXHJVU1/SoJKWUOnc1PSpJz3xWSilVgRaDUkqpCrQYlFJKVaDFoJRSqgItBqWUUhV45FFJIpIJ7D3PxSOxn1znS3SdfYOus/dzdn3bGGOiqpvJI4vBGSKSVJPDtbyJrrNv0HX2fnW1vropSSmlVAVaDEoppSrwxWKYbnUAC+g6+wZdZ+9XJ+vrc/sYlFJKnZ0vfmNQSil1Fl5bDCIyVES2i0iaiEys4vlgEZnteH61iMTWfUrXqsE6PyIiKSKySUSWikgbK3K6UnXrXG6+60XEiIhHH8FSk/UVkRsdf89bROTfdZ3R1Wrw77q1iCwTkfWOf9tXW5HTlUTkIxHJEJHkMzwvIvKm489kk4j0dmkAY4zX3QB/YCfQFggCNgLxlea5H3jXcX8MMNvq3HWwzoOBUMf9+3xhnR3zNQB+BlYBNqtz1/LfcQdgPdDY8Tja6tx1sM7Tgfsc9+OBPVbndsF6DwJ6A8lneP5q4HvsI2EOAFa78v299RtDPyDNGLPLGFMEzAISK82TCMxw3J8LXC5WDrLqvGrX2RizzBiT53i4CvuoeZ6sJn/PAC8CfwcK6jJcLajJ+t4DTDXGHAUwxmTUcUZXq8k6GyDccb8hZxkJ0lMYY34Gcs4ySyLwibFbhX00zOauen9vLYaWwP5yj9Md06qcx9gHEzqOfbAgT1WTdS7vLuy/cXiyatdZRHoBMcaYb+syWC2pyd9xR6CjiPwqIqtEZGidpasdNVnn54FbHIOGLQT+XDfRLHWu/9/PSbVjPnuoqn7zr3z4VU3m8SQ1Xh8RuQWwAZfUaqLad9Z1FhE/4J/A7XUVqJbV5O84APvmpEuxfyP8RUS6GWOO1XK22lKTdR4LfGyMeV1ELgA+daxzWe3Hs0ytfn556zeGdCCm3ONW/PHr5X/nEZEA7F9Bz/bVzd3VZJ0RkSuAp4ERxpjCOspWW6pb5wZAN2C5iOzBvi12gQfvgK7pv+uvjTHFxpjdwHbsReGparLOdwFzAIwxK4EQ7NcU8mY1+v9+vry1GNYAHUQkTkSCsO9cXlBpngXAeMf964GfjGOvjoeqdp0dm1Xew14Knr7tGapZZ2PMcWNMpDEm1hgTi32/yghjjKeOC1uTf9dfYT/IABGJxL5paVedpnStmqzzPuByABHpgr0YMus0Zd1bANzmODppAHDcGHPIVS/ulZuSjDElIvIAsAj7UQ0fGWO2iMgLQJIxZgHwIfavnGnYvymMsS6x82q4zq8B9YEvHPvZ9xljRlgW2kk1XGevUcP1XQRcKSIpQCnwmDEm27rUzqnhOj8KvC8iD2PfnHK7h/+Sh4jMxL45MNKx7+Q5IBDAGPMu9n0pVwNpQB5wh0vf38P//JRSSrmYt25KUkopdZ60GJRSSlWgxaCUUqoCLQallFIVaDEopZSqQItBKaVUBVoMSimlKtBiUEopVcH/A/hIMWCeYyJHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_01 = np.linspace(0, 1, 1000)\n",
    "pt.plot(x_01, f(x_01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.  , 0.25, 0.5 , 0.75, 1.  ])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "degree = 4\n",
    "h = 1\n",
    "\n",
    "nodes = 0.5 + np.linspace(-h/2, h/2, degree+1)\n",
    "nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the gen. Vandermonde matrix and find the coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "V = np.array([\n",
    "    nodes**i\n",
    "    for i in range(degree+1)\n",
    "]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "coeffs = la.solve(V, f(nodes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the interpolant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x_0h = 0.5+np.linspace(-h/2, h/2, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "interp_0h = 0*x_0h\n",
    "for i in range(degree+1):\n",
    "    interp_0h += coeffs[i] * x_0h**i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5f98d180f0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4FNX+x/H32fRACAFCDSEUQ++hIwRDh4QuIFKFCF5ErheEK/4UEa4V2xVEVMpVlKJSRRAwqDQlCCK9SQmdUBJISD2/P2YvlxIkkN2d3ez39Tx5sjs72fnMJvnu7Jkz5yitNUIIIdyLxewAQgghHE+KvxBCuCEp/kII4Yak+AshhBuS4i+EEG5Iir8QQrghKf5CCOGGpPgLIYQbkuIvhBBuyNPsAHdTrFgxHRYWZnYMIYRwKdu2bbugtQ6+13pOW/zDwsKIj483O4YQQrgUpdSx3KwnzT5CCOGGpPgLIYQbkuIvhBBuyGnb/IUQrikjI4OEhASuX79udpR8zdfXl5CQELy8vB7o56X4CyFsKiEhgYCAAMLCwlBKmR0nX9Jak5iYSEJCAuXLl3+g57BJs49SapZS6pxSatddHldKqfeVUoeUUjuVUvVssV0hhPO5fv06RYsWlcJvR0opihYtmqdPV7Zq858DtP+LxzsAD1m/YoEPbbRdIYQTksJvf3l9jW3S7KO1/kkpFfYXq3QB/qONOSO3KKUKK6VKaa1P22L7wjyXLl3i3LlzJCcnk5GRQWZmJllZWbRo0QKLxcLJkydJS0ujcOHCBAYG4uHhYXZkIQSOa/MvA5y46X6CddktxV8pFYvxyYDQ0FAHRRO5kZGRwfHjxzl69Chnzpyhd+/eeHp68uuvv7Jly5Zb1lVK0bJlSwDi4+PZsWPHjeXFihUjJCSE6OhoOToUwkSOKv45/ZffMXO81nomMBMgIiJCZpZ3AseOHePHdeu4/ttvBJ86RdFLl2iQkgIffABJSbRJTKRNcvL/fsEeHhAUhFq8GIoWpUPZsjxcogSXg4M5U6oUR7XmypUrNwr/0qVLUUpRtWpVKlSoIJ8MhHAQRxX/BKDsTfdDgFMO2rb4K/PmwYQJcPw4hIaS/cor7K5enVLHj1MsPp6S339P399/xys9HQBtsaDKloWwMKhaFUtQEBQqBBbr6aPMTLh8GRIT4dw5vFeupMi5cxQBKgBNixWDBg3g3Dl027ZkZ2Wx/8ABtm/fjr+/PzVq1KBevXqUKFHCrFdE5CP/+Mc/+O6774iKiuLf//632XGciqOK/zJgpFJqPtAIuCLt/U5g3jyIjYWUFOP+sWMwcCCVLRa8s7LAwwPvBg1QI0ZAw4ZQty6qQgXw8bm/7SQlwf79sG0b/PorbN4MY8aggG5lypAdHU1Cixb8YrGwbds2/P39KVGiBNnZ2QBYLHItorh/R44cYePGjezZs8fsKE7JJsVfKfUlEAkUU0olAC8BXgBa6xnASqAjcAhIAQbbYrsijyZM+F/ht7JojaePD3ruXFTr1qjChfO+nUKFjKP9Bg1g+HBj2YkTsHo1fPcdlrlzCZ0xg9By5ch49FGySpUCYN++ffzwww+0bNmSGjVqyDkCkWv79++ndevWZGZmUrduXTZs2ECBAgXMjuVUlNEBx/lERERoGdXTTrQmY+lSvLp1y/lxpcB61O0QycmwZInxSWTNGmNZly6c7tmTJZcvc+78eYKDg3nkkUeoXLmyvAk4ub1791K1atUb9+fMmXPHOtWrV6dBgwZkZGQwb968Ox6vU6cOderUISUlhYULF97y2KBBg3KV44UXXiAsLIyhQ4feV35XcvtrDaCU2qa1jrjXz8rnaXeiNdlLlpBStSpe3bqh79ac4uieVgEB0L8/rFoFR47A2LHw44+U6teP4XPnMrB4cbKzsliwYAHLli1zbDbhsv744w9q165tdgznpbV2yq/69etrYUNr1ui0WrW0Bp0YFKR/HDJEX3rvPa39/bWG/335+2v9+edmp9U6JUXrGTO0LldOa9DZDRrove+9p3fv3q211jorK0unp6ebm1HkaM+ePWZH0FprXbFiRX3t2jWzY9hVTq81EK9zUWPlyD+/O3IEunaFNm24dvw4q3r14uTatTz8yScUHjUKZs6EcuWMpp5y5Yz7/fqZnRr8/ODJJ+HAAfj4Y9T581R55hmqjRkDe/eyadMmZsyYwfHjx81OKpxQcnIyXl5e+Pv7mx3FaUnxz6+uX4f/+z+oVg3WruXUyJH8/NFHtJw7l5r16v2v3bxfPzh61GjjP3rUOQr/zby9YehQ2LcP3noLNm2CmjWp9fHHeF67xuzZs1m1ahUZGRlmJxVOZNeuXdSoUcPsGE5NTvjmR5s3o4cMQe3bx8UOHSjy8cdQpozZqWzj/Hl46SX46CN0qVJsGzyYbz09KV68OD179iQ4+J5Tlwo7y+kkpLAPOeErDKmpMGYMulkzUs6f5/PHH+eHIUPQpUubncx2goNh+nTYsgVVtCgRkyfzjy1bUGfPyvjxQtwHGc8/v9i3D3r3hp07+aNJE75r1YqW0dE0atQof3aNbNAA4uPhzTcpOGkST27diqpbF8qWZd++fVSqVAlPT/nzFuJu5L8jP5g7F556imw/PxYNGMCpOnV4/NFHKZNfmnruxssLnn8eundHPfYYdO1K6oABfB0SQskKFejduzcFCxY0O6UQTkmKvytLTTWumP3PfyAyEsu8eTx05gydwsPdq+hVqWIMGfF//4ffm28ypnx5/tO5Mx8nJdGnTx9KWa8YFkL8j7T5u6qEBHj4YfRnn7Hn0Uc5MXs2lC5NvXr13Kvw/5ePD7zxBqxdi8/16wz9+GMqb9vGrFmz2Lt3r9nphHA6Uvxd0ebNEBGB3r+fdU8/zaJq1Th74YLZqZxDVBT89huqQQM6zptH17g4UpKSzE4lhNOR4u9qPv8cIiPJ9vdnwejRbC5WjB49ehARcc+eXe6jZElYtw6eeYbq69ZRf+xYOH+eCxcu4Kxdm4Vt5ebT77vvvkvKbQMb2sOcOXMYOXLkA/3s5cuXmT59uo0TGaT4uwqtjWaN/v3JbNSIj4YM4U9/fx5//HG5mCUnXl7w7rvGYHG//kpWo0YsmjyZFStW3BgqWri3Byn+WVlZdkqTMyn+7i47G559FsaNgz59sKxeTbl69Rg0aBDly5c3O51ze+wxWL8ey9WrDJ01i8Svv2bp0qXyBuAm1q9fT2RkJD179qRKlSr069cPrTXvv/8+p06dolWrVrRq1QqA77//niZNmlCvXj169erF1atXAQgLC2PSpEk0b96cRYsWERkZyejRo2natCk1atTg119/BeDixYt07dqVWrVq0bhxY3bu3HlHnuXLl9OoUSPq1q1L69atOXv2LAATJ05kyJAhREZGUqFCBd5//30Axo8fz+HDh6lTpw5jx4616WsjvX2cXXo6DBwI8+dzLTaW7DfeIMDPj44dO5qdzHU0boz65Re8OnZkwLx5LL1yhcXZ2XTr1k0mirG30aPBOoezzdSpY3yqy6Xt27eze/duSpcuTbNmzdi4cSOjRo3i7bffJi4ujmLFinHhwgUmT57M2rVrKVCgAK+//jpvv/02L774IgC+vr5s2LABgBkzZnDt2jU2bdrETz/9xJAhQ9i1axcvvfQSdevWZcmSJfzwww8MGDDgxvzV/9W8eXO2bNmCUopPPvmEN954g6lTpwLG/BVxcXEkJydTuXJlRowYwWuvvcauXbvueB5bkOLvzNLSoGdPWLGCpAkTmBEQQKlly+jfv7/ZyVxP+fKweTOWnj3ptmQJ36ek8Fu5cnKuxA00bNiQkJAQwJgn4OjRozRv3vyWdbZs2cKePXto1qwZAOnp6TRp0uTG4717975l/b59+wLQokULkpKSuHz5Mhs2bODrr78G4JFHHiExMZErV67c8nMJCQn07t2b06dPk56efssn906dOuHj44OPjw/Fixe/8anAXqT4O6vr16FbN1i1isuvvspHSuHj40Pnzp3NTua6CheGlSuhf3/aLlyIrl8f6tc3RjQV9nEfR+j24nPTtKMeHh5kZmbesY7WmjZt2vDll1/m+By3zwJ2+1XzSqkcOxPcvt7TTz/Ns88+S0xMDOvXr2fixIn3ldOW5DOvM0pJgZgYWL2a5LffZqbFgo+PD4MGDSIoKMjsdK7N2xu++AJiY1GvvkrasGGsWb1azgG4oYCAAJKTkwFo3LgxGzdu5NChQwCkpKRw4MCBu/7sggULANiwYQOBgYEEBgbSokWLG7OSrV+/nmLFilGoUKFbfu7KlSs3rryfO3fufWW0NTnydzYpKRAdDXFxMHs2SywWPM6dY8CAARS2xXy6Ajw8YMYMCArC5/XXKfnLL3w7aRKdu3bNn+MgiRzFxsbSoUMHSpUqRVxcHHPmzKFv376kpaUBMHnyZMLDw3P82aCgIJo2bUpSUhKzZs0CjJO2gwcPplatWvj7++dY3CdOnEivXr0oU6YMjRs35s8///zLjEWLFqVZs2bUqFGDDh068Oabb+Zxr2+SmxlfzPhyy5m80tK0btdOa4tF688+01prfe3aNX3+/HmTg+Vjr76qNeg/qlfXq7/9VmdnZ5udyOU5y0xe9tKyZUu9detWs2NorfM2k5cc+TuLrCxjIpXVq0mbNo0fixcnKisLf39/mY3InsaPR1ss1Bg3Dj1uHD/5+NAyKsrsVELYnbT5OwOtITYWvvqKzNdf5z9eXmzdupVz586ZncwtqOeeQ7/2GjV37aLM+PFkpKaaHUk4sfXr1+eLXmJS/M2mNfzjHzBrFvqFF1hUtiynT5+mZ8+eMhqlA6lx48h+9VUqxcfj9cQTYOeeFvmdlmE07C6vr7EUf7O99Ra88w766adZ2agRBw4coEOHDlSuXNnsZG7HMn48vPoqfPklf7Zrx8mEBLMjuSRfX18SExPlDcCOtNYkJibi6+v7wM8hbf5mWrgQnnsOHn2UpJdfZtfMmTRt2pQGDRqYncx9jR9P+uXLlH/9dbb264f/4sUEFSlidiqXEhISQkJCAufPnzc7Sr7m6+t74+K1ByETuJtl40Zj+OEGDWDNGvD15fLlywQGBkp3Q7NpTerQofjNmsWmzp2ps2CBnHQXLkMmcHdmBw4YF3GVK8fJadPYuG0bWmsKFy4shd8ZKIXfxx9zrUsXmq5Ywe8jRjh8NEch7E2Kv6NduAAdO4KHB8kLFvDl99/z22+/kZGRYXYycTOLhQKLFpH08MM0/uwzrs2ebXYiIWxKir8jZWTAo49CQgKZ33zDF7/8QkZGBn369MHb29vsdOJ2Xl4UWr0a3bgxhf72N7CO6ihEfiDF35HGjIG4OPRHH7HkzBnOnDlDjx49CA4ONjuZuBs/PyzLl6PDwsjo1IkTa9eanUgIm5Di7yizZsH778Pf/87JqCj27NlD69at7zp2iHAiRYuSsWQJGdnZBHXtSlaZMmCxQFiYMVOYEC5Ievs4wpYt0LIltGgB330Hnp6cOXOGEiVKyAleF5Iyfjx+r7/OLb8xf3+YOdMYmkMIJyC9fZzFqVPQvTuEhJA4fTqHjh4FoGTJklL4XYz//Pnc8RtLSYEJE8yII0SeSPG3p4wM6N0bkpLI+OorFqxZw5IlS0hPTzc7mXgQx4/f33IhnJgUf3uaMAE2bEDPnMmyI0e4cOEC3bt3l549rio09P6WC+HEpPjby7Jl8OabMGIEWytVYteuXURGRlKhQgWzk4kHNWWK0cZ/Ew2kPf64OXmEyAMp/vbw558wcCDUq0fihAmsXr2ahx56iIcfftjsZCIv+vUzTu6WKwdKkV26NCn+/mROm0aWDAInXIwUf1tLSzMu5NIaFi2iSOnSdOrUiW7duskJ3vygXz84ehSys7GcPMnJuXPxunaN5Nat4fp1s9MJkWs2Kf5KqfZKqf1KqUNKqfE5PD5IKXVeKbXD+jXUFtt1SmPGQHw82bNmcblIEZRS1KtXDz8/P7OTCTsI79mTP557jsL795PUp4/xpi+EC8hz8VdKeQDTgA5ANaCvUqpaDqsu0FrXsX59ktftOqUVK+CDD2D0aDYGB/Phhx9y6dIls1MJO6s9cSJbo6MptHQp6e++a3YcIXLFFkf+DYFDWusjWut0YD7QxQbP61rOnIEhQ6B2bU787W/ExcURHh5O4cKFzU4m7MzT05MKn37KlRYt8B43zrioTwgnZ4viXwY4cdP9BOuy2/VQSu1USn2llCprg+06j+xsGDwYkpNJmz2br1esIDAwkE6dOkk7v5soGhxM4JIlULYsukcPkPmXhZOzRfHPqbrd3vC5HAjTWtcC1gJzc3wipWKVUvFKqXiXmgXo3/+GVavQU6ey7NAhkpOT6dGjR56mWBMuKCiIU//+N5nnzpHatavMAyycmi2KfwJw85F8CHDq5hW01ola6zTr3Y+B+jk9kdZ6ptY6Qmsd4TIjXe7caUzFGB1Ndmwsvr6+tGrVKk/TqwnXVax1a37s3Ru/zZtJH39H3wchnEaeB3ZTSnkCB4Ao4CSwFXhMa737pnVKaa1PW293A8ZprRv/1fO6xMBuqanGNIyJicabgPUNS2stzT1u7OzZs5zs3Jl68fHob75BdetmdiThRhw2sJvWOhMYCawG9gILtda7lVKTlFIx1tVGKaV2K6V+B0YBg/K6XafwwguwezeZn3zCwrg4zpw5AyCF382VKFGC7Hfe4WTp0mT17w+HD5sdSYg72KSfv9Z6pdY6XGtdUWs9xbrsRa31Muvtf2qtq2uta2utW2mt99liu6basAHeeQdGjOB7pdi7dy9Xr141O5VwEvWbNWPrc8+hlYI+fUAG8xNORq7wfRDXrhm9e8LCOBQby9atW2ncuDGVKlUyO5lwEkopYp5+Gq/PPoP4eJD2f+FkpPg/iOefh0OHuP7hhyxdt47g4GCioqLMTiWcjMViga5dSR40yPiUuHy52ZGEuEGK//1av96YjnHUKDZ6epKSkkL37t3x9PQ0O5lwUj/HxHC6ZEmyBgyAEyfu/QNCOIBM43g/rl6FmjXB0xN27CDL15eEhATKlStndjLhxNLT05n/yiv0efNNPOvXx/Ljj8bfkBB2INM42sNzz8GxY6RMn06qxYKHh4cUfnFP3t7eRA0fzorOnbFs2gQTJ5odSQgp/rm2fj18+CF69GgWnTrFp59+SnZ2ttmphIsoU6YMxUaNYnvduuh//QvWrjU7knBzUvxzIzUVYmOhYkV+jYnh6NGjNGvWzDihJ0QuNW/eHMsHH0CVKsZkP4mJZkcSbkyqV25MngwHD3LptddYs2EDVapUoU6dOmanEi7GYrFQu2lT1BdfoM+fR8fGyvj/wjRS/O/ljz/gjTfIHjCABRcu4OvrS+fOneUqXvHALoaG8mObNqhvvoE5c8yOI9yUFP+/kpUFw4ZBUBBpU6ZQqFAhoqOjKVCggNnJhAsLCgrieM+eHCtfHj1qlAz/IEwhxf+vTJ8Ov/wC776LX0gIffv2pXLlymanEi5OKUV0164s79WL9Oxs9OOPy/DPwuGk+N/NiRPw/PPotm1ZERDApUuXpKlH2ExQUBCNe/VieceOqC1bjPNKQjiQFP+caA1/+xtkZ7Nl0CC2/fYbLjW5jHAJ9evXJzUmhuMtW8Irr8DmzWZHEm5Ein9OvvkGli8naexY1hw8SO3atQkPDzc7lchnlFL07t2bstbpH3n8cUhONjuWcBNS/G939SqMHo2uXZsvg4MpUKAA7dq1MzuVyKe8vb1RhQtzZdo09NGjMHq02ZGEm5Dif7vJkyEhgT+GD+fMhQt07twZPz8/s1OJfO77lBQ2NW8Os2bBihVmxxFuQIr/zfbuhalTYfBgKg8eTJcuXaR3j3CI9u3bs7ltWy6GhKCHDZOrf4XdSfH/L61h5Eh0wYJkTp6Mj4+PXMUrHCYgIICojh1Z2KkT+vx5ePppsyOJfE6K/38tWAA//MCBgQOZuWQJaWlpZicSbqZOnToUbN6cnyMj4csv4auvzI4k8jEZVByMHhbPPkt6rVosCgqieqlS+Pj4mJ1KuBmlFJ07d2ZrkSJkX7qEZcQIaNECihc3O5rIh+TIH2DiRPSZM6zo0AG/ggVp37692YmEmypcuDBtOnbE8p//oJOT4cknZfA3YRdS/Hftgvfe42znzvzh50fHjh2ld48w3ekiRdgaHQ1LlsC8eWbHEfmQexd/rWHUKHRgIOtat6Zq1apUrVrV7FRC4OXlxZqaNbkQHo5++mk4edLsSCKfce/iv2QJxMWhJk2i78iRxMTEmJ1ICACKFStGy0ce4Yt27dBpaTB0qDT/CJty3+KflgZjxpAeHk7qgAFYLBZ8fX3NTiXEDU2bNsW3Rg3i2rWDVauMC8CEsBH3Lf7vvgtHjvBN8+asWLXK7DRC3MFisRATE8OmOnW4VKcOPPusNP8Im3HP4n/mDHryZE7Wr8/hChWIiooyO5EQOSpZsiSDhgwhcOFCyMiA4cOl+UfYhHsW/wkT4Pp1vmnalMjISIoUKWJ2IiHuqmzZslgeeoi0l14yxv358kuzI4l8wP2K/7Zt6NmziW/WDO/q1WnSpInZiYS4p7S0ND5QikuVK8OoUXDunNmRhItzr+KvNTzzDLpoUU4MHEhMTAwWi3u9BMI1+fj4UKN2bb6IijIu/pKxf0QeuVflW7gQNm7E8q9/0X3wYEqVKmV2IiFyrVWrVmSGh7OlbVvjb3nxYrMjCRfmPsU/NRX93HNcLl+ei127mp1GiPvm7e1NdHQ0a+vWJbliRXjqKbh40exYwkW5T/F//33U8eMsbdGCK1evmp1GiAdSoUIFatWvz5Ynn0RfuGB0/xTiAbjHqJ4XLpA9ZQqHwsMJ6t6d8uXLm51IiAfWuXNnPDw84MoVmDIFeveGDh3MjiVcjFsc+We/8gpcvcrGmBjatGljdhwh8sTDwwOAxKee4nqFCsbIn0lJJqcSrib/F//Dh2H6dLbXrUvjIUNkxE6Rb8Rt2sSXbdqgT56EcePMjiNcTP4v/v/8J8rbG69//UtG7BT5Srt27TgbFsbuNm1gxgxYv97sSMKF5Ovir7dsgUWLUGPGUKtdO7PjCGFTAQEBtG3blqX165NWpgzExkJqqtmxhIvIv8Vfa66NGEFKQABXR4wwO40QdlG3bl1CwsP5pn17OHgQJk82O5JwETYp/kqp9kqp/UqpQ0qp8Tk87qOUWmB9/BelVJgttvtXUufPp+COHfzevTsFSpSw9+aEMIVSiujoaIo++ihZ/fvDG2/Azp1mxxIuIM/FXynlAUwDOgDVgL5KqWq3rfYEcElrXQl4B3g9r9u9q3nzoFw5fB97jCyLhRoNG6KUstvmhDBbkSJFaNu2LR7vvANBQTBsGGRlmR1LODlbHPk3BA5prY9ordOB+UCX29bpAsy13v4KiFL2qMjz5hntnsePowCP7GwCxo6VOVCFWzh5/Trru3WDX3+FDz4wO45wcrYo/mWAEzfdT7Auy3EdrXUmcAUoaoNt32rCBEhJuXVZSoqxXIh8ztPTk5/LlOFMvXrG3/yxY2ZHEk7MFsU/pyP422ebyM06KKVilVLxSqn48+fP33+S48fvb7kQ+UiJEiVo1rw581u0IFtrY+wfmfhF3IUtin8CUPam+yHAqbuto5TyBAKBO0ak0lrP1FpHaK0jgoOD7z9JaOj9LRcin2nRogVelSrxU9u2sHIlzJ9vdiThpGxR/LcCDymlyiulvIE+wLLb1lkGDLTe7gn8oLUdDkmmTAF//1uX+fsby4VwA56enkRHR/NTrVokV6sGzzwDiYlmxxJOKM/F39qGPxJYDewFFmqtdyulJimlYqyrfQoUVUodAp4F7ugOahP9+sHMmVCuHChlfJ8501guhJsIDQ2l/6BBFJg3Dy5dgjFjzI4knJCyxwG4LUREROj4+HizYwjh0tLHjsX7rbdgzRpo3drsOMIBlFLbtNYR91ov/17hK4Sbu3btGtOCgkgJCTFG/ry9J5xwa1L8hcinChQoQLkqVfiqTRs4cgReftnsSMKJSPEXIh9r3749Z6tWZX/z5uipU2H7drMjCSchxV+IfMzf35/27duzpFkzMgIDYehQyMw0O5ZwAlL8hcjnatSoQdlatdgxZAj89hu8957ZkYQTkN4+QriBzMxMPD08oEsXWLsWdu2CChXMjiXsQHr7CCFu8PT0BKW4OHkyWRYLDB8uQz+4OSn+QriRb3//nXWtWxv9/j//3Ow4wkRS/IVwIx06dODXevW4EB4Of/87PMgAiiJfkOIvhBspVqwYLVu1YmHr1uikJOMNQLglKf5CuJmmTZtiqVmTLZGRxkRHq1aZHUmYwNPsAEIIx/Lw8CAmJoZdZcqgjx9HDR9u9P4pWNDsaMKBpPgL4YZKly5N6dKljTl/H34YXnwR3n7b7FjCgaTZRwg3dqJcOQ5ERaHfew+2bjU7jnAgKf5CuLHMzEy+adCAtKAgY+iHjAyzIwkHkeIvhBsrX7481Zs2ZUnbtrBzJ7z1ltmRhINI8RfCzbVp04aT9etzuG5d9Msvw8GDZkcSDiDFXwg35+vrS6dOnVjyyCNkeXlBbKwM/eAGpPgLIahSpQrRsbFY3noL1q+HWbPMjiTsTIq/EAKA8PBwLMOGkdW8OXrMGDhzxuxIwo6k+Ashbrh05QqzGzdGp6TAqFFmxxF2JMVfCHFD4cKF8alVi59btoRFi2DZMrMjCTuR4i+EuEEpRefOndncvDmXypZFP/UUJCWZHUvYgRR/IcQtgoKCiGzblq/btYNTp+D5582OJOxAir8Q4g4NGzaERo34Mzoapk+HTZvMjiRsTObwFULkKC0tDZ+MDKheHQoUgO3bwcfH7FjiHmQOXyFEnvj4+EDBglx5/XXYuxdee83sSMKGpPgLIe5Ka83C5GT21q2LnjIF9uwxO5KwESn+Qoi7UkoRHR3NyjZtSPfxQQ8bBtnZZscSNiDFXwjxl0qWLEmj6Gi+a90atWkTfPSR2ZGEDUjxF0LcU9OmTbnYuTN/VqqEHjcOTp40O5LIIyn+Qoh7slgsdO3WjZMvvACZmTBihIz86eKk+AshcqVIkSI0HzgQ9corsHw5fPml2ZFEHkj5cAGvAAAQJElEQVTxF0LclzN9+nCmfHmyR46UkT9dmBR/IcR9KVCoECu6dyc7ORktzT8uS4q/EOK+BAQE0GTwYOIiI1FLlsDChWZHEg9Air8Q4r5Vr16dq7GxnCxThqwRI+DcObMjifskxV8I8UA6REfzQ//+kJwMI0eaHUfcJyn+QogH4uvrS/S4caiXXzYmfvnqK7MjifuQp+KvlCqilFqjlDpo/R50l/WylFI7rF8yNZAQ+UThwoWxPPccWXXrkvXkk3DhgtmRRC7l9ch/PLBOa/0QsM56PyepWus61q+YPG5TCOFMPD1Z268fXLlCWmys2WlELuW1+HcB5lpvzwW65vH5hBAuqMmwYWyKisJn8WKyvv7a7DgiF/Ja/EtorU8DWL8Xv8t6vkqpeKXUFqWUvEEIkc8UKlSI4Lfe4nTJkmQOGwaJiWZHEvdwz+KvlFqrlNqVw1eX+9hOqHVmmceAd5VSFe+yrVjrm0T8+fPn7+PphRBmq1KzJgeefx7PK1e4OmCAXPzl5O5Z/LXWrbXWNXL4WgqcVUqVArB+z7Gzr9b6lPX7EWA9UPcu683UWkdorSOCg4MfcJeEEGZpOnw4h/v3p+DKlTL2j5PLa7PPMmCg9fZAYOntKyilgpRSPtbbxYBmgEwHJEQ+5OXlRfgnn0CTJui//Y2sY8fMjiTuIq/F/zWgjVLqINDGeh+lVIRS6hPrOlWBeKXU70Ac8JrWWoq/EPmVpycZn35KZmoql7t1k5m/nFSeir/WOlFrHaW1fsj6/aJ1ebzWeqj19iatdU2tdW3r909tEVwI4by8qlZl77BhFN2+naygILBYICwM5s0zO5qw8jQ7gBAif6reoAHZSuGRlGQsOHYM/nsdQL9+5gUTgAzvIISwE48XX8Rye4+flBSYMMGcQOIWUvyFEPZx/Pj9LRcOJcVfCGEfoaH3t1w4lBR/IYR9TJkC/v63LNJA4ogR5uQRt5DiL4Swj379YOZMKFcOlEIXL44Gzi5YwJXLl81O5/ak+Ash7KdfPzh6FLKzUWfPkjpuHNW2b+e3v/+djIwMs9O5NSn+QgiHKTBlCikRETSbN4+4jz5Cy/g/ppHiL4RwHA8P/BcvRvn5Ueu117gsc/+aRoq/EMKxQkLwnDOHkidPEvTmm2ancVtS/IUQDqe6dYOnnoKpU9k9dSoJCQlmR3I7UvyFEOZ46y2yq1en/MSJLJ8xg0uXLpmdyK1I8RdCmMPPD8uiRfhlZxP9+efM/+wzUlNTzU7lNqT4CyHMU7Uq6tNPCTl2jLoLFrBgwQIyMzPNTuUWpPgLIczVpw+MHEnjTZso8N13HJexfxxCir8QwnxTp0KjRvRcuZIKcuTvEFL8hRDm8/aGhQtRPj7QowcHduzghx9+MDtVvibFXwjhHEJD4YsvYPduCjz7LD//9BMbNmwwO1W+JcVfCOE82raFSZMoExdHt6NHWbduHVu3bjU7Vb4kxV8I4VwmTIBevaj52WdEpqaycuVKeQOwAyn+QgjnohTMno2qVYsWH35IRMGCXLhwwexU+Y4UfyGE8ylQAJYuRfn40PGjj2jfqBEAaWlpJgfLP6T4CyGcU2goLF6MOnoU1acPVxITmTZtGps2bTI7Wb4gxV8I4byaNYMZM2DNGgJeeIGyISGsWbOGtWvXylwAeeRpdgAhhPhLQ4bA/v1Y3niDHuXK4Ve/Phs3buTatWtER0djscgx7IOQV00I4fxefRX69sXyz3/S6fJlWrZsyY4dO/j555/NTuay5MhfCOH8LBaYPRtOn0YNGULkqlUU7d6d8PBws5O5LDnyF0K4Bh8fWLwYwsOhWzdqZmXh4+NDRkYGCxYs4PTp02YndClS/IUQrqNwYfjuOyhUyLgaeN8+kpKSOHXqFLNnz2bPnj1mJ3QZUvyFEK6lbFlYt85oCoqKouiVKwwbNozixYuzaNEi1qxZQ1ZWltkpnZ4UfyGE6wkPhzVrIDUVWrem4OXLDBo0iPr167Np0ya+/fZbsxM6PTnhK4RwTTVrwurVEBUFrVvjuX49nTt3JiwsjFKlSgGQnZ0tXUHvQl4VIYTratAAvv0WTpyAli3hxAlq1KhB0aJFAfjmm29YunQp169fNzmo85HiL4RwbQ8/DN9/D2fOGLePHAFAa01QUBC///4706dP59ChQyYHdS5S/IUQrq9ZM+MkcHKy8Qawbx9KKaKionjiiSfw9fVl3rx5LFq0iKtXr5qd1ilI8RdC5A8REbB+PWRlGW8AW7YAUKZMGWJjY4mMjOTEiRMopczN6SSk+Ash8o+aNeHnnyEwEFq1gq+/BsDT05OWLVsyatQoChQogNaaL774gq1bt5LpTBPGz5sHYWFGN9awMOO+nUjxF0LkLw89BJs3Q9260KsXTJ0K1hFAPT2NDo6pqamkpaWxcuVKPvjgA+Lj48nIyDAztVHoY2Ph2DEj77Fjxn07vQEoZx0WNSIiQsfHx5sdQwjhqlJTYeBAWLTI+D59Ovj733hYa82RI0eIi4vj5MmT+Pr6MnDgQEqWLGlO3rAwo+Dfrlw5OHo010+jlNqmtY6413p56uevlOoFTASqAg211jlWa6VUe+A9wAP4RGv9Wl62K4QQ9+TnB/PnQ7VqMGkS7NhhNANVrAiAUoqKFStSoUIFTpw4we+//05wcDAA27dvJysri2rVquF/0xuGraWmprJnzx4OHjxI7+PHyfFsxPHjdtl2Xi/y2gV0Bz662wpKKQ9gGtAGSAC2KqWWaa1lEA4hhH1ZLDBxIjRsCI8/DvXrw8cfG81BVkopQkNDCQ0NvbFs7969HDx4kJUrV1K6dGkqVqxI5cqVKV26dJ4jXb58mR07dnDkyBFOnjxJdnY25TMzjaw5DUtxUy5bylPx11rvBe519rwhcEhrfcS67nygCyDFXwjhGB07wrZt0Ls3PPooPPYYfPABBAXluHrfvn05e/Ys+/bt4/Dhw/z8889cvHiRHj16oLVmxYoV+Pv7U7hwYQoVKoS3tzdBQUEUKlSIzMxMTp06RVpaGklJSSQlJXHp0iVq165NxYoVSUpK4scff6R06dI0bdqUOleuUGT4cJSfH2Rmws0XpPn7w5QpdnlJHDG8QxngxE33E4BGDtiuEEL8T/nysGmTMTHMpElGt9Bp06BLF7jtAFYpRcmSJSlZsiSRkZGkpqbeuEr4+vXrHD58mKSkpFumkmzRogWtWrUiNTWV2bNn3/JcAQEBVKhQAYCQkBDGjh2Lv58fvPsujB0LlSrB0qUQHw8TJhhNPaGhRuHv188uL8c9T/gqpdYCOZ0BmaC1XmpdZz0wJqc2f+t5gXZa66HW+/0xzg88ncO6sUAsQGhoaP1jOZ38EEKIvNq2DQYNgl27oH17eP99o5fQfcjOziY5OZnk5GTS09MJDAykaNGiZGVlcfToUby9vSlUqBABAQF3ji+UkABDhxpjE3XrBnPmGMNU24DNTvhqrVvnMUsCUPam+yHAqbtsayYwE4zePnncrhBC5Kx+ffjtN+PI/6WXoHp1eOIJeP55Y8joXLBYLAQGBhIYGHjLcg8PDypaTyrfISMDPvwQXnzRuD1tGowYcccnD0dwRD//rcBDSqnySilvoA+wzAHbFUKIu/PygtGjYf9+4yj800+N5pfhw+GPP2y7rcxMWLDAuAjtmWeMAel27oSnnjKl8EMei79SqptSKgFoAnyrlFptXV5aKbUSQGudCYwEVgN7gYVa6915iy2EEDZSsqRxDcDBg0ZT0Jw5UKuWMUTEp5/C+fMP/twJCfDmm8b8A336GIV++XJjILq7fTpwELnISwghbpaYaLwBzJgBhw4ZXTCbNoXmzaFJE6OJKDTU+ORwM62NoaX37IENG4yB5n75xVjerBmMGQMxMcbz2VFu2/yl+AshRE60Ni4MW7zYmDd4xw6j+QaMAl6smHEhmbe3MZroxYuQnm487uFhNO20b2/01qlUyWGxpfgLIYQtpaQYJ4kPHoQ//zTmD0hLMwp+QIBxzUD58lC1qjGukI1679wvhwzvIIQQbsPf32j6ad7c7CQ2IaN6CiGEG5LiL4QQbkiKvxBCuCEp/kII4Yak+AshhBuS4i+EEG5Iir8QQrghKf5CCOGGnPYKX6XUeSAvA/oXAy7YKI6rcLd9drf9Bdlnd5GXfS6ntQ6+10pOW/zzSikVn5tLnPMTd9tnd9tfkH12F47YZ2n2EUIINyTFXwgh3FB+Lv4zzQ5gAnfbZ3fbX5B9dhd23+d82+YvhBDi7vLzkb8QQoi7cOnir5Rqr5Tar5Q6pJQan8PjPkqpBdbHf1FKhTk+pW3lYp+fVUrtUUrtVEqtU0qVMyOnLd1rn29ar6dSSiulXL5nSG72WSn1qPV3vVsp9YWjM9paLv62Q5VScUqp7da/745m5LQVpdQspdQ5pdSuuzyulFLvW1+PnUqpejYNoLV2yS/AAzgMVAC8gd+Baret8xQww3q7D7DA7NwO2OdWgL/19gh32GfregHAT8AWIMLs3A74PT8EbAeCrPeLm53bAfs8ExhhvV0NOGp27jzucwugHrDrLo93BL4DFNAY+MWW23flI/+GwCGt9RGtdTowH+hy2zpdgLnW218BUUop5cCMtnbPfdZax2mtU6x3twAhDs5oa7n5PQO8ArwBXHdkODvJzT4PA6ZprS8BaK3POTijreVmnzXw37kRA4FTDsxnc1rrn4CLf7FKF+A/2rAFKKyUKmWr7bty8S8DnLjpfoJ1WY7raK0zgStAUYeks4/c7PPNnsA4cnBl99xnpVRdoKzWeoUjg9lRbn7P4UC4UmqjUmqLUqq9w9LZR272eSLwuFIqAVgJPO2YaKa53//3++LKc/jmdAR/e9el3KzjSnK9P0qpx4EIoKVdE9nfX+6zUsoCvAMMclQgB8jN79kTo+knEuPT3c9KqRpa68t2zmYvudnnvsAcrfVUpVQT4DPrPmfbP54p7Fq/XPnIPwEoe9P9EO78GHhjHaWUJ8ZHxb/6mOXscrPPKKVaAxOAGK11moOy2cu99jkAqAGsV0odxWgbXebiJ31z+7e9VGudobX+E9iP8WbgqnKzz08ACwG01psBX4wxcPKrXP2/PyhXLv5bgYeUUuWVUt4YJ3SX3bbOMmCg9XZP4AdtPZPiou65z9YmkI8wCr+rtwPDPfZZa31Fa11Max2mtQ7DOM8Ro7WONyeuTeTmb3sJxsl9lFLFMJqBjjg0pW3lZp+PA1EASqmqGMX/vENTOtYyYIC1109j4IrW+rStntxlm3201plKqZHAaoyeArO01ruVUpOAeK31MuBTjI+GhzCO+PuYlzjvcrnPbwIFgUXWc9vHtdYxpoXOo1zuc76Sy31eDbRVSu0BsoCxWutE81LnTS73+R/Ax0qpv2M0fwxy5YM5pdSXGM12xaznMV4CvAC01jMwzmt0BA4BKcBgm27fhV87IYQQD8iVm32EEEI8ICn+QgjhhqT4CyGEG5LiL4QQbkiKvxBCuCEp/kII4Yak+AshhBuS4i+EEG7o/wGEgxiL1UJz2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pt.plot(x_01, f(x_01), \"--\", color=\"gray\", label=\"$f$\")\n",
    "pt.plot(x_0h, interp_0h, color=\"red\", label=\"Interpolant\")\n",
    "pt.plot(nodes, f(nodes), \"or\")\n",
    "pt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now build the gen. Vandermonde matrix $V'=$`Vprime` of the derivatives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "def monomial_deriv(i, x):\n",
    "    if i == 0:\n",
    "        return 0*x\n",
    "    else:\n",
    "        return i*nodes**(i-1)\n",
    "\n",
    "Vprime = np.array([\n",
    "    monomial_deriv(i, nodes)\n",
    "    for i in range(degree+1)\n",
    "]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the value of the derivative at the nodes as `fderiv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#clear\n",
    "fderiv = Vprime @ la.inv(V) @ f(nodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot vs `df`, the exact derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5f98bc9d68>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlUVmeC5/Hvw74IgqKIIGBporiGSFyyjCl3k0iMFY1pYsVKpzxVM6nTdaame7o7f/Sc6ZOZOTNnupPpSqfGjlXJWKillgE1iSSWWkmMGlExxjXuhbjgBggi2zN/vGipQUF5eZ93+X3O8QD3vd77u6I/L/c+73ONtRYREQkeYa4DiIiId6nYRUSCjIpdRCTIqNhFRIKMil1EJMio2EVEgoyKXUQkyKjYRUSCjIpdRCTIRLjYaUpKis3OznaxaxGRgLVjx47z1tpe7a3npNizs7MpLS11sWsRkYBljDnRkfV0KUZEJMio2EVEgoyKXUQkyHjlGrsxJgl4FxgGWOAVa+0Wb2xbRMRbGhsbKS8vp76+3nWUu4qJiSEjI4PIyMj7+v3eunn6FrDOWvu8MSYKiPPSdkVEvKa8vJyEhASys7MxxriO0yZrLRcuXKC8vJz+/fvf1zY6fSnGGJMI/DtgUWuoBmvt5c5u9zsKCyE7G8LCPB8LC72+CxEJbvX19fTs2dNvSx3AGEPPnj079VOFN66xfw+oBH5jjNlljHnXGBPvhe3+WWEhLFgAJ06AtZ6PCxao3EXknvlzqV/X2YzeKPYI4GHgHWttLlAL/O3tKxljFhhjSo0xpZWVlfe2h9dfh7q6W5fV1XmWi4jILbxR7OVAubV2W+vXK/EU/S2stQuttXnW2rxevdp949StTp68t+UiIiGs08VurT0D/MkYM6h10URgX2e3e4vMzHtbLiISwrw1jv1nQKEx5mvgIeC/eWm7Hm+8AXG3DbSJi/MsFxEJQL/4xS8YMmQIP/vZz7y+ba8Md7TWlgF53thWmwoKPB9ff91z+SUz01Pq15eLiASQo0ePsnnzZvbt8+7FjeucTAJ2XwoKVOQiEvAOHjzIpEmTaGpqIjc3ly+++IL4eO8OJAycYhcR8bL33nvvO8uGDh3KI488QmNjI4VtDKl+6KGHeOihh6irq2P58uW3vDZ//vx29zlo0CBefvllsrOzefXVV+83+l1prhgRER/bs2cPI0eO7LLt64xdRELW3c6wIyMj7/p6XFxch87Q27J3716GDh16X7+3I3TGLiLiQzU1NURGRhJ3+0g/L1Kxi4j40DfffMOwYcO6dB8qdhERHxo3bhwrVqzo0n2o2EVEgoyKXUQkyKjYRSSkWGtdR2hXZzOq2EUkZMTExHDhwgW/LvfrT1CKiYm5721oHLuIhIyMjAzKy8u552dC+Nj1Z57er4Ar9n/7t3/j0qVLREREEB0dTWJiItnZ2TzxxBMA1NbWEhcXFxBPSRER34qMjLzv54gGkoAr9iFDhlBVVUVjYyPXrl2jqqqKqqoqwPMjzNtvvw1AdnY2/fv353vf+x49e/Z0GVlExKcCrtgfe+yxO75mrWXChAmcOnWKY8eOsX//fgCefPJJxo8f76uIIiJOBVyx301YWBh5eXnk5eVhreXSpUt8++23ZLY+aenUqVOUlJQwatQohgwZQmRkpOPEIiLeF1TFfjNjDD169GDMmDE3ltXX11NXV0dRURElJSWMHj2aMWPGEBsb6zCpiIh3GRfDfvLy8mxpaanP9wueyzXHjx9n27ZtHDx4kISEBH7+858TFqaRnyLi34wxO6y17T6tLmjP2O/EGEP//v3p378/586do7KykrCwMKy17Nq1i+HDh+sSjYgEtJA+Te3du/eNOZFPnjzJmjVr+Jd/+Rd2797t129gEBG5m5Au9ptlZWUxf/58EhISKCoqYvHixVy4cMF1LBGRe6Ziv0lWVhavvvoqTz/9NBUVFSxZskRn7iIScELuGnt7jDHk5eUxaNAgqqurMcbQ3NxMdXU1ycnJruOJiLRLZ+x3kJCQQHp6OgBffPEF77zzDjt37tQZvIj4PRV7B+Tm5pKRkcGaNWtYsWIF9fX1riOJiNyRir0DEhMTmTdvHpMnT+bgwYO8++67fj87nIiELhV7BxljePTRR/nhD39Ic3MzTU1NriOJiLTJa8VujAk3xuwyxqz11jb9UVZWFq+99hppaWkAfPvtt7ruLiJ+xZtn7H8F7Pfi9vxWeHg4ACdOnGDJkiWsXLmSxsZGx6lERDy8UuzGmAzgaeBdb2wvUGRmZjJlyhT27dvH+++/z5UrV1xHEhHx2hn7m8DfAC1e2l5AMMYwbtw45syZw9mzZ1m0aJHerSoiznW62I0xzwDnrLU72llvgTGm1BhTGmwjSnJycpg/fz6NjY1UVFS4jiMiIa7T0/YaY/47MA9oAmKARGCVtfalO/0el9P2dqVr164RHR39nc9FRLyho9P2dvqM3Vr7d9baDGttNjAX2HC3Ug9m14v8xIkTvPnmmxw+fNhxIhEJRRrH3gVSUlJISkpi6dKlHDhwwHUcEQkxXi12a+0ma+0z3txmIIqPj+fll18mLS2NFStWcPDgQdeRRCSE6Iy9i8TExPDSSy/Rp08fVqxYwZkzZ1xHEpEQoWl7u1BMTAzz5s2jtLSU1NRU13FEJETojL2LxcTE8Pjjj2OM4fLlyxw7dsx1JBEJcip2H/r4449ZunQp5eXlrqOISBBTsfvQjBkz6NatG4WFhZw7d851HBEJUip2H+rWrRvz5s0jIiKC3/72t1y+fNl1JBHxlcJCyM6GsDDPx8LCLtuVit3HkpOTmTdvHo2Njaxfv951HBHxhcJCWLAATpwAaz0fFyzosnLv9JQC9yNYpxS4F2fOnCE5OVnTDoiEguxsT5nfLisLjh/v8GZ8NqWA3J8+ffoQHR1NQ0MDX3zxBS0tITUxpkhoOXny3pZ3kordsUOHDvGHP/yBdevW6UlMIsEqM/PelneSit2xYcOGMW7cOLZv387WrVtdxxGRrvDGGxAXd+uyuDjP8i6gYvcDkydPJicnh08++YT9+0Pi6YIioaWgABYu9FxTN8bzceFCz/IuoCkF/IAxhueee46amhrWrVvHAw88QESEvjUiQaWgoMuK/HZqDz8RGRnJ3LlzaWhoUKmLSKfoUowfiY+PJzk5GWstO3bsoKmpyXUkEQlAKnY/VF5eztq1a/nwww81UkZE7pmK3Q/169eP8ePHU1ZWxrZt21zHEZEAo2L3U+PHj2fw4MF88sknHDlyxHUcEQkgKnY/dX2kTK9evVi1ahUNDQ2uI4lIgNDwCz8WFRXF3Llzqa6uJioqynUcEQkQOmP3c8nJyWRlZQFw+vRp3UwVkXap2APE0aNHWbhwIbt27XIdRUT8nIo9QGRnZzNgwAA++ugjKioqXMcRET+mYg8QYWFhzJo1i/j4eJYvX87Vq1ddRxIRP6ViDyBxcXHMmTOHmpoaVq1apevtItImjYoJMOnp6Tz11FMqdRG5IxV7ABo1atSNz1taWggL0w9eIvJnnW4EY0w/Y8xGY8x+Y8xeY8xfeSOYtO/AgQO888471NXVuY4iIn7EG6d6TcAvrLU5wFjgPxhjhnhhu9KO7t27c+nSJYqLi3VpRkRu6HSxW2tPW2t3tn5eA+wH0ju7XWlfWloakydP5tChQ3qsnojc4NWLs8aYbCAX0JSEPjJ69GgGDx7M+vXrNb5dRAAvFrsxphvwe+Dn1trqNl5fYIwpNcaUVlZWemu3Ic8YQ35+PgkJCRw6dMh1HBHxA8Yb12aNMZHAWqDEWvtP7a2fl5dnS0tLO71f+bOrV68SGxvrOoaIdCFjzA5rbV5763ljVIwBFgH7O1Lq0jWul/rp06fZu3ev4zQi4pI3xrE/BswD9hhjylqX/b219iMvbFvu0aZNmzh69CipqamkpKS4jiMiDnhjVMwX1lpjrR1hrX2o9ZdK3ZFnnnmGyMhIVq1aRXNzs+s4IuKA3rIYZBISEpgxYwanT59m48aNruOIiAMq9iCUk5NDbm4umzdvpry83HUcEfExzRUTpKZNm0ZKSgppaWmuo4iIj+mMPUhFRUXx6KOPEh4ezrVr11zHEREfUrEHuYsXL/L222/z9ddfu44iIj6iYg9ySUlJJCUl8dFHH3H58mXXcURC1pEjRygrK/PJhH0q9mBVWAjZ2YRFRPDyP/wDObt28cEHH9DS0uI6mUjIqauro6ioiC+//NIn/wZV7MGosBAWLIATJ8BawsvLmbF6NYlr17J582bX6URCTlVVFZGRkcyaNYvw8PAu35+KPRi9/jrc9vCNsPp6pn32GadOndLc7SI+lpaWxmuvvUafPn18sj8VezA6ebLNxXEXLvDCCy/gmd5HRLra5cuX2bhxI83NzT59hKWKPRhlZra52GRmYozh8uXL7Nq1y8ehREKLtZbVq1ezdetWampqfLpvFXsweuMNiIu7dVlcnGc5sGXLFlavXs3JO5zZi0jnbd++nWPHjjF16lSSkpJ8um8VezAqKICFCyErC4zxfFy40LMcmDBhAklJSRQVFdHQ0OA4rEjwuXDhAp9++ikDBw4kNzfX5/tXsQerggI4fhxaWjwfW0sdIDo6mpkzZ3Lp0iU++eQTZxFFgtXatWuJiIggPz/fyT0tzRUTorKyshg3bhxbtmwhJyeHAQMGuI4kEjSmTp1KTU0NCQkJTvavYg9hEyZMICIigoyMDNdRRIJCQ0MDUVFR9OnTx2dDG9uiSzEhLCIiggkTJhAdHa2Hcoh0UnNzM++9955fXN5UsQvV1dUsXLhQz0oV6YTPP/+c06dP069fP9dRVOwC8fHxRERE8OGHH3LlyhXXcUQCTkVFBZ999hkjRowgJyfHdRwVu0B4eDgzZ86ksbGRNWvWaMoBkXvQ1NREUVER3bp1Y9q0aa7jACp2adWrVy8mTpzIoUOHKCsrcx1HJGCcPXuW6upq8vPziY2NdR0HULHLTcaMGUN2dja7d+/WWbtIB6Wnp/Pzn/+cgQMHuo5yg4Y7yg3GGGbPnk10dLQmChNpR0NDA/v27WPkyJHExMS4jnMLnbHLLeLi4ggPD+fq1ascOXLEdRwRv7V+/XqKi4s5c+aM6yjfoWKXNpWUlLBs2TIuXLjgOoqI3zl69Cjbt29n7NixpKWluY7zHSp2adPEiROJiIigqKhIj9MTuUl9fT3FxcX07NmTCRMmuI7TJhW7tCkhIYGnnnqK8vJyvvzyS9dxRPxGSUkJNTU1zJw5k8jISNdx2qSbp3JHw4YNY//+/WzatIkHHniA1NRU15FEnBs2bBi9e/f26zmWvFLsxphpwFtAOPCutfZ/eGO74pYxhqeffpqIiAi/GZ8r4oq1FmMMAwYM8PvZUDt9KcYYEw68DUwHhgAvGmOGdHa74h/i4+OZNWsWiYmJrqOIOFVcXMymTZtcx+gQb1xjHw0cttYetdY2AMuAZ72wXfEjNTU1/Pa3v6WiosJ1FBGf279/P7t37w6Y93d4o9jTgT/d9HV567JbGGMWGGNKjTGllZWVXtit+FJERATnzp3jgw8+oKmpyXUcEZ+pra1l7dq1pKWl8fjjj7uO0yHeKPa2/gv7zvvRrbULrbV51tq8Xr16eWG34kuxsbHk5+dz/vx5NmzY4DqOiE9Ya/nwww+5du0aM2fOJDw83HWkDvFGsZcDN09AnAHo5/UgNHDgQEaNGsWWLVs4efKk6zgiXe78+fMcPHiQJ598kt69e7uO02HeKPbtwAPGmP7GmChgLrDaC9sVPzR58mSSkpIC5iaSSGf06tWLn/zkJzz66KOuo9yTTg93tNY2GWNeA0rwDHf8tbVWj+IJUtHR0bz44ot0797ddRSRLmOt5dSpU2RkZBCIl4698s5Ta+1H1toHrbUDrLVveGOb4r969+5NdHQ0TU1NnD9/3nUcEa/bvXs3ixYt4vDhw66j3BdNKSD3bdWqVSxevJj6+nrXUUS8prq6mnXr1pGZmen3b0S6ExW73LfHHnuMmpoaSkpKXEcR8QprLatXr6alpYVnn302YMat307FLvctPT2dxx9/nLKyMg4ePOg6jkin7dy5kyNHjjBp0iR69OjhOs59U7FLp4wfP57U1FTWrFlDXV2d6zginRIVFUVOTg6PPPKI6yidomKXTgkPD+e5554jOTlZ19ol4A0fPpw5c+YE7CWY6zRtr3Raamoqr7zySsD/Y5DQtWPHDlpaWsjLywuKv8c6YxevMMZw9epV1q5dS01Njes4Ih128eJFSkpKOHTokOsoXqNiF6+pra1l9+7drFmzBmu/M12QiN9paWmhuLiYsLAwZsyYERRn66BiFy9KSUlh4sSJfPvtt5SVlbmOI9Kubdu2cfLkSaZPnx5UzxxQsYtXjRkzhqysLNatW8fly5ddxxG5o9raWjZs2MCgQYMYMWKE6zhepWIXrzLG8Oyznues6I1L4s/i4+OZO3cuzzzzTNBcgrlOo2LE65KTk5k9e3ZATXMqoaWmpoaEhISAnTKgPTpjly4xcOBAEhMTsdZy9epV13FEbjh16hRvvfUW+/fvdx2ly6jYpUutWLGCJUuW0NLS4jqKCA0NDaxatYr4+Hj69+/vOk6XUbFLlxo8eDDl5eVs2bLFdRQRPv30Uy5evMjMmTOJiYlxHafLqNilSw0fPpycnBw2btzIuXPnXMeREHb48GFKS0sZO3ZsUJ+tg4pdupgxhqeffpro6Gg++OADmpubXUeSEHXp0iX69OnDxIkTXUfpcip26XLx8fHMmDGDq1evUlVV5TqOhKhHHnmEH//4x0REBP9gwOA/QvELgwcPZsCAAURGRrqOIiFm//79hIWFMWjQIMLCQuNcNjSOUvxCZGQkTU1NbN68maamJtdxJARUV1ezevVqvvjii5Cav0jFLj5VXl7O+vXr2bBhg+soEuSstRQXF9Pc3Mxzzz0XdO8uvRsVu/hUdnY2o0aNYsuWLRw9etR1HAli27Zt4+jRo0yZMiWgH3N3P1Ts4nNTpkwhJSWFoqIiPU5PusSlS5dYv349gwYNYtSoUa7j+JyKXXwuKiqKWbNmUVdXx0cffeQ6jgShpKQkZsyYQX5+fkhdgrlOo2LEibS0NPLz8+nTp4/rKBJkamtriY+PZ+TIka6jOKMzdnFmxIgRN2aAvHbtmuM0EgwOHDjAW2+9RXl5uesoTqnYxbmSkhJ+85vfaAikdMr1oY0pKSmkpaW5juNUp4rdGPO/jDEHjDFfG2M+MMYkeSuYhI7+/ftz9uxZ1q9f7zqKBChrLUVFRTQ1NTFr1izCw8NdR3Kqs2fsnwLDrLUjgEPA33U+koSaBx98kNGjR7Nt2zYOHz7sOo4EoC+//JJjx44xbdo0UlJSXMdxrlPFbq39xFp7/efnrUBG5yNJKJo8eTK9e/emqKiIK1euuI4jAebKlSvk5OSQm5vrOopf8OY19leAj724PQkhERER/OAHP8AYw/nz513HkQAzdepUnn/++ZAc2tgW0978CcaY9UBbY9Jet9YWt67zOpAHzLJ32KAxZgGwACAzM3PUiRMnOpNbglRTU1NIzL4nnWet5ZNPPmH48OH07dvXdRyfMMbssNbmtbdeu/+CrLWT2tnRy8AzwMQ7lXrrdhYCCwHy8vJCZzYeuScRERFYa9m2bRt9+/YlMzPTdSTxU6WlpWzdupX4+PiQKfaO6uyomGnAfwbyrbV6b7h4RWNjI1999RUrV67UlAPSptOnT1NSUsLAgQN57LHHXMfxO529xv5LIAH41BhTZoz5lRcySYiLiopi9uzZ1NXV8cEHH4TUdKvSvmvXrrFixQri4uJCbtbGjursqJiB1tp+1tqHWn/9xFvBJLSlpaUxdepUDh8+zObNm13HET+ydetWLl++zPPPP09cXJzrOH5Jd6nEb+Xl5XHixAk2btzIiBEjSExMdB1J/MATTzxBdna27r/chYpd/JYxhhkzZpCXl6dSFyorK4mNjaVbt25kZWW5juPXNFeM+LXo6Giys7MBz9OXWlpa3AYSJ+rq6liyZAlLly7VPZcOULFLQDh79iyLFi3SI/VCUEtLC6tWraKmpoannnpKN0s7QMUuASE1NZWHH36YzZs3s2/fPtdxxIc2bdrEkSNHmD59Ounp6a7jBAQVuwSM6dOnk5GRQXFxMZWVla7jiA98++23fP755+Tm5vLwww+7jhMwVOwSMCIiIpg9ezaRkZH87ne/o6GhwXUk6WIZGRmMGzdOl2DukUbFSEBJTExk9uzZVFRUEBkZ6TqOdJGGhgbCwsKIjY1lypQpruMEHBW7BJysrKwbw93q6+uJiYlxnEi86fpDM65cucL8+fMJC9OFhXulPzEJWBUVFbz11lscOHDAdRTxok2bNrF//35ycnJU6vdJf2oSsHr16kWPHj1YtWoVZ86ccR1HvOCbb77hs88+46GHHmLs2LGu4wQsFbsErMjISObOnUtsbCxLly7Vk5cC3KlTpyguLiYzM5NnnnlGN0s7QcUuAS0hIYG5c+dy9epVli1bRmNjo+tIcp+io6PJyspizpw5If8w6s5SsUvAS0tLY9asWaSkpOgsLwA1NzdjrSUlJYWXXnqJ+Ph415ECnkbFSFAYPHgwgwcPBvR4vUDS0tLCypUriYmJIT8/X/8xe4nO2CWoVFdX86tf/YqysjLXUaQd1lrWrVvHgQMHSE1NVal7kYpdgkp8fDyJiYmsXr2aw4cPu44jd7Flyxa2b9/OuHHjNALGy1TsElTCw8N54YUXSE1NZfny5VRUVLiOJG3Ys2cPn376KUOHDmXy5Mmu4wQdFbsEnejoaP7iL/6C+Ph4lixZwsWLF11HktvExcXx4IMPMnPmTF2C6QK6wyRBKSEhgYKCAkpKSoiOjnYdR1pdnwJiwIABDBgwwHWcoKUzdglaKSkpFBQUEB8fT1NTE/X19a4jhbTKykp++ctf6sa2D6jYJehZa1m5ciWFhYWa6teRS5cusXjxYgD69evnOE3wU7FL0DPGMHLkSE6dOsXSpUv17lQfq6mpYfHixTQ2NjJv3jx69uzpOlLQU7FLSMjJyWHmzJkcP36cFStW0NzcfOsKhYWQnQ1hYZ6PhYUuYgadpqYmFi9eTG1tLQUFBaSmprqOFBJ081RCxogRI2hsbGTt2rV8+OGH5Ofne14oLIQFC6CuzvP1iROerwEKCtyEDRIRERHk5ubSp08fMjIyXMcJGcZa6/Od5uXl2dLSUp/vVwSgtLSU9PR00tLSPAuysz1lfrusLDh+3JfRgsaVK1eorq6mb9++rqMEFWPMDmttXnvr6VKMhJy8vLwbpb57927syZNtr3in5XJXNTU1vP/++yxbtoympibXcUKSV4rdGPOfjDHWGJPije2J+MLZs2cpKiqitkePtlfIzPRtoCBQVVXF+++/T1VVFT/4wQ80GZsjnS52Y0w/YDKg0xsJKKmpqcyaNYuS8eNpjIq69cW4OHjjDTfBAlRlZSW//vWvuXLlCi+99NKN59KK73njjP2fgb8BfH+xXqSThg8fTs4//iNr8/OpTk7GGuO5tr5woW6c3qNt27bR3NzM/PnzydRPO0516uapMSYfmGit/StjzHEgz1p7vr3fp5un4m+OHTvG8uXLyc/PJycnx3WcgNLS0kJYWBjNzc1cuXKF7t27u44UtDp687TdYjfGrAf6tPHS68DfA1OstVXtFbsxZgGwACAzM3PUibZGIYg4dPXqVWJjYwGoq6sjLi7OcSL/t2PHDr766itefvll/Xn5gNdGxVhrJ1lrh93+CzgK9Ad2t5Z6BrDTGNPWfwJYaxdaa/OstXm9evW6t6MR8YHrpX7ixAnefPNNdu7c6TiR/2ppaWHdunWsXbuWhIQEPaPUz9z3LWtr7R6g9/Wv7+VSjIg/6927N/369WPNmjWcPXuWKVOmqLhucu3aNVauXMnhw4cZM2YMU6ZMISxMI6f9ib4bIreJjY2loKCAsWPH8tVXX1FYWEjd9XelCuvWrePo0aM8/fTTTJs2TaXuh7w2yNRam+2tbYm4FhYWxtSpU+nTpw9r1qxhz549jBkzxnUsp64/JPz73/8+I0eOJDs723UkuQO9e0DkLkaOHEnfvn1JSfG89+7ChQv06NEjpJ7609TUxMcff0xVVRUFBQUkJiaSmJjoOpbchYpdpB3Xb/bX1dWxaNEi0tLSmDlzJgkJCY6Tdb2KigqKioqorKzksccew1obUv+pBSpdHBPpoNjYWCZNmsTJkyf513/9V8rKynAxiZ4vNDc388c//pFFixZRX19PQUEBkyZN+vP1dE1z7Nd0xi7SQcYYHn74YbKysli9ejXFxcXs3buXF154IejmRGlsbGTnzp0MHTqU6dOn3xgKCmia4wCgaXtF7oO1lq+++oozZ87w7LPP3lgWyJcp6uvr2bp1K0888QTh4eHU1tYSHx//3RU1zbEzHX2DUnCdZoj4iDHmllEy58+fZ/ny5UyaNIkHHnggoAq+ubmZXbt2sWnTJurq6ujXrx8DBgxou9ThztMZa5pjv6FiF/GCa9eu0dLSwtKlS8nKyuLJJ5/0++GA1lr27dvHhg0buHjxIpmZmRQUFPz5ASR3kpnZ9hm7Jv7yG7p5KuIF6enp/PSnP2X69OlcuHCB999/nyVLlvjlzdWbM3355ZeEh4fz4osvMn/+/PZLHTzTGd8+L4ymOfYrOmMX8ZLw8HBGjx5Nbm4uO3bsoKGhAWMM1lr27t3LoEGDiIyMdJavtraWsrIydu3axSuvvEJcXBxz584lPj7+3t49ev0G6euvey6/ZGZ6Sl03Tv2Gil3EyyIjIxk7duyNr0+dOsXvf/97oqOjGTJkCMOHDycrK8snb8VvamriyJEjlJWVcejQIVpaWsjKyroxe+V9j8UvKFCR+zEVu0gXS09P50c/+hE7duxg79697Nq1i27dujFv3jx69+7t1dE01lqqqqpobGykV69e1NbWsmzZMuLj4xkzZgy5ublodtXgp2IX6WLGGDIzM8nMzKSxsZFDhw5x4MABkpOTAdi0aRN79+6lb9++9O3blx49epCcnNzhAj506BCnT5+msrKSP/3pT1RXV5NLkZ4UAAAFAElEQVSTk8OcOXPo3r07P/rRj0hPT9cMlSFExS7iQ5GRkQwdOpShQ4feWNajRw9SUlI4fvw4e/bsASAuLo6//uu/BqCoqIiTrUMJW1paaG5upnv37rz66qsA/PGPf6SiooKkpCT69etHZmYm/fv3v7F9PaYu9KjYRRwbOXIkI0eOBODKlStcunSJ+vr6G6/37t2b5uZmjDGEhYURFhZGjx49brw+e/Zs4uLiiLr9gdwSslTsIn6kW7dudOvW7ZZljz766F1/T1JSUldGkgCkcewiIkFGxS4iEmRU7CIiQUbFLiISZFTsIiJBRsUuIhJkVOwiIkFGxS4iEmScPBrPGFMJtDFTf4ekAOe9GCcQ6JhDg445NHTmmLOste1OIuSk2DvDGFPakWf+BRMdc2jQMYcGXxyzLsWIiAQZFbuISJAJxGJf6DqAAzrm0KBjDg1dfswBd41dRETuLhDP2EVE5C78ttiNMdOMMQeNMYeNMX/bxuvRxpjftb6+zRiT7fuU3tWBY/6Pxph9xpivjTF/MMZkucjpTe0d803rPW+MscaYgB5B0ZHjNcbMaf0+7zXGLPF1Rm/rwN/rTGPMRmPMrta/20+5yOlNxphfG2POGWO+ucPrxhjzf1r/TL42xjzs1QDWWr/7BYQDR4DvAVHAbmDIbev8e+BXrZ/PBX7nOrcPjvn7QFzr5z8NhWNuXS8B+AzYCuS5zt3F3+MHgF1AcuvXvV3n9sExLwR+2vr5EOC469xeOO5/BzwMfHOH158CPgYMMBbY5s39++sZ+2jgsLX2qLW2AVgGPHvbOs8C77d+vhKYaLz1qHc32j1ma+1Ga21d65dbgQwfZ/S2jnyfAf4R+J9AfRuvBZKOHO+PgbettZcArLXnfJzR2zpyzBZIbP28O1Dhw3xdwlr7GXDxLqs8C/w/67EVSDLGpHlr//5a7OnAn276urx1WZvrWGubgCqgp0/SdY2OHPPN/hLP//iBrN1jNsbkAv2stWt9GayLdOR7/CDwoDFmszFmqzFmms/SdY2OHPN/AV4yxpQDHwE/8000p+713/s98ddnnrZ15n378J2OrBNIOnw8xpiXgDxgfJcm6np3PWZjTBjwz8B8XwXqYh35HkfguRzzJJ6fyD43xgyz1l7u4mxdpSPH/CLwnrX2fxtjxgGLW4+5pevjOdOl/eWvZ+zlQL+bvs7guz+e3VjHGBOB50e4u/3o4+86cswYYyYBrwP51tprPsrWVdo75gRgGLDJGHMcz7XI1QF8A7Wjf6+LrbWN1tpjwEE8RR+oOnLMfwksB7DWbgFi8MynEsw69O/9fvlrsW8HHjDG9DfGROG5Obr6tnVWAy+3fv48sMG23pUIUO0ec+tlif+Lp9QD/dortHPM1toqa22KtTbbWpuN575CvrW21E3cTuvI3+siPDfJMcak4Lk0c9SnKb2rI8d8EpgIYIzJwVPslT5N6XurgR+2jo4ZC1RZa097beuu7x7f5a7yU8AhPHfUX29d9l/x/MMGzzd/BXAY+Ar4nuvMPjjm9cBZoKz112rXmbv6mG9bdxMBPCqmg99jA/wTsA/YA8x1ndkHxzwE2IxnxEwZMMV1Zi8c81LgNNCI5+z8L4GfAD+56fv8duufyR5v/73WO09FRIKMv16KERGR+6RiFxEJMip2EZEgo2IXEQkyKnYRkSCjYhcRCTIqdhGRIKNiFxEJMv8f/aQ/gQH+T1cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pt.plot(x_01, df(x_01), \"--\", color=\"gray\", label=\"$f$\")\n",
    "pt.plot(nodes, fderiv, \"or\")\n",
    "pt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Why don't we hit the values of the derivative exactly?\n",
    "* Do an accuracy study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8560489763862222\n"
     ]
    }
   ],
   "source": [
    "#clear\n",
    "print(np.max(np.abs(df(nodes) - fderiv)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*   Can we assign a meaning to the entries of the matrix $D=V'V^{-1}$?\n",
    "*   What happens to the entries of $D$ if we...\n",
    "    * change $h$?\n",
    "    * shift the nodes?\n",
    "* Using this, how would you construct methods for finding $f''$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -8.333,  16.   , -12.   ,   5.333,  -1.   ],\n",
       "       [ -1.   ,  -3.333,   6.   ,  -2.   ,   0.333],\n",
       "       [  0.333,  -2.667,   0.   ,   2.667,  -0.333],\n",
       "       [ -0.333,   2.   ,  -6.   ,   3.333,   1.   ],\n",
       "       [  1.   ,  -5.333,  12.   , -16.   ,   8.333]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(\n",
    "    Vprime @ la.inv(V)\n",
    ").round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.  , 0.25, 0.5 , 0.75, 1.  ])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}