ketch/NLS double relaxation.ipynb

## NLS double relaxation.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8b645e7c",
   "metadata": {},
   "source": [
    "This is an implementation of the FD algorithm for NLS from Herbst, Morris, and Mitchell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "33f9e058",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "07923e74",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 20 # # of grid points\n",
    "\n",
    "It = np.eye(2*m)\n",
    "It[0,0] = 0.5; It[1,1] = 0.5\n",
    "It[-1,-1] = 0.5; It[-2,-2] = 0.5\n",
    "Itinv = np.diag(1./np.diag(It))\n",
    "\n",
    "A = np.array([[0,1.],[-1.,0]])\n",
    "S = np.zeros((2*m,2*m))\n",
    "for j in range(1,m-1):\n",
    "    S[2*j:2*j+2,2*j:2*j+2] = -2*A[:,:]\n",
    "    S[2*j:2*j+2,2*j-2:2*j] = A[:,:]\n",
    "    S[2*j:2*j+2,2*j+2:2*j+4] = A[:,:]\n",
    "S[:2,:2] = -A[:,:]\n",
    "S[-2:,-2:] = -A[:,:]\n",
    "S[:2,2:4] = A\n",
    "S[-2:,-4:-2] = A\n",
    "\n",
    "IinvS = np.dot(Itinv,S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "f609e5c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5 0.  0.  ... 0.  0.  0. ]\n",
      " [0.  0.5 0.  ... 0.  0.  0. ]\n",
      " [0.  0.  1.  ... 0.  0.  0. ]\n",
      " ...\n",
      " [0.  0.  0.  ... 1.  0.  0. ]\n",
      " [0.  0.  0.  ... 0.  0.5 0. ]\n",
      " [0.  0.  0.  ... 0.  0.  0.5]]\n"
     ]
    }
   ],
   "source": [
    "print(It)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "1ba1ddf9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0. -1.  0.  1.  0.  0.]\n",
      " [ 1. -0. -1.  0.  0.  0.]\n",
      " [ 0.  1. -0. -2.  0.  1.]\n",
      " [-1.  0.  2. -0. -1.  0.]\n",
      " [ 0.  0.  0.  1. -0. -2.]\n",
      " [ 0.  0. -1.  0.  2. -0.]]\n"
     ]
    }
   ],
   "source": [
    "print(S[:6,:6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "bf5603d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Domain: [-20,20]\n",
    "dx = 1./8\n",
    "m = int(40/dx)\n",
    "x = np.linspace(-20,20,m)\n",
    "\n",
    "It = np.eye(2*m)\n",
    "It[0,0] = 0.5; It[1,1] = 0.5\n",
    "It[-1,-1] = 0.5; It[-2,-2] = 0.5\n",
    "Itinv = np.diag(1./np.diag(It))\n",
    "\n",
    "A = np.array([[0,1.],[-1.,0]])\n",
    "S = np.zeros((2*m,2*m))\n",
    "for j in range(1,m-1):\n",
    "    S[2*j:2*j+2,2*j:2*j+2] = -2*A[:,:]\n",
    "    S[2*j:2*j+2,2*j-2:2*j] = A[:,:]\n",
    "    S[2*j:2*j+2,2*j+2:2*j+4] = A[:,:]\n",
    "S[:2,:2] = -A[:,:]\n",
    "S[-2:,-2:] = -A[:,:]\n",
    "S[:2,2:4] = A\n",
    "S[-2:,-4:-2] = A\n",
    "\n",
    "# Initial condition\n",
    "U = np.zeros((2*m))\n",
    "for i, xc in enumerate(x):\n",
    "    U[2*i] = 1./np.cosh(xc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "379edf23",
   "metadata": {},
   "outputs": [],
   "source": [
    "It = np.eye(2*m)\n",
    "It[0,0] = 0.5; It[1,1] = 0.5\n",
    "It[-1,-1] = 0.5; It[-2,-2] = 0.5\n",
    "Itinv = np.diag(1./np.diag(It))\n",
    "\n",
    "A = np.array([[0,1.],[-1.,0]])\n",
    "S = np.zeros((2*m,2*m))\n",
    "for j in range(1,m-1):\n",
    "    S[2*j:2*j+2,2*j:2*j+2] = -2*A[:,:]\n",
    "    S[2*j:2*j+2,2*j-2:2*j] = A[:,:]\n",
    "    S[2*j:2*j+2,2*j+2:2*j+4] = A[:,:]\n",
    "S[:2,:2] = -A[:,:]\n",
    "S[-2:,-2:] = -A[:,:]\n",
    "S[:2,2:4] = A\n",
    "S[-2:,-4:-2] = A\n",
    "\n",
    "IinvS = np.dot(Itinv,S)\n",
    "\n",
    "\n",
    "def udot(u,dx,q,ret,F):\n",
    "    ret[:] = 0.\n",
    "    ret += 1./dx**2 * np.dot(IinvS,u)\n",
    "    for j in range(len(u)//2):\n",
    "        F[2*j:2*j+2] = (u[2*j]**2+u[2*j+1]**2) * np.dot(A,u[2*j:2*j+2])\n",
    "    ret += q*F\n",
    "    return -ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "e467c1ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc0c2f0b970>]"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgCUlEQVR4nO3da3Bb533n8e+fAMELeL/pQlGiLMlOpMSOHVqJk7bjre36ksReT9rE3nqadDP1i60z2U324mwy2azzYuumzUyzdZM6iTeXSeJx0zRVu/Labi5tGl8i2fFNsmXTsmSRkngXLyBBkMSzL4CjQDQvIAngHAC/zwxH4MEhzl8Y6KeH//Oc55hzDhERKX4VfhcgIiK5oUAXESkRCnQRkRKhQBcRKREKdBGREhH268BtbW2uu7vbr8OLiBSlp59+etg5177Uc74Fend3N4cPH/br8CIiRcnMTi73nFouIiIlQoEuIlIiFOgiIiVCgS4iUiIU6CIiJWLVQDezB8xs0MxeXOZ5M7Mvm1mvmT1vZlfkvkwREVlNNiP0bwI3rPD8jcCe9NedwFc2XpaIiKzVqoHunPsXYHSFXW4Bvu1SngSazGxLrgoUKaSnjo/wpUeP8cwbY36XIrJmueihdwKnMr7vS297EzO708wOm9nhoaGhHBxaJHecc3zqb57jyz/p5b//8AW/yxFZs4KeFHXO3e+c63HO9bS3L3nlqohvXjozSd/YDLs76nj57CSnRqf9LklkTXIR6P1AV8b329LbRIrKY0cHMIN7P3gpAI8eHfC5IpG1yUWgHwD+ID3b5d3AuHPuTA5eV6Sgnjg+zNs7G3nnjmZ2tUd5vHfY75JE1mTVxbnM7PvA1UCbmfUB/wOoBHDOfRU4CNwE9ALTwB/mq1iRfOodnOKat2wC4C2bG3jx9LjPFYmszaqB7py7fZXnHfDHOatIxAdjsQTDUwl2d9QBsLujjodfPEN8boHqypDP1YlkR1eKigC9Q1MA7N7060BPOjg+FPOzLJE1UaCLkGq3AOxuTwX6nnSwe0EvUgwU6CLAqwNT1FSG6GyqAWBnW5QK+3XQixQDBboIcHIkRndblIoKA6AqHKKzuYYTw2q5SPFQoIsA/edm6GyqvmDb1sYazozP+FSRyNop0EWAM+NxtqbbLZ7OphpOn4v7VJHI2inQpezFZucZn5ljS+OFgb6lqZqzE3EWks6nykTWRoEuZc9rq2xd3HJpqmEh6RianPWjLJE1U6BL2etPt1UWt1y2pkfs/efUR5fioECXsncmHdhbGt88Qgd0YlSKhgJdyt7p8ThmsKnhwkDfkm7BnNGJUSkSCnQpe2fOzdBRX0Vl6MJ/Dg3VldRVhTmtEboUCQW6lL2hqVk66quXfK6jvkonRaVoKNCl7A1NztJeX7Xkc20KdCkiCnQpe0OTs7TXLR3o7fVVDE0p0KU4KNClrC0kHSOxxLIj9PY6jdCleCjQpayNTSdYSDo6GpYfoU/G54nPLRS4MpG1U6BLWfNG3yu1XDL3EwkyBbqUtUEv0JdruaS3DyrQpQgo0KWsDa0W6HUaoUvxUKBLWfOCum2ZlkuH13LRTBcpAgp0KWtDk7PURkJEq8JLPt8SjWCmEboUBwW6lLWx6QStdZFlnw+HKmiqqWQ0pkCX4FOgS1kbiSVoqV0+0AGaoxFGY4kCVSSyfgp0KWtjsQTN0ZUDvTUaYWRKgS7Bp0CXsjaaxQi9JRphbFqBLsGnQJeyNjadoGWVEXpLtEotFykKCnQpW/G5BaYTC6u2XFqilYxNz5HUzaIl4BToUra8UXc2I/SFpGMiPleIskTWTYEuZcsL9OZVeuit6cAfUdtFAi6rQDezG8zsmJn1mtndSzy/3cx+ama/MrPnzeym3Jcqklveic6V5qED51sy6qNL0K0a6GYWAu4DbgT2Areb2d5Fu30WeMg5dzlwG/BXuS5UJNfWOkJXoEvQZTNC3w/0OueOO+cSwIPALYv2cUBD+nEjcDp3JYrkR/Y9dAW6FIdsAr0TOJXxfV96W6bPA3eYWR9wEPj4Ui9kZnea2WEzOzw0NLSOckVyZyyWwAwaaypX3E+BLsUiVydFbwe+6ZzbBtwEfMfM3vTazrn7nXM9zrme9vb2HB1aZH1GpxM010YIVdiK+1VXhqiNhBToEnjZBHo/0JXx/bb0tkwfAx4CcM49AVQDbbkoUCRfRmMJmmtXHp17WrSeixSBbAL9ELDHzHaaWYTUSc8Di/Z5A7gGwMzeSirQ1VORQBuNrX6VqKc1GtG0RQm8VQPdOTcP3AU8ArxEajbLETO7x8xuTu/2KeCPzOw54PvAR51zuqxOAm0sNrfqDBdPasVFLaErwbb0qv6LOOcOkjrZmbntcxmPjwLvzW1pIvk1Op3gih1NWe3bEo3w6sBUfgsS2SBdKSplyTmXWjo3yxF6quWiEboEmwJdytJEfJ75pMu6h94cjRCfSzKdmM9zZSLrp0CXsjSW5VWiHl0tKsVAgS5laTS9jkvLKuu4eFqiVamfU6BLgCnQpSyNpm8pt9rdijwt0dR8dU1dlCBToEtZOj9Cz7KH7o3QxxToEmAKdClL53voWQe6eugSfAp0KUuj0wkioQqikVBW+zdUhwlVmAJdAk2BLmVpdCp12b/ZygtzecyM5trI+ZtiiASRAl3K0th0Iut2i6clWqkRugSaAl3KUmphruxWWvS0RCOMxXSjaAkuBbqUpbHp7Bfm8rREI+dnx4gEkQJdytJoLHH+6s9sNddGNG1RAk2BLmVnbiHJ+MzcOnroqZOiyaRWhpZgUqBL2Tk3neqDZ3tRkae5NkLSwURcfXQJJgW6lB1v6uF6euigy/8luBToUna8qYdr7qGn91cfXYJKgS5lZ3SNl/17vIW8NBddgkqBLmXHC+S19tC9pXZ1tagElQJdyo7XMmmqXeOFRedH6DopKsGkQJeyMzqdoL4qTFU4u4W5PDWRENWVFRqhS2Ap0KXsjMbWvo6Lp6U2oh66BJYCXcrORgK9OaqrRSW4FOhSdsamE7SssX/u0XouEmQKdCk7Y7G587eUWyut5yJBpkCXsjMSm13z0rmelmhEV4pKYCnQpazMJBaIzyXXf1I0GmEyPs/cQjLHlYlsnAJdyorX/25Z4zounvOX/6uPLgGkQJeyMrbOq0Q93n8EunORBJECXcrKyAYDvTnde9dcdAkiBbqUlbF1LszlaVHLRQIsq0A3sxvM7JiZ9ZrZ3cvs8yEzO2pmR8zse7ktUyQ3zo/Q19lD14qLEmTh1XYwsxBwH3Ad0AccMrMDzrmjGfvsAT4NvNc5N2ZmHfkqWGQjRmOzhCqMxpr1TVtsqtWa6BJc2YzQ9wO9zrnjzrkE8CBwy6J9/gi4zzk3BuCcG8xtmSK5MTKVoLk2QkWFrevnI+EK6qvCulpUAimbQO8ETmV835feluli4GIz+4WZPWlmNyz1QmZ2p5kdNrPDQ0ND66tYZANGYgna6tbXbvG01OlqUQmmXJ0UDQN7gKuB24GvmVnT4p2cc/c753qccz3t7e05OrRI9kZjiXXPcPE01+pqUQmmbAK9H+jK+H5belumPuCAc27OOfc68AqpgBcJlFwEeks0olkuEkjZBPohYI+Z7TSzCHAbcGDRPj8iNTrHzNpItWCO565MkdwYnpqlrW59C3N5Ugt06cIiCZ5VA905Nw/cBTwCvAQ85Jw7Ymb3mNnN6d0eAUbM7CjwU+C/OOdG8lW0yHok5pNMxudzMEKv1LRFCaRVpy0COOcOAgcXbftcxmMHfDL9JRJIXptkwz30aISZuQVmEgvURNZ2GzuRfNKVolI2hqdmATY+y6VWV4tKMCnQpWyMnl/HZYM99KiuFpVgUqBL2Rjd4MJcnlYFugSUAl3KxvBUKoA33HJJB/pIbHbDNYnkkgJdyoa3jktD9frWcfG01adaNsOTGqFLsCjQpWyMxja2jounvipMJFzBsEboEjAKdCkbw1MbX8cFwMxoi0Y0QpfAUaBL2cjFZf+etvqq89MgRYJCgS5lI6eBXqdAl+BRoEvZGMnBOi6e1mhEgS6Bo0CXspCYTzKRg3VcPG31VYxMJUiteiESDAp0KQu5WsfF01ZXxXzSMT6jVRclOBToUhZGcnRRkcd7HbVdJEgU6FIWvKs6N7qOi8frxQ9p6qIEiAJdykKuVlr0eIGuy/8lSBToUhYGJ1LB29FQnZPXO99ymVSgS3Ao0KUsDE7OUhsJUVeV1T1dVtVUG6HCfr3gl0gQKNClLAxNztJen5v+OUCowmiJ6uIiCRYFupSFwck4HTkMdEi1XTRClyBRoEtZyPUIHaBd67lIwCjQpSwMTs7SUZ+bE6IeXf4vQaNAl5IXn1tgMj6f8xG6t0CXLv+XoFCgS8kbSk8tzHmg11cRn0sSSyzk9HVF1kuBLiVvcDIOkPOTot7NokfUdpGAUKBLycvXCN17vUFdXCQBoUCXkucFbq5Pim5uTL3ewEQ8p68rsl4KdCl5gxOzVFjuls71bE4vI3B2XIEuwaBAl5I3NJm6U1GownL6uo01lUTCFRqhS2Ao0KXkDU7G6WjIbf8cwMzY3FDN2Qn10CUYFOhS8gYnZ2nP0b1EF9vcUM2AWi4SEAp0KXlDebhK1LOpsZqzarlIQGQV6GZ2g5kdM7NeM7t7hf0+aGbOzHpyV6LI+i0kHcNTuV/HxbO5oYqBibiuFpVAWDXQzSwE3AfcCOwFbjezvUvsVw98Angq10WKrNdoLEHSkZceOsCmhmpm55O6WbQEQjYj9P1Ar3PuuHMuATwI3LLEfl8A7gX0+6cEhneVaL566Ju8qYtqu0gAZBPoncCpjO/70tvOM7MrgC7n3P9d6YXM7E4zO2xmh4eGhtZcrMhanTmXCtotTTV5eX3v4iLNRZcg2PBJUTOrAL4EfGq1fZ1z9zvnepxzPe3t7Rs9tMiqzozPALC1MT8nRb2LizQXXYIgm0DvB7oyvt+W3uapB94G/MzMTgDvBg7oxKgEQf+5OJUhoy1PLRevN392XHPRxX/ZBPohYI+Z7TSzCHAbcMB70jk37pxrc851O+e6gSeBm51zh/NSscganBmfYXNjNRU5vkrUUxUO0RKNMDCpEbr4b9VAd87NA3cBjwAvAQ85546Y2T1mdnO+CxTZiNPnZtjamJ/+uWeTLi6SgAhns5Nz7iBwcNG2zy2z79UbL0skN06fi7N/Z0tej7G5oUqzXCQQdKWolKyFpOPsRJytTfk5IerZ1FCtk6ISCAp0KVlDk7MsJB1bCtByGZ5KkJhP5vU4IqtRoEvJ6j+XmrLYmac56J4tutGFBIQCXUrW6XSgb8lzy2Vbcy0AfWMzeT2OyGoU6FKyzl9UlOcReldL6vVPjU3n9Tgiq1GgS8k6fS5OfVWYhurKvB5nS2MNZtA3qkAXfynQpWSdPjeT93YLQCRcwZaGak6p5SI+U6BLyTo9PpP3dotnW0stfWq5iM8U6FKyzpyL533KoqeruZZToxqhi78U6FKS4nMLjMQSdBag5QKwrbmGgck4s/MLBTmeyFIU6FKSvCmEhWq5dLXU4hz0q48uPlKgS0k6ORIDoLstWpDjdTWn/uPQXHTxkwJdStKJkdQJyu7WwgT6tpbUxUWaiy5+UqBLSToxHKOhOkxzbX7noHs2N1RTGTKdGBVfKdClJJ0YidHdFsUsPze2WCxUYWxtqtEIXXylQJeSdGIkxo4CtVs8Xc216qGLrxToUnIS80n6x2bY2Vpb0ONua67R5f/iKwW6lJy+sWmSjsKP0FtqGYklmIzPFfS4Ih4FupScEwWesujZ1V4HwGtDsYIeV8SjQJeSc2LYm7JY2JbL7o50oA9OFfS4Ih4FupScEyMx6qvDtEQjBT3ujtZawhVG75ACXfyhQJeSc2Jkmu7Wwk1Z9FSGKuhui9KrEbr4RIEuJedkeg66H3a316nlIr5RoEtJScwn6RubKXj/3LOrI8rJ0WkS80lfji/lTYEuJeX14RgLSXf+BGWh7e6oYyHpzi8OJlJICnQpKa8MTAKwp6Pel+Pvbk8dV3108YMCXUrKKwOThCqMi9r96aHv6kgdV4EuflCgS0l5ZWCSHa21VFeGfDl+bSRMZ1MNr2nqovhAgS4l5dWBKS72qd3iuag9qrno4gsFupSMmcQCJ0ZiXLzZ30C/eFM9vYNTzC9oposUVlaBbmY3mNkxM+s1s7uXeP6TZnbUzJ43sx+b2Y7clyqysmMDkyQd7N3S4Gsdb+tsID6X1JouUnCrBrqZhYD7gBuBvcDtZrZ30W6/Anqcc5cCPwD+NNeFiqzmyOlxAPZt9TfQ397ZBMDzfed8rUPKTzYj9P1Ar3PuuHMuATwI3JK5g3Pup845byHoJ4FtuS1TZHVHT0/QUB1mW/qGzX65qC1KNBLixf5xX+uQ8pNNoHcCpzK+70tvW87HgIeXesLM7jSzw2Z2eGhoKPsqRbJw5PQEe7c2FHwNl8UqKox9Wxt5XoEuBZbTk6JmdgfQA3xxqeedc/c753qccz3t7e25PLSUubmFJC+fnWDf1ka/SwHg7dsaOXp6QidGpaCyCfR+oCvj+23pbRcws2uBzwA3O+dmc1OeSHaOnZ0kPpfksq4mv0sB4O2djczOJ3lVFxhJAWUT6IeAPWa208wiwG3AgcwdzOxy4K9Jhflg7ssUWdmzp84B8I5tTb7W4XlbZ+o3hRfUdpECWjXQnXPzwF3AI8BLwEPOuSNmdo+Z3Zze7YtAHfA3ZvasmR1Y5uVE8uLZU+doiUboavH3hKhHJ0bFD+FsdnLOHQQOLtr2uYzH1+a4LpE1efbUOS7b1uj7CVFPRYWxr7OR5/oU6FI4ulJUit5YLEHv4BQ93S1+l3KBnh3NHOkfJzY773cpUiYU6FL0Dp0YBWD/zmAF+lW7WplPuvP1ieSbAl2K3i9fHyUSruDSbcGYsujp2dFCZch44rURv0uRMqFAl6L3yxOjvGNbE1Vhf5bMXU5NJMTlXc08cVyBLoWhQJeidm46wQv941y1q9XvUpb07l2tvNg/zvjMnN+lSBlQoEtRe/y1EZyD37q4ze9SlnTVRa0kXaotJJJvCnQpaj9/dZj6qjCXBeSCosUu395EVbhCfXQpCAW6FC3nHP98bJCrdrUSDgXzo1xdGeLK7hZ+9sogzjm/y5ESF8x/BSJZOHJ6gtPjca7du8nvUlZ0/b5NHB+KaV0XyTsFuhStx44OYAa//ZYOv0tZ0fX7NmMGD79w1u9SpMQp0KVoPXLkLO/c3kxbXZXfpayoo6Gad25v5uEXz/hdipQ4BboUpVcGJnn57CTvv3SL36Vk5Ya3bebls5O8Pqz7jEr+KNClKP3Dc6epMLipiAId0Chd8kqBLkVnIen44TP9vGdXGx311X6Xk5VtzbVcvr2JHzzdp9kukjcKdCk6/9o7TP+5GT58ZdfqOwfIHe/awfGhGI9rTrrkiQJdis73njpJc20lv7Mv2NMVF3vfpVtorq3k20+c8LsUKVEKdCkqJ0diPHp0gNv3bw/cYlyrqa4M8aEru3js6ABnxmf8LkdKkAJdiso3/vV1whXGR97T7Xcp63LHu3bggG89ftLvUqQEKdClaJwdj/PgoVPcenknmxqK42ToYl0ttXzg0q186/ETDE7E/S5HSowCXYrGX/z4VZJJx8d/e4/fpWzIJ6+7mLmFJF/+yat+lyIlRoEuReG5U+d48NAb/MFV3XS11PpdzoZ0t0W5ff92HvzlKV1oJDmlQJfAW0g6PvujF2mvq+I/XVfco3PPx6/ZTSRcwWd/9ALJpOalS24o0CXw/s8vXueF/nE+8763Ul9d6Xc5OdFRX81n37eXX/SO8J0ndYJUckOBLoH29MlR/uThl7n2rZu4+bKtfpeTU7fv7+LqS9r5Xw+/RO/gpN/lSAlQoEtgDU7G+Q/ffYatTTX8+Ycuw8z8LimnzIx7P3gp0UiYP/zmIYYmZ/0uSYqcAl0CaXhqlt//2lNMzMzzlTuuoLGmNFoti21qqOYbH72SoclZPvatQ0zNzvtdkhQxBboEztnxOP/ua09yamyaBz56Jfu2NvpdUl69o6uJv7z9Co6cnuDDf/2E5qfLuinQJVCeeG2E9//vn9M3NsMDH72Sq3a1+l1SQVy7dxNf/0gPrw/HuPWvHudXb4z5XZIUIQW6BMJkfI57/uEov//1J2moqeTv//i9vGdXm99lFdS/uaSDB+98NwAf/Mrj/Nkjx5hJLPhclRQT82tt5p6eHnf48GFfji3BEZud53tPvcH9Pz/O8NQst+/fzqdvfEvJTE9cj4n4HP/zwFH+9pk+Ouqr+Pg1e/jdK7ZREymuxcgkP8zsaedcz5LPKdCl0BLzSZ4+OcaB507zj8+fZjI+z1UXtfLfbnwL7+hq8ru8wDh0IjVl8+mTYzRUh7n18k6u37eZK3e2UBnSL9flasOBbmY3AH8BhICvO+f+ZNHzVcC3gXcCI8CHnXMnVnpNBXr5GJ6a5aUzExw9PcEvXx/lyeMjxBIL1FSGuPFtm7njqh1csb3Z7zIDyTnHoRNjfPuJEzx6dIDEfJKG6jDv2dXGO7Y3cdm2Ji7ZXE9zbWXJTeuUpa0U6OEsfjgE3AdcB/QBh8zsgHPuaMZuHwPGnHO7zew24F7gwxsvXYJgIemYW0gyO58kMZ9kbiH15+x8ksn4HOMzc0zE5xifnmN8Zp6x6QR9YzP0n5uhf2yaifivp+J1t9Zy6xWd/Oaedn5jdxvRqlU/gmXNzNi/s4X9O1uYTszz81eHeezoAIdOjPL/jpw9v199VZiullp2tNbSUV9FczRCSzRCc23qqyYSorqygurKEDWVofN/RsIVVBj6z6BEZPOvaT/Q65w7DmBmDwK3AJmBfgvw+fTjHwB/aWbm8tDPeejQKe7/+XGAN92b0S3zTeb2zJ9ZXJy74Gfc0ttX+Bst99rLve6bn1vuONn+TJZ/t+X2W/S6C0lHYiHJwhrXGqmvCrO1qYbO5hp6djSzo7WWvVsaeOuWBpqjkTW9lvxabSTM9fs2c/2+1A2nR2MJnu87x2tDMd4YiXFydJpjA5M8/toI4zNza3rtUIURMqOiAkJmqe/TXxVmhCtsydDP3OQ9NmyJbZn72Zu2scJ+pegT1+zhA3m48jmbQO8ETmV83we8a7l9nHPzZjYOtALDmTuZ2Z3AnQDbt29fV8HN0QiXbKrPeNELn1/uA3Hh9mV/fNmf4YKfufCnlns9W+ZnFn9ObZkDrVzn6q+90t8tm9cKh4zKkBEJpUZylSGjKlxBxPsKhaivDtNYU0lDTSWNNZXUV4fV3y2QlmiEqy/p4OpL3vzc/EKSczNzjMUSjE3PMTO3QDz9NZNI/zmX+k1rwTkWkkkWkpB0joVkxpdzJJOO+eSbBzhuiVHTUgOMpQYqq+33phFJicnXhXIF/X3XOXc/cD+keujreY3r9m7iur3FdS9JkUILhypoq6uira7K71KkgLIZSvUDmbdX35betuQ+ZhYGGkmdHBURkQLJJtAPAXvMbKeZRYDbgAOL9jkAfCT9+HeBn+Sjfy4iIstbteWS7onfBTxCatriA865I2Z2D3DYOXcA+AbwHTPrBUZJhb6IiBRQVj1059xB4OCibZ/LeBwHfi+3pYmIyFpoOoKISIlQoIuIlAgFuohIiVCgi4iUCN9WWzSzIWC9tztvY9FVqAGhutZGda1dUGtTXWuzkbp2OOfal3rCt0DfCDM7vNxqY35SXWujutYuqLWprrXJV11quYiIlAgFuohIiSjWQL/f7wKWobrWRnWtXVBrU11rk5e6irKHLiIib1asI3QREVlEgS4iUiKKKtDN7Itm9rKZPW9mf2dmTRnPfdrMes3smJldX+C6fs/MjphZ0sx6MrZ3m9mMmT2b/vpqEOpKP+fb+7Wojs+bWX/Ge3STX7Wk67kh/Z70mtndftaSycxOmNkL6ffIt7urm9kDZjZoZi9mbGsxs8fM7NX0nwW/4/cydfn+2TKzLjP7qZkdTf9b/ER6e37eM+dc0XwBvwOE04/vBe5NP94LPAdUATuB14BQAet6K3AJ8DOgJ2N7N/Cij+/XcnX5+n4tqvHzwH/2+7OVriWUfi8uAiLp92iv33WlazsBtAWgjt8Crsj8XAN/Ctydfny39+8yAHX5/tkCtgBXpB/XA6+k//3l5T0rqhG6c+5R55x3C/knSd09CVI3qX7QOTfrnHsd6CV1c+tC1fWSc+5YoY6XrRXq8vX9CrDzN0R3ziUA74bokuac+xdS9zzIdAvwrfTjbwH/tpA1wbJ1+c45d8Y590z68STwEql7MOflPSuqQF/k3wMPpx8vdSPrzoJXtLSdZvYrM/tnM/tNv4tJC9r7dVe6jfaAH7+uZwja+5LJAY+a2dPpm60HySbn3Jn047NAkG76G5TPFmbWDVwOPEWe3rOC3iQ6G2b2T8DmJZ76jHPu79P7fAaYB74bpLqWcAbY7pwbMbN3Aj8ys33OuQmf6yqolWoEvgJ8gVRgfQH4c1L/WcuFfsM5129mHcBjZvZyelQaKM45Z2ZBmQsdmM+WmdUBfwv8R+fchJmdfy6X71ngAt05d+1Kz5vZR4H3A9e4dAOK7G5knde6lvmZWWA2/fhpM3sNuBjI2Umt9dRFAd6vTNnWaGZfA/4xX3VkoaDvy1o45/rTfw6a2d+Rag8FJdAHzGyLc+6MmW0BBv0uCMA5N+A99vOzZWaVpML8u865H6Y35+U9K6qWi5ndAPxX4Gbn3HTGUweA28ysysx2AnuAX/pRYyYzazezUPrxRaTqOu5vVUCA3q/0h9lzK/DicvsWQDY3RC84M4uaWb33mNTkAD/fp8UybxL/ESAovxn6/tmy1FD8G8BLzrkvZTyVn/fMzzPA6zhj3Euqx/ls+uurGc99htQMhWPAjQWu61ZS/dZZYAB4JL39g8CRdK3PAB8IQl1+v1+LavwO8ALwfPpDvsXnz9hNpGYivEaqbeVbLRk1XURqxs1z6c+Tb3UB3yfVSpxLf7Y+BrQCPwZeBf4JaAlIXb5/toDfINXyeT4jt27K13umS/9FREpEUbVcRERkeQp0EZESoUAXESkRCnQRkRKhQBcRKREKdBGREqFAFxEpEf8fmeceMhWiNCgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,U[::2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "id": "9e9fc0ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc0c7752070>]"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEDCAYAAAAsr19QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg6ElEQVR4nO3deXxU9b3/8dcHkH3f1xBQ9k0xIri11g2pFavQq13Erdje/rrd21ZAAbXtrYJdb3trvd1stdaKQAAVQevS9losqCQTNgERCIGAQMKW/fP7Yw52TAFDZpKZM/N+Ph7zyJnvnJnzySF588l3zpxj7o6IiIRPk2QXICIi9aMAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkGr0ADezX5tZsZlFEvR61Wb2VnBbchrPG2pmr5lZuZl94xTrfczM3jCziJk9ambNgvFOZrbIzPLM7HUzGxmMtwzurzWzAjO7L+a1/hJT6y4zWxyMm5n9xMw2B683NuY508zs7eA2LWb8XDPLD57zEzOzYLyzma0M1l9pZp3quw0RSXHu3qg34BJgLBBJ0OsdrsM6204w1h04D/gu8I2TPK8JsAMYHNy/H7g9WJ4PzA2WhwIvBssGtA2WzwBWAeNP8NpPAzcHy5OA54LnjgdWBeOdga3B107BcqfgsdeDdS147tXB+DxgRrA8A3iwvtvQTTfdUvvW6B24u78K7I8dM7MzzWy5ma0JutShjVBHsbv/A6g8xWpdgAp33xTcXwncECwPB/4cvNYGINvMenjU4WCdM4LbBz4tZWbtgY8Bi4OhycDvguf+HehoZr2Aq4CV7r7f3Q8E258YPNbe3f/u7g78Drgu5rUeDZYfrTVe522ccueJSEpIlTnwR4Avu/u5wDeA/zmN57Y0s9Vm9nczuy7Bde0DmplZTnB/CtAvWF4LXA9gZuOA/kDf4H5TM3sLKCYajqtqve51RDv20uB+H6Kd/nE7g7FTje88wThAD3cvCpZ3Az3quQ0RSXHNkl2AmbUFLgCeCqZxAVoEj11PdNqitkJ3vypY7u/uhWY2EPizmeW7+xYz+xlwYbBO7yBQAZ5y9+/WpTZ3dzO7EfihmbUAVgDVwcMPAD8OXjcfePP4Y+5eDZxtZh2BRWY20t1j5/xvAn5ZlxriEdSvcyWIpKmkBzjRvwIOuvvZtR9w94XAwlM92d0Lg69bzexl4Bxgi7t/6fg6ZrbtRK9fF+7+GnBx8DpXAoOD8VLg1mDcgHeIzh/HPvegmb1EdEoiEqzbFRgHfDJm1UL+2dlDtJMvDG4frTX+cjDe9wTrA+wxs17uXhRMkRTXcxsikuKSPoUSBOE7ZjYV3j9aYkxdnhscCXK8W+9KtONel8j6zKx78LUFcBfwcHC/o5k1D1a7A3jV3UvNrFvQeWNmrYArgA0xLzkFWObuZTFjS4Cbg+99PFASTIM8D1wZfJ+dgCuB54PHSs1sfPCfx81AbsxrHT+SZFqt8TpvI769JiKNorHfNQWeAIqIvnm4E7gdGAAsJzqvvA6YU8fXuoDo9MXa4OvtJ1lv2wnGegbbLwUOBsvtg8eeBXoHy/OB9cBG4Gsxz58AbArGF/LPo0NGE51OySPadc+ptd2XgYm1xgz4GbAl+D5yYh67Ddgc3G6NGc8JXn8L8FPAgvEuwIvA28ALQOf6bkM33XRL7dvxX3oREQmZpE+hiIhI/TTqm5hdu3b17OzsxtykiEjorVmzZp+7d6s93qgBnp2dzerVqxtzkyIioWdm755oXFMoIiIhpQAXEQkpBbiISEgpwEVEQkoBLiISUgpwEZGQUoCLiISUAlxEpAEVHjzGfUsLqKyuSfhrp8LpZEVE0o678+Q/dvCdZ9ZT48715/RlVN8OCd2GAlxEJMF2HTzGjIX5vLppLxMGdmHelNH069w64dtRgIuIJIi786fVO/jOsvVUu/PtySP4zPn9adLEPvzJ9aAAFxFJgNiue/zAzsyfMqZBuu5YCnARkTjU7rrvnzyCzzZg1x1LAS4iUk+7Dh5j5sJ8Xgm67nk3jCGrS8N23bEU4CIipym2666qadyuO5YCXETkNBSVHGPG09Gu+/wB0bnuxuy6YynARUTqwN15avVOvr1sHVU1zn3XjuBz4xu/646lABcR+RBFJdG57pc3RrvueVNG079Lm2SXpQAXETkZd+epNUHXXZ0aXXcsBbiIyAnsLilj5sI8Xtq4l3EDOjM/RbruWApwEZEYtbvuez8xnJsnZKdM1x1LAS4iEghD1x1LAS4iGc/dWbBmJ/eHoOuOpQAXkYz2ga47uzPzp6Z21x1LAS4iGcndefqNwvcvtjD3E8OZFoKuO5YCXEQyzu6SMmYtyufPG4oZlx09rju7azi67lgKcBHJGMe77vuXFlBRXcOca4ZzywXh6rpjKcBFJCPsKS1j5sJo131edifmTxkTyq47lgJcRNKau7MwmOtOh6471ocGuJn9GrgGKHb3kcFYZ+BJIBvYBnzK3Q80XJkiIqdvT2kZsxbm82LQdc+bMoYBIe+6YzWpwzq/BSbWGpsBvOjug4AXg/siIikh2nXv5IofvMLftuxj9jXD+eP0CWkV3lCHDtzdXzWz7FrDk4GPBsuPAi8DdyWyMBGR+ojtunP6d2L+1PTqumPVdw68h7sXBcu7gR4JqkdEpF7cnUVvFnLvkuhc9+xgrrtpGsx1n0zcb2K6u5uZn+xxM5sOTAfIysqKd3MiIv+iuDR6XPcL69O/645V3wDfY2a93L3IzHoBxSdb0d0fAR4ByMnJOWnQi4icrtiuu7yqhns+PoxbLxyQ1l13rPoG+BJgGvBA8DU3YRWJiNRBtOuO8ML6PZzbvxPzp4xmYLe2yS6rUdXlMMIniL5h2dXMdgJziQb3n8zsduBd4FMNWaSIyHHuzuK3Crl3yTrKKqszruuOVZejUG46yUOXJbgWEZFTqt11z5symjMzrOuOpU9iikjKc3dy39rF3CUFGd91x1KAi0hKKz5Uxt2LIqxct4exWR2ZP3VMRnfdsRTgIpKSanfdd08axm0XqeuOpQAXkZSjrrtuFOAikjLcnSVro133sQp13R9GAS4iKaH4UBn3LIqwYt0ezsnqyEPquj+UAlxEkiq26z5aUc2sSUO5/aKB6rrrQAEuIkmz91A5dy/Kf7/rnj9lDGd1V9ddVwpwEWl0tbvumVcP5Y6L1XWfLgW4iDSqvYfKuWdxPs8XqOuOlwJcRBqFuu7EU4CLSIPbe6ic2YsjLC/Yzdn9OvLQ1NGc1b1dsssKPQW4iDQYd2dpXhFzcyMcqahmxtVD+by67oRRgItIg4jtusf068j31XUnnAJcRBLK3VmWV8ScmK77josG0Kxpk2SXlnYU4CKSMPsOR7vu5yLRrvuhKaMZ1ENdd0NRgItIQizL28XsxRGOlFdz18ShfP5idd0NTQEuInHZd7icObkRns3fzZi+HXho6hh13Y1EAS4i9bYsbxdzcgs4XFalrjsJFOAictrUdacGBbiInJZn8oqYnRvhcFkV35o4hOkXD1TXnSQKcBGpk/cOlzM76LpHB133YHXdSaUAF5EPFdt1f/OqIdx5ibruVKAAF5GTeu9wOXNyC3gmv0hddwpSgIvICT2bX8TsxREOqetOWQpwEfmA9w6XM2dJAc/kFTGqT7TrHtJTXXcqiivAzeyrwOcBA/7X3X+UiKJEJDmeyy/insURSssq1XWHQL0D3MxGEg3vcUAFsNzMlrn75kQVJyKNo3bX/Yep49V1h0A8HfgwYJW7HwUws1eA64F5iShMRBqHuu7wiifAI8B3zawLcAyYBKyuvZKZTQemA2RlZcWxORFJpP1HKpiTG2GZuu7QqneAu/t6M3sQWAEcAd4Cqk+w3iPAIwA5OTle3+2JSOLEdt3fuHIwd37kTM5Q1x06cb2J6e6/An4FYGb/BexMRFEi0jD2H6lg7pIClq7dxcg+7Xl86vkM7dk+2WVJPcV7FEp3dy82syyi89/jE1OWiCTa8ki06y45Vsl/XjGYL3xUXXfYxXsc+NPBHHgl8CV3Pxh/SSKSSLW77t/ffj7DeqnrTgfxTqFcnKhCRCTxlkd2c8/ifHXdaUqfxBRJQweCrnvJ2l2M6K2uO10pwEXSTGzX/R9XDOaL6rrTlgJcJE2o6848CnCRNPB8wW7uXhSh5FiFuu4MogAXCbEDRyq4d2kBuW/tYniv9vz+9nHqujOIAlwkpFYU7GbWoggHj1bw9csH8++XquvONApwkZCp3XX/7rZxDO+trjsTKcBFQiS26/7a5YP40qVnqevOYApwkRA4eLSCe5cUsFhdt8RQgIukOHXdcjIKcJEUdfBoBfctXceiNwsZ1qs9j952HiN6d0h2WZJCFOAiKWjluj3MWpTPgSMVfPWyaNfdvJm6bvkgBbhIComd6x7asx2/vVVdt5ycAlwkRdSe6/73j6rrllNTgIskWexx3ZrrltOhABdJouPnMNGnKaU+FOAiSbD/SHSu+/iZA3Vct9SHAlykkel83ZIoCnCRRqJrU0qiKcBFGoGuCC8NQQEu0oDeO1zO3CUFLMsrYmSf9jx2x/kM7amuWxJDAS7SQJ7NL2L24gilZZV848rB3PkRdd2SWApwkQR773A5c3ILeCa/iFF9OvCHqeMZ0rNdssuSNKQAF0mgZ/KKmJ0b4XBZFd+8agh3XjKQZuq6pYEowEUSYN/hcubkRng2fzej+3Zg/pQx6rqlwSnAReK0LG8Xc3ILOFxWxbcmDmH6xeq6pXHEFeBm9nXgDsCBfOBWdy9LRGEiqS626x7TtwPzp45hcA913dJ46h3gZtYH+Aow3N2PmdmfgBuB3yaoNpGU5O4syytiTm6EI+XV3DVxKJ+/eIC6bml08U6hNANamVkl0BrYFX9JIqlr76FyZi+OsLxgN2P6deShKaMZpK5bkqTeAe7uhWb2ELAdOAascPcVtdczs+nAdICsrKz6bk4kqdydpXlFzM2NcKSimhlXD+WOi9R1S3LV+6fPzDoBk4EBQG+gjZl9tvZ67v6Iu+e4e063bt3qX6lIkhQfKuMLj63hK0+8Sf8ubXj2KxfxhY+cqfCWpItnCuVy4B133wtgZguBC4DHElGYSLK5O0vW7mLukgKOVlQz8+qh3HHxQJo2sWSXJgLEF+DbgfFm1proFMplwOqEVCWSZMWHyrhnUYQV6/ZwTlZH5k8Zw1nd2ya7LJEPiGcOfJWZLQDeAKqAN4FHElWYSDLU7rpnTRrK7Rep65bUFNdRKO4+F5iboFpEkqq4tIy7F0dYuW4PY7M6Mk9dt6Q4fRJTMp67s/itQu5dso6yymrunjSM2y4aoK5bUp4CXDJacWkZsxZFeGH9Hs7t34l5U0ZzZjd13RIOCnDJSO7OojcLuXdJAeVVNdzz8WHceqG6bgkXBbhknD2lZcxamM+LG4rJCbrugeq6JYQU4JIx3J2FbxRy39ICKqprmH3NcG65IFtdt4SWAlwyQlHJMWYtzOeljXvJ6d+J+VPHMKBrm2SXJRIXBbikNXfnyX/s4LvPrKeqxplzzXCmqeuWNKEAl7S1Y/9RZi7M56+b9zF+YGcevGE0/buo65b0oQCXtFNT4zy26l0eeG4DBnznupF8elwWTdR1S5pRgEta2bbvCHc9nceqd/Zz8aCufO/6UfTt1DrZZYk0CAW4pIXqGuc3f3uHh1Zs5IymTZh3w2im5vTFTF23pC8FuITe5uLDfGvBWt7YfpDLhnbnu58cRc8OLZNdlkiDU4BLaFVV1/C/f3mHH76widbNm/KjfzubyWf3VtctGUMBLqG0YXcp31qQR97OEiaO6Mn9142gezt13ZJZFOASKpXVNfzPS1v46Utv077lGfzs02P5+OheyS5LJCkU4BIakcISvrkgj/VFpVw7pjf3XjuCzm2aJ7sskaRRgEvKK6+q5r9f3MzPX9lClzbNeeRz53LliJ7JLksk6RTgktLe2nGQbz61lreLD3PD2L7MuWY4HVqfkeyyRFKCAlxSUlllNT9cuYn//ctWerRvyW9uPY9Lh3RPdlkiKUUBLiln9bb9fGtBHlv3HeGmcVnMnDSU9i3VdYvUpgCXlHG0oop5yzfy6Gvb6NOxFY/fcT4XntU12WWJpCwFuKSE/9uyj7uezmPH/mPcckE237xqCG1a6MdT5FT0GyJJdaiskgee28Djq7aT3aU1f7pzAuMGdE52WSKhoACXpHll015mPp3H7tIyPn/xAP7jiiG0at402WWJhIYCXBpdybFKvvvMOv60eidndW/Lgi9ewNisTskuSyR0FODSqFYU7OaexRHeO1LBly49ky9/bBAtz1DXLVIf9Q5wMxsCPBkzNBCY4+4/ircoST97D5Vz79ICnskrYliv9vxq2nmM6tsh2WWJhFq9A9zdNwJnA5hZU6AQWJSYsiRduDuL3izk/mXrOFpezTevGsL0SwZyRtMmyS5NJPQSNYVyGbDF3d9N0OtJGig8eIy7F+Xz8sa9jM3qyLwpozmre7tklyWSNhIV4DcCT5zoATObDkwHyMrKStDmJJXV1DiPBxcVduDeTwzncxOyaaqLCosklLl7fC9g1hzYBYxw9z2nWjcnJ8dXr14d1/YktW3de5gZT+fz+rboRYX/65Oj6NdZFxUWiYeZrXH3nNrjiejArwbe+LDwlvQWe3mzls2aMH/KaKacq4sKizSkRAT4TZxk+kQyQ8GuEu56Oo9IYakubybSiOIKcDNrA1wB3JmYciRMyiqr+e8/v83Dr2ylU+vm/PwzY7l6lC5vJtJY4gpwdz8CdElQLRIia96NnvJ1y94j3DC2L7OvGUbH1rq8mUhj0icx5bQcKa9i/vPRU7727tCKR28bx0cGd0t2WSIZSQEudfbqpr3MXJjPrpJjTJugU76KJJt+++RDHTxawXeeWc+CNTsZ2K0NT905gZxsnfJVJNkU4HJKz+UXMTu3gANHdfIpkVSjAJcTKi4tY05uAcsLdjOid3seve08RvTWyadEUokCXD7A3VmwZiffXraOsqoa7po4lM9fPIBmOvmUSMpRgMv7tr93lFmL8vnr5n2cl92JB24YzZnd2ia7LBE5CQW4UFVdw6//9g4/WLmJZk2a8O3JI/jM+f1popNPiaQ0BXiGixSWMGNh9GPwVwzvwf2TR9CrQ6tklyUidaAAz1DHKqr50Qub+OVf36Fzm+jH4CeO7KmTT4mEiAI8A/317X3MWpTP9v1HuWlcP2ZMHEaH1mckuywROU0K8Axy4Ej0AzlPv7GTAV3b8Mfp4xk/UKeyEQkrBXgGcHeWrN3F/UvXUXKsUh/IEUkTCvA0t/PAUe5ZHOHljXsZ068jj10/imG92ie7LBFJAAV4mqqucR79v208tGIjAHOuGc60C3RdSpF0ogBPQ+uLSpmxMJ+1Ow5y6ZBufPu6kfTtpOtSiqQbBXgaOX6FnF+8spUOrc7gxzeezbVjeuvQQJE0pQBPE3/f+h4zF+bzzr4jTDm3L3dPGkanNrpCjkg6U4CHXMnRSr733Hr++I8dZHVuzWO3n89Fg7omuywRaQQK8JByd56L7GbukgL2H6ngzksG8rXLB9OquQ4NFMkUCvAQ2nngKHNyC/jzhmJG9G7Pb245j5F9dK5ukUyjAA+RquoafvO3bfxg5SbM4J6PD+OWC7J1rm6RDKUAD4m3dhxk1sJ81hWVcvmw7tw3eSR9OuqsgSKZTAGe4g6VVfLQ8xv53d/fpXu7Fjz82bFcNUJnDRQRBXjKcneWR3Zz79ICig+VM21CNv955WDatdRZA0UkKq4AN7OOwC+BkYADt7n7awmoK6MVHjzG3NwIL6wvZniv9vziczmc3a9jsssSkRQTbwf+Y2C5u08xs+aAPq8dh6rqGn77f9E3Kd31JqWInFq9A9zMOgCXALcAuHsFUJGYsjLP2h0HmRm8SXnZ0O7cN3mEzl8iIqcUTwc+ANgL/MbMxgBrgK+6+5HYlcxsOjAdICsrK47NpadDZZV8f8UmHn1tG93attClzUSkzuL527wZMBb4ubufAxwBZtReyd0fcfccd8/p1q1bHJtLL8ffpLziB6/y6GvbuHl8f174z49w9aheCm8RqZN4OvCdwE53XxXcX8AJAlz+VfRNygJeWL+HoT3b8fPPjuWcrE7JLktEQqbeAe7uu81sh5kNcfeNwGXAusSVln5qv0k5a9JQbr1wAGfoTUoRqYd4j0L5MvB4cATKVuDW+EtKT29uP8A9iyMU7Crl0iHduH/ySPp11puUIlJ/cQW4u78F5CSmlPRUcrSSB5/fwBOvb6d7uxb87NNjmTRKb1KKSPz0ScwG4u48/UYh33t2PQePVXLbhQP4+hWDadtCu1xEEkNp0gA27TnEPYsjvP7OfsZmdeT3141ieG9dCV5EEksBnkBHK6r4yYub+eVfttK2ZTMeuH4Un8rpRxNdCV5EGoACPEFWFOzmvqXrKDx4jKnn9mXG1UPp0rZFsssSkTSmAI/Tjv1HuW9pAS+sL2ZIj3Y89YUJnJfdOdlliUgGUIDXU0VVDb/861Z+8uLbNDHTMd0i0ugU4PXw2pb3mJ0bYXPxYa4a0YO5nxhBb10dR0QamQL8NOw9VM73nl3PwjcL6de5Fb++JYePDe2R7LJEJEMpwOugusb5w+vbmb98A8cqq/l/l57Fly49i1bNmya7NBHJYArwDxEpLOHuRfms3VnCBWd24dvXjeTMbm2TXZaIiAL8ZEqOVvLQio08vupdOrdpwY9vPJtrx/TWR+BFJGUowGupqXGeWrODB5dv5ODRCm6ekM3XrxhMh1a6mLCIpBYFeIz8nSXMzo3w1o6DnJfdifuuPV8fgReRlKUABw4cqWD+io088fp2urRpwQ8+NYZPntNH0yUiktIyOsBrapwnV+9g3vINlJZVcesFA/jaFYNo31LTJSKS+jI2wNfuOMic3Ahrd5YwLrsz9183gqE9NV0iIuGRcQG+/0gF85/fwB//sYOubVvwo387m8ln6+gSEQmfjAnw6hrnide389CKjRwqq+L2Cwfw1csH0U7TJSISUhkR4G9uP8Cc3ALyC0sYP7Az908eyeAe7ZJdlohIXNI6wN87XM685Rt5cvUOerRvwU9uOodPjO6l6RIRSQtpGeDVNc4fVr3L/Oc3crSimumXDOQrlw3S9ShFJK2kXaKtefcAc3IjFOwq5YIzu3DftSMYpOkSEUlDaRPg+w6X88BzG1iwZic927fkp58+h4+P0nSJiKSv0Ad4VXUNj/39Xb6/chNlldV84SNn8uWPnUUbTZeISJoLdcqt3raf2bkFrC8q5aKzunLvtSM4q7tO9SoimSGUAb73UDnfe249C98opFeHlvzPZ8Zy9ciemi4RkYwSV4Cb2TbgEFANVLl7TiKKOpmq6hp+99q7/HDlJsqqqvniR6PTJa2bh/L/IRGRuCQi+S51930JeJ1TWrX1PeYuKWDD7kNcPCg6XaIr44hIJgtF6zpzYT5PvL6dPh1b8fBnx3LVCE2XiIjEG+AOrDAzB37h7o/UXsHMpgPTAbKysuq1kewurXUhYRGRWszd6/9ksz7uXmhm3YGVwJfd/dWTrZ+Tk+OrV6+u9/ZERDKRma050XuMTeJ5UXcvDL4WA4uAcfG8noiI1F29A9zM2phZu+PLwJVAJFGFiYjIqcUzB94DWBS8mdgM+IO7L09IVSIi8qHqHeDuvhUYk8BaRETkNMQ1By4iIsmjABcRCSkFuIhISCnARURCKq4P8pz2xsz2Au/W8+ldgQY/50oDCnP9Ya4dwl1/mGuHcNefSrX3d/dutQcbNcDjYWarG/pshw0pzPWHuXYId/1hrh3CXX8YatcUiohISCnARURCKkwB/i9nOgyZMNcf5toh3PWHuXYId/0pX3to5sBFROSDwtSBi4hIDAW4iEhIhSLAzWyimW00s81mNiPZ9dRmZv3M7CUzW2dmBWb21WC8s5mtNLO3g6+dgnEzs58E30+emY1N7ncAZtbUzN40s2XB/QFmtiqo8Ukzax6Mtwjubw4ez05q4dGaOprZAjPbYGbrzWxCyPb914Ofm4iZPWFmLVN1/5vZr82s2MwiMWOnva/NbFqw/ttmNi3J9c8PfnbyzGyRmXWMeWxmUP9GM7sqZjw1MsndU/oGNAW2AAOB5sBaYHiy66pVYy9gbLDcDtgEDAfmATOC8RnAg8HyJOA5wIDxwKoU+B7+A/gDsCy4/yfgxmD5YeCLwfK/Aw8HyzcCT6ZA7Y8CdwTLzYGOYdn3QB/gHaBVzH6/JVX3P3AJMBaIxIyd1r4GOgNbg6+dguVOSaz/SqBZsPxgTP3Dg7xpAQwIcqhpKmVS0n5wT2OHTwCej7k/E5iZ7Lo+pOZc4ApgI9ArGOsFbAyWfwHcFLP+++slqd6+wIvAx4BlwS/cvpgf6vf/DYDngQnBcrNgPUti7R2CALRa42HZ932AHUGYNQv2/1WpvP+B7FoBeFr7GriJ6DV0OdF6jV1/rcc+CTweLH8ga47v+1TKpDBMoRz/AT9uZzCWkoI/ac8BVgE93L0oeGg30YtgQOp9Tz8CvgXUBPe7AAfdvSq4H1vf+7UHj5cE6yfLAGAv8JtgCuiXwRWiQrHvPXpZwoeA7UAR0f25hvDsfzj9fZ1S/wa13Eb0rwYIQf1hCPDQMLO2wNPA19y9NPYxj/5XnXLHbJrZNUCxu69Jdi311Izon8Q/d/dzgCNE/4x/X6rue4Bgvngy0f+IegNtgIlJLSoOqbyvP4yZ3Q1UAY8nu5a6CkOAFwL9Yu73DcZSipmdQTS8H3f3hcHwHjPrFTzeCygOxlPpe7oQuNbMtgF/JDqN8mOgo5kdv2JTbH3v1x483gF4rzELrmUnsNPdVwX3FxAN9DDse4DLgXfcfa+7VwILif6bhGX/w+nv61T7N8DMbgGuAT4T/CcEIag/DAH+D2BQ8K58c6Jv3CxJck0fYGYG/ApY7+4/iHloCXD8HfZpROfGj4/fHLxLPx4oifkTtFG5+0x37+vu2UT37Z/d/TPAS8CUYLXatR//nqYE6yet43L33cAOMxsSDF0GrCME+z6wHRhvZq2Dn6Pj9Ydi/wdOd18/D1xpZp2Cv0CuDMaSwswmEp1CvNbdj8Y8tAS4MTjyZwAwCHidVMqkZEy81+NNh0lEj+zYAtyd7HpOUN9FRP9szAPeCm6TiM5Nvgi8DbwAdA7WN+BnwfeTD+Qk+3sI6voo/zwKZSDRH9bNwFNAi2C8ZXB/c/D4wBSo+2xgdbD/FxM9siE0+x64D9gARIDfEz3qISX3P/AE0bn6SqJ//dxen31NdK55c3C7Ncn1byY6p338d/fhmPXvDurfCFwdM54SmaSP0ouIhFQYplBEROQEFOAiIiGlABcRCSkFuIhISCnARURCSgEuIhJSCnARkZD6/zDfJ9mzBOPTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initial condition\n",
    "U = np.zeros((2*m))\n",
    "for i, xc in enumerate(x):\n",
    "    U[2*i] = 1./np.cosh(xc)\n",
    "    \n",
    "# Integrate with RK4\n",
    "\n",
    "dt = 0.05*dx*dx\n",
    "T = 1.\n",
    "t = 0.\n",
    "UU = []\n",
    "UU.append(U)\n",
    "E = []\n",
    "q = 0.5\n",
    "\n",
    "ret = np.zeros_like(U)\n",
    "F = np.zeros_like(U)\n",
    "\n",
    "\n",
    "while t<T:\n",
    "    f1 = udot(UU[-1],dx,q,ret,F)\n",
    "    y2 = UU[-1] + 0.5*dt * f1\n",
    "    f2 = udot(y2,dx,q,ret,F)\n",
    "    y3 = UU[-1] + 0.5*dt * f2\n",
    "    f3 = udot(y3,dx,q,ret,F)\n",
    "    y4 = UU[-1] + 0.5*dt * f3\n",
    "    f4 = udot(y4,dx,q,ret,F)\n",
    "    Unew = UU[-1] + dt/6 * (f1+2*f2+2*f3+f4)\n",
    "    UU.append(Unew)\n",
    "    t+=dt\n",
    "    E.append(dx*np.sum(UU[-1]**2))\n",
    "    \n",
    "plt.plot(E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "ef8d2bb1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc0c3ab5b50>]"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD4CAYAAAAQP7oXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuy0lEQVR4nO3deXxV5bX/8c+XQBgigwgqUyAggkwqBqItaqWKOIEUe+twkWpb5Xehva29ZXCoFrWV2lpba1VsudUWS1sGpU6gtVWplYJeCGGOGAREJpkDZFq/P84Tu42BHEjIOSdZ79eLF3uvPZy1N8le7L2f8zwyM5xzzrl4NEh0As4551KHFw3nnHNx86LhnHMubl40nHPOxc2LhnPOubg1THQCx1ObNm2sS5cuiU7DOedSyjvvvLPdzNpWtqxOF40uXbqwePHiRKfhnHMpRdL6wy3zx1POOefi5kXDOedc3LxoOOeci5sXDeecc3HzouGccy5uXjScc87FzYuGc865uHnRcM65OsTMmPGvD3h1xZbjsv+4ioakoZJWS8qXNLGS5WMkLZO0RNICSb1CfGCILZG0VNKIyDYFkW0+8w08Sd+VZJLahHlJ+kXIIVdS/2M/bOecq3vW79jPDb9eyMTZy3h2yabj8hlVfiNcUhrwKHAJsBFYJGmuma2IrPaMmT0e1h8GPAQMBfKAbDMrkdQOWCrpL2ZWEra7yMy2V/KZnYAhwAeR8GVA9/AnB3gs/O2cc/VaaZkxbcH7/PSV1TRq0IAfjujLtQM6HZfPiqcbkYFAvpmtA5A0AxgOfFI0zGxPZP0MwEK8MBJvUh6Pw8+A8cBzkdhw4GmLDTX4tqRWktqZ2eY49+mcc3XOqo/2MGFmLks37ubiM07mvqv7cmrLJsft8+IpGh2ADZH5jVTyP3xJY4HbgHRgcCSeA0wDOgOjIncZBsyXZMATZjY1rD8c2GRmSyVVlUcHwIuGc67eOVRSyqOv5fOrv79Hy6aNeOS6s7myXzsqXDdrXI11WGhmjwKPSroeuBMYHeILgd6SzgCekvSSmR0EBpnZJkknA69IWgUsBm4n9mjqmEi6BbgFIDMzs1rH5Jxzyeid9TuZMCuX/K37+NLZHbjryl6cmJFeK58dT9HYBEQfjnUMscOZQex9w6eY2UpJ+4A+wGIz2xTiWyXNIfYYbCeQRezdR/lnvStpYLx5hDuWqQDZ2dnxPg5zzrmkt/9QCQ/OW81T/yygfcum/PamAXyhx8m1mkM8RWMR0F1SFrGL9LXA9dEVJHU3s7Vh9gpgbYhnARvCi/DOQE+gQFIG0MDM9obpIcBkM1sGnBzZbwGxF+nbJc0FxoV3KjnAbn+f4ZyrL95Ys41Js5fx4e4D3HhuZ743tCcnNK790S2q/MRwwR8HzAPSgGlmtlzSZGJ3DOUX84uBYmJ3C6PD5oOAiZKKgTLgv0IB6ArMCXcTDYm1vnq5ilReBC4H8oFC4KajPFbnnEs5uwqLuPf5lcx6dyPd2mbw51vPI7tL64Tlo1hjpLopOzvbfBAm51wqMjNeXPYRd8/NY1dhMWMu7Ma4wafRpFHacf9sSe+YWXZly+r0yH3OOZeKtuw5yF3P5jF/xRb6dmjJ0zfn0Kt9i0SnBXjRcM65pGFm/HHRBu5/cSVFJWVMuqwnXxuURcO05OnxyYuGc84lgfU79jNp9jLeem8HOVmtmTKyH13aZCQ6rc/wouGccwl0uC5AGjQ4vl/SO1ZeNJxzLkFquwuQmuBFwznnalmiugCpCV40nHOuFiWyC5Ca4EXDOedqQTJ0AVITvGg459xxlixdgNSE1MzaOedSQLJ1AVITvGg451wNq9gFyLiLTqu1LkCONy8azjlXg5K5C5Ca4EXDOedqQCp0AVITvGg451w1pUoXIDXBi4Zzzh2jVOsCpCZ40XDOuWOQil2A1IS4HrZJGipptaR8SRMrWT5G0jJJSyQtkNQrxAeG2BJJSyWNiGxTENlmcSR+r6TcEJ8vqX2It5T0l7Cf5ZJ85D7nXK07VFLKQ/NXc+UvFrBx5wEeue5snrwxu14UDIhj5D5JacAa4BJgI7Exw68zsxWRdVqY2Z4wPYzYsK5DJTUDisKQse2ApUD7MF9AGP+7wudF9/UtoJeZjZF0O9DSzCZIagusBk41s6LD5e4j9znnalKqdwESr+qO3DcQyDezdWFnM4DhwCdFo/wiH2QAFuKFkXiT8viRHG5f4e/mivXodQLwMVASR/7OOVctdaULkJoQT9HoAGyIzG8EciquJGkscBuQDgyOxHOAaUBnYJSZlV/oDZgvyYAnzGxqZJv7gRuB3cBFIfxLYC7wIdAc+IqZlVWSxy3ALQCZmZlxHJ5zzh1eXeoCpCbUWANiM3vUzLoBE4A7I/GFZtYbGABMklT+4G+QmfUHLgPGSrogss0dZtYJmA6MC+FLgSVAe+As4JeSPvONGTObambZZpbdtm3bmjo851w9s6uwiO/+aSk3TvsXTRo14M+3nscPhvep1wUD4isam4BOkfmOIXY4M4CrKwbNbCWwD+gT5jeFv7cCc4g9BqtoOjAyTN8EzLaYfOB9oGcc+TvnXNzMjBdyN3PxQ6/z3JJNjLvoNF741vkp32dUTYmnaCwCukvKkpQOXEvsMdEnJHWPzF4BrA3xLEkNw3RnYhf5AkkZkpqHeAYwBMirZF/DgVVh+gPgi2GdU4AewLr4D9U5545sy56D3Pq7dxj7zLu0a9mUueMG8T+X9qgTfUbVlCrvs0JLp3HAPCANmGZmyyVNBhab2VxgnKSLgWJgJzA6bD4ImCipGCgj1qpqu6SuwJwwSlVD4Bkzezls84CkHmH99cCYEL8X+K2kZYCACRVbXjnn3LGoL12A1IQqm9ymMm9y65yrSsUuQB4Y2Y+sOtoFSLyq2+TWOefqnPrYBUhN8KLhnKt36msXIDXBi4Zzrt44VFLKo6/l86u/v0fLpo145LqzubJfO8L7VRcHLxrOuXqhvnQBcrx50XDO1WneBUjN8qLhnKuzvAuQmudnzzlX5+wqLOLe51cy692NdGubwZ9vPc+/0V1DvGg45+oMM+PFZR9x99w8dhUWM+6i0xg3+DT/RncN8qLhnKsTtuw5yF3P5jF/xRb6dmjJ0zfn0Kv9Z/o0ddXkRcM5l9K8C5Da5UXDOZeyvAuQ2udFwzmXcrwLkMTxouGcSyneBUhiedFwzqUE7wIkOXjRcM4lPe8CJHnE1bxA0lBJqyXlS5pYyfIxkpZJWiJpgaReIT4wxJZIWippRGSbgsg2iyPxeyXlhvh8Se0jy74Q4sslvV69Q3fOJbv9h0q4Z+5yrnn8LQ4UlfLbmwbw0FfO8oKRQFUOwiQpDVgDXAJsJDb863VmtiKyTgsz2xOmhxEboW+opGZAURj9rx2wFGgf5guA7Iqj71XY17eAXmY2RlIr4C1gqJl9IOnkML74YfkgTM6lLu8CJHGqOwjTQCDfzNaFnc0gNnb3J0Wj/CIfZAAW4oWReJPy+JEcbl/A9cBsM/sgrHfEguGcS03eBUhyi6dodAA2ROY3AjkVV5I0FrgNSAcGR+I5wDSgMzDKzErCIgPmSzLgCTObGtnmfuBGYDdwUQifDjSS9HegOfBzM3u6kjxuAW4ByMzMjOPwnHPJwLsASQ019pVJM3vUzLoBE4A7I/GFZtYbGABMklTeNm6QmfUHLgPGSrogss0dZtYJmA6MC+GGwDnAFcClwF2STq8kj6lmlm1m2W3btq2pw3POHUdb9hzk1t+9w9hn3qVdy6bMHTeI/7m0hxeMJBTPncYmoFNkvmOIHc4M4LGKQTNbKWkf0AdYbGabQnyrpDnEHoO9UWGz6cCLwN3E7nB2mNl+YL+kN4Azib1vcc6lIDPjz4s3cu8LK7wLkBQRz7/MIqC7pCxJ6cC1wNzoCpK6R2avANaGeJakhmG6M9ATKJCUIal5iGcAQ4C8SvY1HFgVpp8DBklqGF6w5wArj+ZgnXPJY8PHhYz6zb8YPyuXM9q14OVvX8CtF3bzgpHkqrzTCC2dxgHzgDRgmpktlzSZ2B3DXGCcpIuBYmAnMDpsPgiYKKkYKCPWqmq7pK7AnPClnIbAM2b2ctjmAUk9wvrrgTEhj5WSXgZyw7Jfm1leDZwD51wtKi0znv5nAT9+eTVpDcR9V/fh+oGZ3gVIiqiyyW0q8ya3ziWX/K17GT8zl3c/2MUXerTlhyP60r5V00Sn5SqobpNb55yrluLSMp54/T1+8dd8mjVO42dfOZOrz+rgXYCkIC8azrnjKm/TbsbPzGXF5j1c0a8dPxjWmzYnNE50Wu4YedFwzh0XB4tL+flf1zL1jXW0zkjniVHncGnvUxOdlqsmLxrOuRq3qOBjJszMZd32/fxHdkfuuLwXLZs1SnRargZ40XDO1Zh9h0p48OVVPP32ejq0asrvvjaQ87v7l2zrEi8azrkaEe1gcPR5XfjepT3I8A4G6xz/F3XOVcuuwiLue2ElM9+JdTA4c8x5nNPZOxisq7xoOOeO2ct5m7nz2eXsLCxi7EXd+Obg7t5fVB3nRcM5d9S27j3I3c8t56W8j+jdvgVP3TyA3u1bJjotVwu8aDjn4mZmzHp3E/c+v4IDxaWMH9qDb5zflUbeX1S94UXDOReXTbsOcPvsZby+ZhvZnU/kgZH9OO3kExKdlqtlXjScc0dUVmb8fuF6pry0CgN+MKw3o87t7B0M1lNeNJxzh7Vu2z4mzMplUcFOzu/ehh+O6Eun1s0SnZZLIC8azrnPKCkt48k33+dnr66hScMGPHhNP645p6N3MOi8aDjnPm3Fh3sYP2speZv2MLT3qUy+ujcnN29S9YauXoiryYOkoZJWS8qXNLGS5WMkLZO0RNICSb1CfGCILZG0VNKIyDYFkW0WR+L3SsoN8fmS2lf4rAGSSiRdc+yH7Zyr6FBJKT+dv5phv1zAR7sP8qsb+vP4qHO8YLhPqXIQJklpxMbhvoTYON2LgOvMbEVknRZmtidMDyM2Qt/QMCxrURj9rx2wFGgf5guAbDPbXuHzovv6FtDLzMZEcnkFOEhsBMGZR8rdB2FyLj7vrN/JhFm55G/dx5f6d+CuK3pxYkZ6otNyCVLdQZgGAvlmti7sbAaxsbs/KRrlF/kgA7AQL4zEm5THj+Rw+wq+CcwCBsSRt3OuCoVFJTw4bzW/fauAdi2a8L83DeCiHicnOi2XxOIpGh2ADZH5jUBOxZUkjQVuA9KBwZF4DjAN6AyMMrOSsMiA+ZIMeMLMpka2uR+4EdgNXBRiHYARYd6LhnPV9I/87UycncuGjw9w43mdGT+0Jyd4B4OuCjX2NU4ze9TMugETgDsj8YVm1pvYhX6SpPIHpIPMrD9wGTBW0gWRbe4ws07AdGBcCD8MTDCzsiPlIekWSYslLd62bVtNHZ5zdcbuA8VMnJXLDb9eSMMGDfjTrecxeXgfLxguLvH8lGwCOkXmO4bY4cwAHqsYNLOVkvYBfYDFZrYpxLdKmkPsMdgbFTabDrwI3A1kAzNCk782wOWSSszs2QqfMxWYCrF3GnEcn3P1xisrtnDns8vYtvcQt17Yle9cfLp3MOiOSjxFYxHQXVIWsWJxLXB9dAVJ3c1sbZi9Algb4lnAhvDiuzPQEyiQlAE0MLO9YXoIMLmSfQ0HVgGYWVbk834LPF+xYDjnKrd93yHumbuc53M30/PU5jx5Yzb9OrZKdFouBVVZNMIFfxwwD0gj1mppuaTJxO4Y5gLjJF0MFAM7gdFh80HAREnFQBmxVlXbJXUF5oS7hobAM2b2ctjmAUk9wvrrgTE1dbDO1TdmxnNLPuQHf1nO/kOlfPeS07n1wm6kN/QOBt2xqbLJbSrzJreuPtu8+wB3zMnjtVVbOTuzFT8e2Y/upzRPdFouBVS3ya1zLoWUlRl/WPQBP3pxFaVlxl1X9uKrn+tCmncw6GqAFw3n6pCC7fuZODuXt9d9zOe6ncQDX+pH5knewaCrOV40nKsDSsuMaQve56evrKZRgwY88KW+fGVAJ+9g0NU4LxrOpbjVH+1l/MylLN24m4vPOIX7ru7DqS29vyh3fHjRcC5FFZWU8au/5/Po3/Jp0aQRj1x3Nlf2a+d3F+648qLhXApaumEX42fmsnrLXoaf1Z67r+pNa+9g0NUCLxrOpZADRaU89MpqfrPgfU5u3oTfjM7mi2eckui0XD3iRcO5FPH2uh1MnJVLwY5Crs/JZOJlPWnRpFGi03L1jBcN55Lc3oPFPPDSKqYv/IDOJzXjmW/k8LlubRKdlqunvGg4l8T+tmort89ZxpY9B/n6oCy+O6QHTdO9g0GXOF40nEtCH+8vYvJflvPskg85/ZQT+NUNn+PszBMTnZZzXjScSyZmxvO5m7ln7nJ2Hyjmv7/YnbEXneYdDLqk4UXDuSSxZc9B7nw2j1dWbKFfx5ZM/0YOPU9tkei0nPsULxrOJZiZ8afFG7jvhZUUlZRx++U9ufnzWTRM87sLl3y8aDiXQBs+LmTi7Fz+kb+DgVmtmTKyH1ltMhKdlnOH5UXDuQQoLTOeequAB+etJq2BuO/qPlw/MJMG3n25S3Jx3f9KGipptaR8SRMrWT5G0jJJSyQtkNQrxAeG2BJJSyWNiGxTENlmcSR+r6TcEJ8vqX2I3xDiyyS9JenM6h++c7Uvf+tevvz4W0x+fgXndm3N/O9cwH+e29kLhksJVY7cJykNWANcAmwkNmb4dWa2IrJOCzPbE6aHERvWdaikZkBRGDK2HbAUaB/mC4BsM9te4fOi+/oW0MvMxkj6HLDSzHZKugy4x8xyjpS7j9znkklxaRlPvP4ev/hrPhmN07j7qt4MP6u9dzDokk51R+4bCOSb2bqwsxnAcOCTolF+kQ8yAAvxwki8SXn8SI6wr7ci8beBjnHk7lxSyNu0m+/NzGXl5j1c0a8dPxjWmzYnNE50Ws4dtXiKRgdgQ2R+I/CZ/+FLGgvcBqQDgyPxHGAa0BkYZWYlYZEB8yUZ8ISZTY1scz9wI7AbuKiSnL4GvFRZspJuAW4ByMzMjOPwnDt+DhaX8vCra3nyzXWclJHOE6PO4dLepyY6LeeOWY216TOzR82sGzABuDMSX2hmvYEBwCRJ5aPDDDKz/sBlwFhJF0S2ucPMOgHTgXHRz5F0EbGiMeEweUw1s2wzy27btm1NHZ5zR21Rwcdc/vM3efz197imf0de+c6FXjBcyounaGwCOkXmO4bY4cwArq4YNLOVwD6gT5jfFP7eCswh9hisounAyPIZSf2AXwPDzWxHHLk7V+v2HSrh+8/l8eXH/0lRaRm//1oOU67pR8tm3iOtS33xPJ5aBHSXlEWsWFwLXB9dQVJ3M1sbZq8A1oZ4FrAhvPjuDPQECiRlAA3MbG+YHgJMrmRfw4FVIZ4JzCb2iGvNMR+xc8fR62u2cfvsZXy4+wA3fb4L/zOkBxmNvWW7qzuq/GkOF/xxwDwgDZhmZsslTQYWm9lcYJyki4FiYCcwOmw+CJgoqRgoI9aqarukrsCc0GqkIfCMmb0ctnlAUo+w/npgTIh/HzgJ+FXYruRwb/edq227Cou49/mVzHp3I93aZjBzzHmc07l1otNyrsZV2eQ2lXmTW1cbXlq2mbueW87OwiL+34XdGDf4NJo08u7LXeqqbpNb51wltu49yN3PLeelvI/o3b4FT908gN7tWyY6LeeOKy8azh0lM2PWu5u49/kVHCguZfzQHtxyflfvYNDVC140nDsKG3cWcvucPN5Ys43szicy5Zp+dGt7QqLTcq7WeNFwLg5lZcbvF65nykurMOAHw3ozyvuLcvWQFw3nqvDetn1MnJXLooKdnN+9DT/6Ul86ntgs0Wk5lxBeNJw7jJLSMp58831+9uoamjZK4ydfPpOR/Tt4B4OuXvOi4VwlVny4h/GzlpK3aQ9De5/K5Kt7c3LzJlVv6Fwd50XDuYhDJaU88td8Hn/9PVo1S+exG/pzWd92iU7LuaThRcO54J31O5kwK5f8rfsY2b8jd115Bq2apSc6LeeSihcNV+8VFpXw4LzV/PatAtq3bMpvbxrAF3qcnOi0nEtKXjRcvbZg7XYmzs5l484D3HheZ8YP7ckJ3sGgc4flvx2uXtp9oJgfvrCSPy7eQFabDP5063kMzPIOBp2rihcNV+/MX/4Rdz6bx479RYy5sBvfvri7dzDoXJy8aLh6Y/u+Q9wzdznP527mjHYt+M3oAfTt6B0MOnc0vGi4Os/MeG7Jh/zgL8vZf6iU715yOmO+0I1G3sGgc0ctrt8aSUMlrZaUL2liJcvHSFomaYmkBZJ6hfjAEFsiaamkEZFtCiLbLI7E75WUG+LzJbUPcUn6RcghV1L/6h++q+s+3HWArz21mG//cQld2mTwwrcG8c0vdveC4dwxqnIQJklpwBrgEmAjseFfrzOzFZF1WpjZnjA9jNgIfUMlNQOKwuh/7YClQPswXwBkm9n2Cp8X3de3gF5mNkbS5cA3gcuBHODnZpZzpNx9EKb6q6zM+MOiD/jRi6soLTO+d2kPRn+uC2newaBzVaruIEwDgXwzWxd2NoPY2N2fFI3yi3yQAViIF0biTcrjR3K4fYXPfNpiVe5tSa0ktTOzzXEcg6tHCrbvZ8KsXBa+/zGfP+0kfjSiH5kneQeDztWEeIpGB2BDZH4jsf/pf4qkscBtQDowOBLPAaYBnYFRZlYSFhkwX5IBT5jZ1Mg29wM3AruBi46QRwfgU0VD0i3ALQCZmZlxHJ6rK0pKy5j2j/f56fw1pDdswJSRffmP7E7ewaBzNajGHuya2aNm1g2YANwZiS80s97AAGCSpPJe3waZWX/gMmCspAsi29xhZp2A6cC4o8xjqpllm1l227Ztq3lULlWs+mgPIx97ix++uIrzu7fl1dsu5CsDMr1gOFfD4ikam4BOkfmOIXY4M4CrKwbNbCWwD+gT5jeFv7cCc4g9BqtoOjDyGPNw9UBRSRk/e2UNVz2ygI07D/DIdWfz5I3ncEoL75HWueMhnqKxCOguKUtSOnAtMDe6gqTukdkrgLUhniWpYZjuDPQECiRlSGoe4hnAECCvkn0NB1aF6bnAjaEV1bnAbn+fUb8t2bCLqx5ZwM//upYr+rbjldsu5Koz2/vdhXPHUZXvNEJLp3HAPCANmGZmyyVNBhab2VxgnKSLgWJgJzA6bD4ImCipGCgj1qpqu6SuwJzwy90QeMbMXg7bPCCpR1h/PTAmxF8k1nIqHygEbqrmsbsUdaColIdeWc1vFrzPyc2bMO2r2QzueUqi03KuXqiyyW0q8ya3dc8/39vBxNm5rN9RyPU5mUy8rCctmjRKdFrO1SnVbXLrXMLtOVjMAy+t4pmFH9D5pGb84Rvncl63kxKdlnP1jhcNl/ReW7WF22fnsXXvQb5xfha3XdKDpunewaBzieBFwyWtj/cXMfkvy3l2yYf0OKU5j486h7M6tUp0Ws7Va140XNIxM/6Su5l75i5n78Fi/vuL3Rl70WmkN/T+opxLNC8aLqls2XOQO+bk8erKLZzZsSVTrsmh56ktEp2Wcy7wouGSgpnxp8UbuO+FlRSVlHHH5Wdw86As72DQuSTjRcMl3Ac7Cpk0J5d/5O8gJ6s1U0b2o0ubjESn5ZyrhBcNlzClZcZv3yrgJ/NWk9ZA3D+iD9cNyKSB3104l7S8aLiEWLtlL+Nn5fJ/H+zioh5tuX9EX9q3aprotJxzVfCi4WpVcWkZj//9PR55LZ+Mxmk8/JWzGH6W9xflXKrwouFqzbKNu/nezKWs+mgvV/Zrxz3DetPmhMaJTss5dxS8aLjj7mBxKQ+/upYn31zHSRnpTB11DkN6n5rotJxzx8CLhjuu/vX+x0yclcu67fv5SnYnbr/iDFo29Q4GnUtVXjTccbHvUAlTXlrF795eT6fWTZn+9Rw+f1qbRKflnKsmLxquxr2+Zhu3z17Gh7sPcPPns/ifS0+nWbr/qDlXF8TVmY+koZJWS8qXNLGS5WMkLZO0RNICSb1CfGCILZG0VNKIyDYFkW0WR+IPSlolKVfSHEmtQryRpKfCNislTar20bsatauwiNv+tITR0/5F0/Q0Zo75HN+/qpcXDOfqkCp/myWlAY8ClwAbgUWS5prZishqz5jZ42H9YcBDwFBiQ7hmh9H/2gFLJf3FzErCdheZ2fYKH/kKMClsMwWYBEwAvgw0NrO+kpoBKyT9wcwKjvHYXQ16cdlmvv9cHrsKi/nm4NMYN/g0Gjf07sudq2vi+S/gQCDfzNYBSJpBbOzuT4qGme2JrJ8BWIgXRuJNyuNHYmbzI7NvA9eULwIywpjjTYEiYA8uobbuPcj3n13Oy8s/ok+HFjx180B6t2+Z6LScc8dJPEWjA7AhMr8RyKm4kqSxwG1AOjA4Es8BpgGdgVGRuwwD5ksy4Akzm1rJZ98M/DFMzyRWrDYDzYDvmNnHleRxC3ALQGZmZhyH546FmTHznY3c98JKDhSXMmFoT75xfhYN07z7cufqshp72GxmjwKPSroeuBMYHeILgd6SzgCekvSSmR0EBpnZJkknA69IWmVmb5TvT9IdQAkwPYQGAqVAe+BE4E1Jr5bfAUXymApMhdgY4TV1fO7fNu4s5PY5ebyxZhsDupzIAyP70a3tCYlOyzlXC+IpGpuATpH5jiF2ODOAxyoGzWylpH1AH2CxmW0K8a2S5hArCm8ASPoqcCXwRTMrv/BfD7xsZsXAVkn/ALKBdbhaUVZm/O7t9Ux5eRUCJg/vzX/mdPYOBp2rR+J5lrAI6C4pS1I6cC0wN7qCpO6R2SuAtSGeFd5BIKkz0BMokJQhqXmIZwBDiL00R9JQYDwwrMI7kQ8Ij73CNucCq47ucN2xem/bPv7jiX9y99zlZHdpzbzvXMCN53XxguFcPVPlnUZoxTQOmAekAdPMbLmkycTuGOYC4yRdDBQDOwmPpoBBwERJxUAZ8F9mtl1SV2BO6KSuIbHWVy+HbX4JNCb2yArgbTMbQ6wF1/9KWg4I+F8zy62Bc+COoKS0jKlvruPhV9fStFEaP/nymYzs38E7GHSuntK/n/7UPdnZ2bZ48eKqV3SVWv7hbibMyiVv0x4u63MqPxjem5ObN0l0Ws6540zSO2aWXdky/9aV+4yDxaU88tpaHn99HSc2S+exG/pzWd92iU7LOZcEvGi4T3ln/ceMn5nLe9v2c805HbnzijNo1Sw90Wk555KEFw0HwP5DJTw4bzVP/bOA9i2b8vTNA7ng9LaJTss5l2S8aDjeWLONSaGDwdHndeF7l/Ygo7H/aDjnPsuvDPXYrsIi7nthJTPf2UjXthn8+dbzyO7SOtFpOeeSmBeNeuqlZZu567nl7CwsYuxF3fjm4O40aeQdDDrnjsyLRj2zde9B7n5uOS/lfUTv9i146uYB3sGgcy5uXjTqico6GPz6+Vk08g4GnXNHwYtGPbDh40Jun7OMN9du9w4GnXPV4kWjDisrM57+ZwE/nrcaAfcO780N3sGgc64avGjUUflb9zJh1jLeWb+TC09vy/0j+tDxxGaJTss5l+K8aNQxxaVlTH1jHT9/dS3NGqfx0H+cyYizvYNB51zN8KJRh+Rt2s33ZuaycvMerujXjnuu6k3b5o0TnZZzrg7xolEHHCwu5eFX1/Lkm+tonZHOE6PO4dLepyY6LedcHeRFI8X96/2PmTgrl3Xb9/OV7E7cfvkZtGzWKNFpOefqKC8aKWrfoRKmvLSK3729no4nNuX3X8thUPc2iU7LOVfHxfXNLklDJa2WlC9pYiXLx0haJmmJpAWSeoX4wBBbImmppBGRbQoi2yyOxB+UtEpSrqQ5klpFlvWT9E9Jy8O29XJEoL+t3sqQh17n9wvXc/Pns5j/nQu8YDjnakWVI/dJSgPWAJcAG4mNGX6dma2IrNPCzPaE6WHEhnUdKqkZUBSGjG0HLAXah/kCINvMtlf4vCHAa2GdKQBmNiGMNf4uMMrMlko6CdhlZqWHy72ujdy3c38R9z6/gtn/t4nTTj6BKSP7cU7nExOdlnOujqnuyH0DgXwzWxd2NgMYDnxSNMoLRpABWIgXRuJNyuNHYmbzI7NvA9eE6SFArpktDevtiCP3OsHMeGHZZu5+bjm7DxTzrcGnMXbwaTRu6B0MOudqVzxFowOwITK/EcipuJKkscBtQDowOBLPAaYBnYndJZSERQbMl2TAE2Y2tZLPvhn4Y5g+HTBJ84C2wAwz+3EledwC3AKQmZkZx+Elty17DnLXs3nMX7GFvh1a8vuv53BGuxaJTss5V0/V2ItwM3sUeFTS9cCdwOgQXwj0lnQG8JSkl8zsIDDIzDZJOhl4RdIqM3ujfH+S7gBKgOmRXAcBA4BC4K/hFuqvFfKYCkyF2OOpmjq+2mZm/GnxBu57YSVFJWVMuqwnXxuURUPvYNA5l0DxFI1NQKfIfMcQO5wZwGMVg2a2UtI+oA+w2Mw2hfhWSXOIPQZ7A0DSV4ErgS/av1+6bATeKH8HIulFoD/wqaJRF3ywo5BJc3L5R/4OBma1ZsrIfmS1yUh0Ws45F1frqUVAd0lZktKBa4G50RUkdY/MXgGsDfGs8AIbSZ2BnkCBpAxJzUM8g9j7irwwPxQYDwyr8E5kHtBXUrOwzwuJvFepC0rLjN8seJ9LH36DpRt2c9/VfZjxjXO9YDjnkkaVdxqhFdM4YhftNGCamS2XNJnYHcNcYJyki4FiYCfh0RSxx0kTJRUDZcRaVW2X1BWYE/pDagg8Y2Yvh21+CTQm9sgK4G0zG2NmOyU9RKyIGfCimb1QEychGazZspfxM3NZsmEXF/Voy/0j+tK+VdNEp+Wcc59SZZPbVJYKTW6LSsp4/PX3eOS1tZzQuCH3DOvNsDPbeweDzrmEqW6TW3ecLN2wiwmzcln10V6uOrM991zVi5NO8A4GnXPJy4tGAhwoKuXhV9fw5JvraNu8MU/emM0lvU5JdFrOOVclLxq17J/v7WDS7FwKdhRy3cBOTLr8DFo08Q4GnXOpwYtGLdlzsJgHXlrFMws/ILN1M575eg6fO837i3LOpRYvGrXgtVVbuH12Hlv3HuQb52dx2yU9aJruXYA451KPF43jaMe+Q0x+fgXPLfmQHqc05/FR53BWp1aJTss5546ZF43jwMz4S+5m7pm7nL0Hi/n2xd35ry+cRnpD7wLEOZfavGjUsM27D3DXs3m8unIrZ3ZqxY9H9qPHqc0TnZZzztUILxo1pKzMmLFoAz96cSXFZWXcecUZ3PT5LNIa+Jf0nHN1hxeNGlCwfT8TZ+fy9rqPOa/rSTwwsi+dT/L+opxzdY8XjWooLTOmLXifn76ymkYNGvCjL/Xl2gGdvAsQ51yd5UXjGK3+aC/jZy5l6cbdXHzGydx3dV9ObVkvhyx3ztUjXjSO0qGSUn71t/f41d/zadGkEY9cdzZX9mvndxfOuXrBi8ZR+L8PdjJhVi5rtuzj6rPa8/2retM6Iz3RaTnnXK3xohGHwqISfjp/DdP+8T6ntmjCtK9mM7indzDonKt/4vq2maShklZLypc0sZLlYyQtk7RE0gJJvUJ8YIgtkbRU0ojINgWRbRZH4g9KWiUpV9IcSa0qfFampH2S/ueYj/oovJW/naEPv8lvFrzP9QMzmf+dC7xgOOfqrSoHYZKUBqwBLiE2Tvci4DozWxFZp4WZ7QnTw4iN0DdUUjOgKIz+1w5YCrQP8wVAdvmY35F9DQFeC+tMATCzCZHlM4mN3LfQzH5ypNyrMwjT7gPF/OjFlcxYtIEuJzXjgZH9OLfrSce0L+ecSyXVHYRpIJBvZuvCzmYAw4mMz11eMIIMYhd1Kozx3aQ8fiRmNj8y+zZwTfmMpKuB94H9ceR9zHI37uIbTy9m295D3HphV75z8ek0aeQdDDrnXDxFowOwITK/EcipuJKkscBtQDowOBLPAaYBnYFRZlYSFhkwX5IBT5jZ1Eo++2bgj2E/JwATiN3xHPbRlKRbgFsAMjMz4zi8z8ps3YzTT2nOkzdm069jq2Pah3PO1UU11oOemT1qZt2IXdjvjMQXmllvYAAwSVL5lxkGmVl/4DJgrKQLovuTdAdQAkwPoXuAn5nZvirymGpm2WaW3bZt22M6llbN0vnd13K8YDjnXAXx3GlsAjpF5juG2OHMAB6rGDSzlZL2AX2AxWa2KcS3SppD7DHYGwCSvgpcCXzR/v3SJQe4RtKPgVZAmaSDZvbLOI7BOedcDYjnTmMR0F1SlqR04FpgbnQFSd0js1cAa0M8S1LDMN0Z6AkUSMqQ1DzEM4AhQF6YHwqMB4ZF34mY2flm1sXMugAPAz/0guGcc7WryjuN0IppHDAPSAOmmdlySZOJ3THMBcZJuhgoBnYCo8Pmg4CJkoqBMmKtqrZL6grMCd+ibgg8Y2Yvh21+CTQGXgnL3zazMTV0vM4556qhyia3qaw6TW6dc66+OlKTWx9KzjnnXNy8aDjnnIubFw3nnHNx86LhnHMubnX6RbikbcD6auyiDbC9yrWSUyrnDqmdfyrnDqmdfyrnDsmTf2czq/Tb0XW6aFSXpMWHa0GQ7FI5d0jt/FM5d0jt/FM5d0iN/P3xlHPOubh50XDOORc3LxpHVlnPu6kilXOH1M4/lXOH1M4/lXOHFMjf32k455yLm99pOOeci5sXDeecc3HzolEJSUMlrZaUL2liovOpSFInSX+TtELSckn/HeKtJb0iaW34+8QQl6RfhOPJldQ/sUcQIylN0v9Jej7MZ0laGPL8Y+iKH0mNw3x+WN4lwXm3kjRT0ipJKyWdl0rnXtJ3ws9NnqQ/SGqSzOde0jRJWyXlRWJHfb4ljQ7rr5U0urLPqqXcHww/O7mS5khqFVk2KeS+WtKlkXjyXJPMzP9E/hDr/v09oCuxoWuXAr0SnVeFHNsB/cN0c2AN0Av4MTAxxCcCU8L05cBLgIBzgYWJPoaQ123AM8DzYf5PwLVh+nHg/4Xp/wIeD9PXAn9McN5PAV8P0+nEBgVLiXNPbPjm94GmkXP+1WQ+98AFQH8gLxI7qvMNtAbWhb9PDNMnJij3IUDDMD0lknuvcL1pDGSF61Basl2TEvbDm6x/gPOAeZH5ScCkROdVRc7PERs7fTXQLsTaAavD9BPAdZH1P1kvgTl3BP5KbDz558Mv+fbIL9Mn/w7ExnI5L0w3DOspQXm3DBddVYinxLkPRWNDuHg2DOf+0mQ/90CXChfeozrfwHXAE5H4p9arzdwrLBsBTA/Tn7rWlJ/7ZLsm+eOpzyr/pSq3McSSUnhccDawEDjFzDaHRR8Bp4TpZDymh4mN0FgW5k8CdplZSZiP5vhJ/mH57rB+ImQB24D/DY/Wfh1Gn0yJc2+xYZZ/AnwAbCZ2Lt8hNc591NGe76T6d4i4mdidEaRI7l40UpikE4BZwLfNbE90mcX+S5KU7aklXQlsNbN3Ep3LMWhI7HHDY2Z2NrCf2OORTyT5uT8RGE6s+LUHMoChCU2qmpL5fB+JpDuAEmB6onM5Gl40PmsT0Cky3zHEkoqkRsQKxnQzmx3CWyS1C8vbAVtDPNmO6fPAMEkFwAxij6h+DrRSGFOeT+f4Sf5heUtgR20mHLER2GhmC8P8TGJFJFXO/cXA+2a2zcyKgdnE/j1S4dxHHe35Tqp/B0lfBa4EbghFD1Ikdy8an7UI6B5ak6QTe/k3N8E5fYokAb8BVprZQ5FFc/n3+Oyjib3rKI/fGFqWnAvsjtza1zozm2RmHc2sC7Hz+5qZ3QD8DbgmrFYx//Ljuiasn5D/WZrZR8AGST1C6IvAClLk3BN7LHWupGbh56g8/6Q/9xUc7fmeBwyRdGK42xoSYrVO0lBij2aHmVlhZNFc4NrQYi0L6A78i2S7JiXqZUoy/yHWAmMNsRYLdyQ6n0ryG0TsdjwXWBL+XE7sWfNfgbXAq0DrsL6AR8PxLAOyE30MkWP5Av9uPdWV2C9JPvBnoHGINwnz+WF51wTnfBawOJz/Z4m1xkmZcw/8AFgF5AG/I9ZaJ2nPPfAHYu9fiond6X3tWM43sfcH+eHPTQnMPZ/YO4ry393HI+vfEXJfDVwWiSfNNcm7EXHOORc3fzzlnHMubl40nHPOxc2LhnPOubh50XDOORc3LxrOOefi5kXDOedc3LxoOOeci9v/B5l2/CT5HN+UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skip = 1\n",
    "Umag = np.zeros((len(UU),m))\n",
    "E1 = []\n",
    "E2 = []\n",
    "for j, U in enumerate(UU[::skip]):\n",
    "    V = U[::2]\n",
    "    W = U[1::2]\n",
    "    Umag[j,:] = np.sqrt(V**2+W**2)\n",
    "    E1.append(dx*np.sum(V**2+W**2))\n",
    "    E2.append(np.sum((np.diff(V)**2+np.diff(W)**2)/dx - 0.5*q*dx*(V[1:]**2+W[1:]**2)**2))\n",
    "\n",
    "plt.plot(E2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "id": "f5111944",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fc0d3898640>"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAojElEQVR4nO2df6xsV3XfP2vm3fcekBTzI7KobRXTWI3cKA3WE3FEhKK4TcCJYioRZFQVN7VktYWWlFaJaaSStoqU9EdokFIiN3ZjKsSPOkRYFSkxhiiqKpw8fhkbh/gFQmzL4CSAQ0v8fN89q3/MPnP3nNn7/L7z8/uR7psz++xzzt5n5n3PmrXXXtvcHSGEEPvBZN0NEEIIsTok+kIIsUdI9IUQYo+Q6AshxB4h0RdCiD3i1LobUMdpO+Nned66myHGwmz20qJqHFNm8T8W/U2Ot31iYCy8ugETcAOfHJ/MJ4uvmMMEzBwzZxK9TsyZEPbhZReO2+ngoUfuRgEUbhRus/fh1d2gIDQGrGDp1RwIr1b4/CaU27PX8Bdvu89O4on71vY+K4pvp/gmX/8zd/+O1L6NFv2zPI/vsxvW3QzRF5uEF6t9P8eLaPNYhGxiMJ1ip07BwSns4AAODuDMafz0AX76FH72FEdnphRnphydnXDpzISjs8alM8bRGTg6MzvX0Vnm74uzTnG6gLMFdvqIgzOXODg44uzpQ55zcMhzTh3y3FOHnJ4ccXZ6yKnJEVNzJkEqC4wjNy4VUwCeOTrg2WLKty4d8JeXDvjLwwOeefaAw8MphxdP4c9O4ZkJk2cnTJ4xpheZ/T3DfPvURWf6jHPqYsH0mdn9mFw8YnrxCHvmEvbsJezZQ7j4LBwe4oeHcHgJv3QJjo6W71vb+xzKqu/FdvJRv+fLuX1y74jtxR0rLVT3maUMM6s4GL99mFQOnJgv7JtYMfurqZc6TysSbZ+997DtssrFICT6YnOwybJVGhMJfCx8FotgRQ+tKpCbppcL4p5oa25/fA/qHgJ191PsJfpGiJOj6jLoc7wXs+PLc/Q+F50Fv2q5T235BKmy6nGNDHkQRffFC5/fsz7ItbMfSPTF9lJE1m5B5AKh8ksgc3yivLNgN5A8X6JooVrsqgrvKeLtTfu5IrYJib5YHWNZkFW/dsbVM0S/l/zz4f3EiqismO/L1e9KfJjFYxaw3OexfPuy7PcKib7YGrwougtdxqpuq8kp8U6JfF39FMk29NFw99l9EaIlEn2x+eT81F7MXR3mztwILzwr6n0M8Ni6H6NeTKo9s3j96BdLEf2CKf32VQb48sV+IdEXK2f0AcPY+i/S5UZ+slLbHw+pEMwuYZm561iuUkG2P72JBtd7D7CLrUaiL06WMa3PBWu3g2Btg7Z1bWM8oDtm3L5+Lew8jaJvZneZ2VNm9lBU9h/M7A/M7EEz+00zuyza9zYzu2BmXzCzH4nKXx3KLpjZ7aP3RGw0vazKVIx57L+uxqrPt6NTjCz4E3OmzP7GjvRZjNiJtzP9hMX7MT9Rd1tOVv/+0Obb8evAqytl9wHf7e7fA/wh8DYAM7sWuBn4m+GY/2JmUzObAr8CvAa4FnhDqCtEknkKgUCtKKWid0bQsIWZuDX++nKGbuq4XoTDa6N3qodU7k/1/glR0ij67v67wNcqZb/t7pfC208AV4btm4D3uftFd/8ScAF4Rfi74O5fdPdngfeFukI0UhU09yJt4VZYspzjv5bUReoMqbvUpkDrw4tidh/i08laFy0Yw6f/D4HfCttXAI9F+x4PZbnyJczsNjM7b2bnD7k4QvPERjHEZ1z15cfvZ+kuE9fLCGnNA6CaSbMkZe3X/QJInidz3dpcQfFks+rkrKGTteTD3zsGib6Z/SxwCXjPOM0Bd7/D3c+5+7kDzox1WrHLLES4dIvDz5GLyplG5dOanwy9kq0F5u2vPgikz2IEeqdWNrN/APwYcIP73Px6ArgqqnZlKKOmXOw6XgxP/OUFMEth3Dp654S8HaV1X/gJBL+1bfOSxT/CE0FW/17Q61trZq8Gfhr4cXf/VrTrXuBmMztjZlcD1wC/B/w+cI2ZXW1mp5kN9t47rOli7wmRLBZy0yxlqZzXG279pxKrxfvq9jfRlAZ6oX9KrSwG0mjpm9l7gR8EXmxmjwNvZxatcwa4z2aOy0+4+z9y94fN7APA55m5fd7k7kfhPG8GPsLMXLvL3R8+gf6IDcYLn0eVzLdTvwJskp116kVRvyLUgjDaOKkOArmUDIN+UaQGcTuI+yw1RcZCzyyesjQwrgHgvaJR9N39DYniO2vq/zzw84nyDwMf7tQ6IVK0FMTsGGvHGbht/fNd6ycHkdt6WGTxi55oRq7YCpLWaLRW7EmtKLW8OlZxvHZuKprnJAYS4v5lxjNkrYu2bPQauUIs0SDspe8+XjqxMX/9Gki6nfx4X+sxCFn7oiOy9MXqGbqyEyzHpqfELzE5q6/Y11nwE7y3hb80iJuN1a/siPrfy8pXVs69RaIvVsOYAlPOxnUnTq/cNY696wMgFZdfF6s/1nXn/Zovh1j2d8xkdnoA7Aty74jtII7Tr63njakW2oiuLeTdWdyOwzNjIzuuZy0uUruMYzwLtw1jzIUQe4G+JWJjWUoaFln11bwz88HO6jlys1vnxy1eo41Yd2HpfJ4IOG3w4acGqef9TyyqomRrog6Jvlgpnf3PKes1dm8kBLEU/8GpGKIEam2zbHZOulahPDQ50WwewRO5t5ZO0O2/tKJ+9g+JvlgbJyE4pTbXuk7i1+o2I6RGzp0ndc2ayCIrOsTtt0QiLyT6Yv2kBhG7+qdzcfoN/v2Feg1MrFjw508zsfq9zl/Xzj5zEDKzcYWQ6IuNZ8FHHfuyy9elcMaybjSwWo3d70BV2CcUTBKhQp0XRo/i8o/LaiKR4iRriXWG5csXbZDoi62hdE14nU/bK6JYe8KRGtaVpuvGa9/W9LG8D3LZiC5I9MXqiF0MY7sbCpICWXWr956cVXPgkMHb+LD8OISPm0t/KfJJrp99QqIvtpPcQuG+HL2TjoRpd5mmtW+rZX19/Mk2L83W9cV+C9EDib7YDsLA5IIro074qrvaDugGcpb7hCjhWo353Xmt3Lr3S/WPK8zvhyZmiZbomyI2muTgZNKNk3GBNLh86rS5r8um7rjaa+ceYrkFYlJ902CuaECiL1bOKBO05ifzhUHbuVsEWsXkt6Htqli9Vs/KtTW4qRbEvmmZyA7WvheuAeA9RaIv1kqt8NSJWFHURui0S0vcok6gLrHaFO+WeK2hamPbC69Ptpa5bxJ5ARJ9sSU0ui1qcu+Ur0Mn2k6DL39KwTT49oesjZtqY3p/8+QsuXVEW5RlU2wf1cHchVDQY4HMimmHB0BXv37b+tU8+kttTQl9nFYZ2s1FEKKCLH2xWnIx4W1jxaN6C5k2q/7unGW8aTqZas+S2C+OW3jX+Q51dRSjv3dI9MX6GCI4XePU60I2U+mOI7oujJ5vQ+Y6fR5EfeP0JfJ7j9w7YuOxieFHmZ1F1eUR59D3RkGvY3HxlFm+ndKXD4AXFNhybp6+fv7Q5nn7q/ocp1VOIL++aIMsfbHZlJEoNqkXtbrBzjoNTuwbayGV5HmyvzZqTtQwkGsTW7hPQtShb4hYC4PCB6uZNiMsMZC7sB/S2S0rDM2p32qClof2JPf5cpx+SSLDZhcUurnfyL0jthIviiDgno5iKYW9aVGVBtr68/vWh8VFU4C81R8NVvuYi6KLvaLR0jezu8zsKTN7KCp7oZndZ2aPhtcXhHIzs3ea2QUze9DMrouOuSXUf9TMbjmZ7oi9pcYFMnzZxCjfjjlTK5hWyiZdFlTp2sY+i6gIkaGNe+fXgVdXym4H7nf3a4D7w3uA1wDXhL/bgHfB7CEBvB34PuAVwNvLB4UQo7ob3Bct5kSM/kZRTRERW/0jir1cOqKkUfTd/XeBr1WKbwLuDtt3A6+Nyt/tMz4BXGZmLwF+BLjP3b/m7l8H7mP5QSJ2jZMeVKxJM7zkC48FtS58c1X4om8/JunHh9WlVdZg8E7T99O93N2fDNtfAS4P21cAj0X1Hg9lufIlzOw2MztvZucPudizeWLXWIhQiYmFMOdd6RuffwJZNrPXrTskXiAmJfhNkU1CRAx+pLv7qHaTu9/h7ufc/dwBZ8Y6rVgXSaHuOKM0Pke5XWcNR9E7bXXbW8bzT61ggjNh5ttvS5vzL67j6w1hqOXPlsS9yTaict/HWJBebB19P+GvBrcN4fWpUP4EcFVU78pQlisXYjjJsMZqHbI5d5pCN8vEatPET4kpxWxfNKCbI3vtann2F8u6fVJiF+gr+vcCZQTOLcCHovI3hiie64GngxvoI8APm9kLwgDuD4cyIY7JWZ854vTKRRTTvuQjr16ndwvHJWpHXRstDkttSqu8dA2FdopFGuP0zey9wA8CLzazx5lF4fwC8AEzuxX4MvD6UP3DwI3ABeBbwE8CuPvXzOzfAb8f6v1bd68ODotdw4uZu8Am/cRnYt0zSc7j8j2Ka192rYw06bY3yeuHNs8fXn31us99g8h1pgfFLtMo+u7+hsyuGxJ1HXhT5jx3AXd1ap3YabzwTgOQNrF06GGLqJZcpEyO/Bq5kS/f8z+UWw8Ctxl3aOhfl3uo0E2hURux2ZSC1nKQciF1QU7oOz4A4HhB9OXmzRZIr1skfbmtietXysp+HPel4fzl/VEUj2hAoi9Ww0hRIVWr1rP5+UP9EQzbtkshdloyMUGbXyPV/o4Wqqmonb1Bn7TYCHq5Harujsw55hEyCcu6NnKH7qkVJiGkM0dS2JtcPNV+dYzikUtHxCjhmth8mgYm45z6DVk2h0buzHLuHD8IDn067IS5SKNqjH5DLn1Arh3RCln6YnUMdSF0OT4Tvgn9XT6pyVhdJmg1tiMW/E4nWuF9FVuPLH2xfsrQzsHnyYQ6htWoZhxbwzkfem7ANkc20qdannEpWdWqh8U+jDEpS2GYIqBHvFgtY1uV8aQlWHTtLIg9x+XQLXKnpTXfyf+feAAstTl+GFT7ORay8vcOWfpi9fQRmtKvX5d0bQVMrDiO0rEChvr0u1DtZ3wvuvrzJfZ7i0Rf7A5zK39xIHcpv83Cti2XRUzjxdETlapl05wbaG69W6Ks0lYq6ZWVc0eMiERf7BbV2avVNAw58U8Q++RTydZy+5pTK+e2U20XYlz0G0+slXRahYzA1rkwOgpkky7XD9wWjf77JuHvHEFU17/cfancR8XrC5Doi5NmFVEj1YVUyiUTK+6TqgulLYsWf3q7Wq8NSfeTR0slFjQvoDI2ivLZeST6YrtpEMJs9M7COSycqn4wtG+c/vy8ifMvRe80n6xFJSHySPTFxuCFN7sg2kSp7IIwtu1Dw/2QS0dU0UCu2D5SQre0EHrIUDl3l8T7OI6YaXrGhAopF09p408aInwW25XIBRS1MTlRK/RnuXFKuyC6I0tf7A7xAiQ5Bhq+8QBu12RsSWraM0+tLGtdjIhEX2wHVmPVVkWxqLwO1Mwy9j4l8nFZNka/DeWh1bbH1Il/3f0RIkKiLzaDNlEjbYStkmgtjpDpHrXTz5LvetxSVFEqDUPtCRruiyJyRIR8+mJzKTzvt64KnXsrcSv96W38+TFllM40zrHvEw7pkWnTK+1orF8si79Z/oEgd5CoQZa+2HxigW9h7S8sNVh4ehGVknn5SbtHrLYN5szHJBaWSWw8bbd7I4REX2wfKxS3BZ99QrGXJ2it0JUikRc9kHtH7BaxayP2jXOcWM3ifS2oLnpe5topKjZT18XRF9qRWvRFbhpxAsjSF7tDYsnEeEA3u1JVgly8fduQzdp4/dws4dwA7i5MNhMbg0RfnDx9oke6HtO0fmyGthE9dSI+wZsnZXW83hJd+7eKey62Eom+2Ho8FquEhbwgtB2jdmKa4vQ7E7XFqm6ecjt67xJlMQKDRN/M/rmZPWxmD5nZe83srJldbWYPmNkFM3u/mZ0Odc+E9xfC/peO0gOxl9ik5qtbZqik4jqpin9qO6Ju7ds+a+Y2XX9+WJxZM0Nt/4Woofc3x8yuAP4ZcM7dvxuYAjcDvwi8w92/E/g6cGs45Fbg66H8HaGeEEvUJQmzySQveB0HPtu4WqoCHy+Yklo8pU165c4unky/au8FSrYm0gw1F04BzzGzU8BzgSeBHwLuCfvvBl4btm8K7wn7bzBTzJkYmchCLhcaj63q7MzchsmvdROwplbU7veMS2mhLaGti66oljNyhehAb9F39yeA/wj8CTOxfxr4JPANd78Uqj0OXBG2rwAeC8deCvVfVD2vmd1mZufN7PwhF/s2T+wgS1ZtzmZITG6qiv1SxEwLypj8KcXcyk/F7tdSdeWw/BDKZtqcV1jst1w9ogtD3DsvYGa9Xw38VeB5wKuHNsjd73D3c+5+7oAzQ08n9oWadMRLA7kDWEijPCTBWkndQO5CPVn8YhyGmAh/G/iSu/+pux8CHwReCVwW3D0AVwJPhO0ngKsAwv7nA38+4PpC5KkKfiJ2f+1UhL5R+IUYgSGi/yfA9Wb23OCbvwH4PPBx4HWhzi3Ah8L2veE9Yf/H3PXtFiNSCnubyMYO1v+U+oXQy4XSqwO7Q645j+CRX1+MzBCf/gPMBmQ/BXwunOsO4GeAt5rZBWY++zvDIXcCLwrlbwVuH9BusSNYnEXTi/YThNokXitXzFoY3K0IaKPgx+6c4MePBm6rD4NGH3/Kfw/zAefWYf9tYyCi+2laaUswMPeOu78deHul+IvAKxJ1nwF+Ysj1xJ7i3j25WGIg93hf5bUD1ZDNQ6bdT1K5/nJ6iJ6WvX4RiBZo2F9sHx0fALkQzdx7K1fK6un8L4+z+MI118u2sQlFPIseSPTF9jKZJIUvdplYPHgb/tpOyoqZ4nMrP861Mwvf9ESK5ZaTtKIY/jhUM5tP32zWbyF6om+P2Gy6ziot67cdGnDA6y3mtqGZjfXc2ln083VyvX//hcgg0Re7Q26QNjERqs+aubA4M7fzMolRG/Lr4iYOkK9ejIhEX2wuqUieLoujw3HitSKKjMnOf1o+d9csmqn6qfPG7ZgfEh/aZ7nE+bHKxinyaOUssZ108GuPMXF2Fosf/PhBpQufhiRrw0W2cxsnk95rCIj9Rpa+WBtjxY2bLX+NlxYYrwyaHr+3Vn79kjiip1N0T3ydVDtSba72KdHPrihWX8jSF7tFC9dGV6u6re++q4+/sR1y04gTQJa+2F5y/v1qBMvIES0Tim6LoLehS5sVny8GINEXG8ngBUDmeWtmb8uImYWomcQl4kHX5Vj9ZaFPlcXHJQdxPRNFVKaMGCFaRwuoiBwSfbHd9PFRd0l8Fpjn0k+EbKYmZ4157QXkkxcDkeiL7aPrjNTYgm6TkiG+VMLxPjVnWrN+bm07sm3LH5ZvnP77iu7oWyM2B5vM/sYis8B46dpZSI1TJ/wtFbmuXnz+1PXnlcYeux3zfoqdQN8IsXtU/OJ1xvdi9s1618nUFgdwJ9SvjZs6b64tWhtXrAqFbIq9IJv2oCVNi6h0ZmA6CCH6IktfrI1BESY1i6IvVKvz4XdIxxATL4xeR1P6hcb0ynUWf8+wTUX1CIm+2H4yES3VVakWZsL2uQyLA7hT89459+fM2zU7TzalMihyR4yCRF9sL3UiWCeeDhZtL7xWmEXqtJ+Rm4vqWVotq+aajf58ib8YgERfbDWN+WhqBLTJlx776mN3TrkQeryvur+O7HVbDN6OkX9H7Df6BondIJliuHSZLM/IXTh0RW7u3LUXZ+WO68MXoopEX+weDRZznIKhbaw+BPdNZWH0Nq6fpRh9aF62USGb4oSQ6IudxgoqqQ9ile93zrbRO1lS7fFoMRUhThDF6YvdpCqgc7G3TDqEmfukCK9xOoV4AZUU5b7Yl1+WledLTvzKuXQk/uIEkaUvdpc2g7hOhwVUioqwd0yx7JZ0Ky3WkVtHnCyy9MV20FcM56tRzd4upjtoPjwl6mVY5lHi+FYPgYqPf/Y3MPWCHhaiJYMsfTO7zMzuMbM/MLNHzOz7zeyFZnafmT0aXl8Q6pqZvdPMLpjZg2Z23ThdELvIwrJ+TbNIO0a2WN+slmNSCn0XmvoZ3SctiyhyDHXv/DLwv9z9u4C/BTwC3A7c7+7XAPeH9wCvAa4Jf7cB7xp4bbFn+BgLgSeyXc4jaWo0uBy4rUbvlK+1A7ueuFaiPX0Z5b6IvaG36JvZ84FXAXcCuPuz7v4N4Cbg7lDtbuC1Yfsm4N0+4xPAZWb2kr7XF6IVTkixzIBondmBqdm2C2kZ+l4gbqMQJ8wQS/9q4E+B/2ZmnzazXzOz5wGXu/uToc5XgMvD9hXAY9Hxj4eyBczsNjM7b2bnD7k4oHliL4jdGFX3R84AHpyDZ4Q1cnNtSJ027pfcNmIgQ0T/FHAd8C53fznw/zh25QDg7p3/e7n7He5+zt3PHXBmQPPEXlFZRSrnL1+IhW+YIJXNo9NA3XFVN1IqNj/Zdq2SJUZiyDfpceBxd38gvL+H2UPgq6XbJrw+FfY/AVwVHX9lKBPiRCmjd5KLkUMrs2S+Hm4ly2a5r91iKnGbqu2Rb0esht6i7+5fAR4zs78Rim4APg/cC9wSym4BPhS27wXeGKJ4rgeejtxAQpwIrcW0nBhbqT5Z8NmXwl9Ei6IX2foL52vRDAm/WAVD4/T/KfAeMzsNfBH4SWYPkg+Y2a3Al4HXh7ofBm4ELgDfCnWF6MZQYayx8HOLnrRdGStXb+m8Y0bv6EEhOjJI9N39M8C5xK4bEnUdeNOQ6wnRSC6WPY71r/rVO+hmdY3cTk2rin35PjcPQZk1xQmgGbliZzH3ma72FPiYCT635EuXzhGTkIqhv7VdXbNXLh5x0igkQOwdCzns56J7wlb1PPFafW5/IU4aib7YfurcIFHq4vzyhO0vFUfptF1GsfE6uWybVeTuESMg0Rf7RV0+/YzmVmfaTkgviJ6dkVt3HblzxIqRT1/sPlVhjaz+qk/d3Y5z4FMJ00xl3KysnhXXKdyOI3dS15P4izUg0Re7Q437o35pwvaXSCVca00uSKfp+nLriBGRe0fsBi2FsToTtg+tZ+DWtKHaluaDJPxiHCT6Yi+punUgPzmrpCnLZo6F8w584AgxFIm+2D665I8vytBMT/v2ac6lHy+LGA/glttl1s3ambueftDM3oeVvbr+eFAefdED+fTFfuEsDqL2sLjLKJ2jvtdPtUOIFSFLX4gMfWfaDpmhK8RJI9EXu0tTGGQ8GzdY3XV+/Uk1PNOaF1Px6NxUwjd7tVmIgUj0xV4Qpz6IffhNA6rVVMlQv3JWqn61HYDSMYi1IZ++2C0sYceMYD3XRen0XWFrgVQb533RgK0YD1n6YvdoEdNuDvNaDXn141m5VXKzcVPnmf+6QBOyxPqQpS/2i3noZsW/7kb0GFhIxRAzxedCX3jeZlo83tL+/FQYqRAnjERf7AcpbU2UjS3ByfO1TPQmxEkg0Rf7ScdUCHV5dtrm4Fm4loRerAmJvtg7UiJv3qzD0xC1czw5K515s8s1hVg1GsgVO011+cFYaJdENyPCTWGYjfUqxbk2aKlEsQok+mLnSaXEsTYratF9daza+gvzAyoPI0VlihUh947YP2oGdXMzcssF0KdWzC36KQUT99pEa9VFVBrbIcQJI9EXe0HSdZIR3Vy4Zley51ly90j9xeqQe0fsLqkJTtXomVVE1KRSOOeuqUlZ4oSRpS92Em8jnh1FfrZOrs8jdgqs3wpaTXngzJD0i5NisKVvZlMz+7SZ/c/w/moze8DMLpjZ+83sdCg/E95fCPtfOvTaQiwwablkYqTTTQuoQFj8vCFOv6kOsPirgg6Dty37JUQbxnDvvAV4JHr/i8A73P07ga8Dt4byW4Gvh/J3hHpCrIVkdGWZBrmBmbXf8mdCnFK56fpCrIBBom9mVwI/CvxaeG/ADwH3hCp3A68N2zeF94T9N4T6QqyF2pj9DHEsftv4/b7XEuIkGGrp/2fgpznO/foi4Bvufim8fxy4ImxfATwGEPY/HeovYGa3mdl5Mzt/yMWBzRNikSV3TsK9kwumSVn3dRb/0nkS19YDQKya3qJvZj8GPOXunxyxPbj7He5+zt3PHXBmzFMLcUznQdxFCz+28qfm3XPqS+zFmhgSvfNK4MfN7EbgLPBXgF8GLjOzU8GavxJ4ItR/ArgKeNzMTgHPB/58wPWFaE9q0HTB6raQ6XjmccynVs7nz1+4XDjew3kX/PpVwddsXLFCelv67v42d7/S3V8K3Ax8zN3/HvBx4HWh2i3Ah8L2veE9Yf/H3DUrRayZFi6WMiwztURiXNYUvtkmUkiIk+YkJmf9DPBWM7vAzGd/Zyi/E3hRKH8rcPsJXFuIUZlk/Pg5X36qvhCbxCiTs9z9d4DfCdtfBF6RqPMM8BNjXE+IISQXJG9phcdi3y1ss+H6QqwIzcgVe01VeN1nM22rpFw7TfsKbCmCR0Iv1o1y74idZikdQ9FfdUuffeze6ZWGIdOeVqkjhBiIRF/sL179axbdXpOzypm+8Z8Qa0KiL/aT3CJXI6VVbjyfhF+sCfn0xe5QeG8zps5oLxdPKV068/xnHlw985DO/EkG+fIHuKSEqCJLX+wfmXVzy1d3Sw7mQvvondkgriXPn2uHEKtAoi/2lnWETSpUU6wbib4QMS18+t1SKysiR2wWEn2xlyQnZjXo+MR8yb3TKoIncX5Z+2JdaCBX7Dctxbdp5ayxryfESSHRF/tL1doP5DJslrT9eVzkMmtK+MUakXtHiEiE44CaMoInnnUbL1dbXbq2rBdH/rjEXmwYEn0hSiLLPGXtp9w4OdfOopWvwVyxOUj0xX6Ti6GvYWqzvy7Mzy9rX6wZ+fTF3lG7mImn50yVs21Tk7NyM3G94TqK4BHrQJa+2F+qVn7TClpz4Tempb+/jene8TpCnCQSfSGE2CMk+mKvmbvmT3qwtczDc7JXEaIR+fSFKPFZsjV3y8bqTzi2lOospiKcx8tc+kJsCBJ9IdqsjWtFMuXCxLz96lkSf7EByL0j9ouaaJqdvK4QFWTpC1GSCbGMFz6fxehHkTuertd0TiHWhSx9sbfk4vXLxU+Oeg7ulscll0pUfL5YMxJ9sfe0FeFUTH7rvPodriPESSL3jhAsC3IuemeKMZlPzMr/EqgeL8EXm0JvS9/MrjKzj5vZ583sYTN7Syh/oZndZ2aPhtcXhHIzs3ea2QUze9DMrhurE0J0JSnCbmmXDPUWfd1KWrOQzeVz6iEg1sUQ984l4F+4+7XA9cCbzOxa4Hbgfne/Brg/vAd4DXBN+LsNeNeAawtxYjTl01/1eYQYk96i7+5PuvunwvY3gUeAK4CbgLtDtbuB14btm4B3+4xPAJeZ2Uv6Xl+IdRHn3hFi2xhlINfMXgq8HHgAuNzdnwy7vgJcHravAB6LDns8lFXPdZuZnTez84dcHKN5QmSZR/C0maCFL/yHmdBhIDdcQ24dsW4Gi76ZfRvwG8BPuftfxPvca5PLJnH3O9z9nLufO+DM0OYJ0Zo4hHOoa2Z+vIRebBiDRN/MDpgJ/nvc/YOh+Kul2ya8PhXKnwCuig6/MpQJsVH4XPj7/fcoj0vl5Rdi3QyJ3jHgTuARd/+laNe9wC1h+xbgQ1H5G0MUz/XA05EbSIitYGqzkM0JNp+ZK8Q2MSRO/5XA3wc+Z2afCWX/CvgF4ANmdivwZeD1Yd+HgRuBC8C3gJ8ccG0hVsqEYmkhdJgtjp5MvyDEhtJb9N39f5NPD35Dor4Db+p7PSFGJ3K/HPvzx7bebdmvL7ePWCOakStEgtwkLThOuFYn3nXHC7FOlHtHiIp41wn2JPolMGn4VbB0Hln4YgOQ6AshxB4h944QMZE1Xo3VLydiTea20lFyctbCcbLuxYYhS1+IGoqOA7td6wuxaiT6Yr/JWOI58Y5z7tSmVs7tk+Uv1oxEX4gKY0XeKIJHbCLy6Yv9I5lLv/3hTVE7J3FNIcZCoi9EC6aWn3Vbt0+ITUPuHSFiaqzvCctx+o3/gWTNiw1Doi+EEHuE3DtCZMhG8FiwlTpG/gixCcjSF3tPdZETZ5xFVKrPBC2mIjYBWfpiL7AWK5qYn4wLvkns27RNiLGQpS+EEHuERF+IDsSrZWnlLLGNSPSFEGKPkE9fiBqOEoujT2pspVR9ITYJfUOFaMk04c1JlQmxyUj0hYgZO0makq6JDUOiL4QQe4REXwg4+Rw5CsUXG4JEX4gEQ3PhK5e+2FQk+mI/aTELtm8qhlbHaRauWBMSfSGE2CNWLvpm9moz+4KZXTCz21d9fSGa0HKJYpdZqeib2RT4FeA1wLXAG8zs2lW2QYghxBOz6iZpCbGprPpb+wrggrt/0d2fBd4H3LTiNgghxN6y6jQMVwCPRe8fB74vrmBmtwG3hbcXP+r3PLSitq2KFwN/tu5GjEi7/uTGLTdvedkXA3/2KPB/Oh/6pdEbMxL7+Z3bHk6iP38tt2Pjcu+4+x3AHQBmdt7dz625SaOya31SfzafXeuT+jOMVbt3ngCuit5fGcqEEEKsgFWL/u8D15jZ1WZ2GrgZuHfFbRBCiL1lpe4dd79kZm8GPgJMgbvc/eGaQ+5YTctWyq71Sf3ZfHatT+rPAMw1M1AIIfYGBRoLIcQeIdEXQog9YmNFfxfSNZjZH5vZ58zsM2Z2PpS90MzuM7NHw+sL1t3OOszsLjN7ysweisqSfbAZ7wyf2YNmdt36Wp4m05+fM7Mnwuf0GTO7Mdr3ttCfL5jZj6yn1XnM7Coz+7iZfd7MHjazt4TyrfyMavqzzZ/RWTP7PTP7bOjTvwnlV5vZA6Ht7w/BLZjZmfD+Qtj/0lEb5O4b98dskPePgJcBp4HPAteuu109+vHHwIsrZf8euD1s3w784rrb2dCHVwHXAQ819QG4EfgtwIDrgQfW3f6W/fk54F8m6l4bvntngKvDd3K67j5U2vgS4Lqw/e3AH4Z2b+VnVNOfbf6MDPi2sH0APBDu/QeAm0P5rwL/OGz/E+BXw/bNwPvHbM+mWvq7nK7hJuDusH038Nr1NaUZd/9d4GuV4lwfbgLe7TM+AVxmZi9ZSUNbkulPjpuA97n7RXf/EnCB2XdzY3D3J939U2H7m8AjzGa+b+VnVNOfHNvwGbm7/9/w9iD8OfBDwD2hvPoZlZ/dPcANZjZa9r5NFf1Uuoa6D35TceC3zeyTIb0EwOXu/mTY/gpw+XqaNohcH7b5c3tzcHfcFbnctqo/wQ3wcmaW5NZ/RpX+wBZ/RmY2NbPPAE8B9zH7RfINd78UqsTtnvcp7H8aeNFYbdlU0d8VfsDdr2OWVfRNZvaqeKfPfr9tdczsLvQBeBfw14HvBZ4E/tNaW9MDM/s24DeAn3L3v4j3beNnlOjPVn9G7n7k7t/LLAvBK4DvWldbNlX0dyJdg7s/EV6fAn6T2Yf91fLndHh9an0t7E2uD1v5ubn7V8N/ygL4rxy7B7aiP2Z2wEwg3+PuHwzFW/sZpfqz7Z9Ribt/A/g48P3MXGvlBNm43fM+hf3PB/58rDZsquhvfboGM3uemX17uQ38MPAQs37cEqrdAnxoPS0cRK4P9wJvDBEi1wNPRy6GjaXi0/67zD4nmPXn5hBNcTVwDfB7q25fHcHXeyfwiLv/UrRrKz+jXH+2/DP6DjO7LGw/B/g7zMYqPg68LlSrfkblZ/c64GPh19o4rHtku2bE+0ZmI/d/BPzsutvTo/0vYxZV8Fng4bIPzHxz9wOPAh8FXrjutjb0473Mfk4fMvM73prrA7MohV8Jn9nngHPrbn/L/vz30N4Hw3+4l0T1fzb05wvAa9bd/kR/foCZ6+ZB4DPh78Zt/Yxq+rPNn9H3AJ8ObX8I+Neh/GXMHlAXgP8BnAnlZ8P7C2H/y8Zsj9IwCCHEHrGp7h0hhBAngERfCCH2CIm+EELsERJ9IYTYIyT6QgixR0j0hRBij5DoCyHEHvH/AbTnaN7deozWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolor(Umag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "d5acec66",
   "metadata": {},
   "outputs": [],
   "source": [
    "from nodepy import rk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "44097c30",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classical RK4\n",
      "The original four-stage, fourth-order method of Kutta\n",
      " 0   |\n",
      " 1/2 | 1/2\n",
      " 1/2 |      1/2\n",
      " 1   |           1\n",
      "_____|____________________\n",
      "     | 1/6  1/3  1/3  1/6\n"
     ]
    }
   ],
   "source": [
    "print(rk.loadRKM('RK44'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5524d9d5",
   "metadata": {},
   "source": [
    "# Double-relaxation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acd9cf3b",
   "metadata": {},
   "source": [
    "Here I have tried to conserve both invariants, by using an embedded method and \"relaxing\" along the increment directions for both methods.  However, fsolve is unable to find a solution most of the time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "73456b44",
   "metadata": {},
   "outputs": [],
   "source": [
    "def E_1(U, dx):\n",
    "    V = U[::2]\n",
    "    W = U[1::2]\n",
    "    return dx*np.sum(V**2+W**2)\n",
    "\n",
    "def E_2(U,dx):\n",
    "    V = U[::2]\n",
    "    W = U[1::2]\n",
    "    return np.sum((np.diff(V)**2+np.diff(W)**2)/dx - 0.5*q*dx*(V[1:]**2+W[1:]**2)**2)\n",
    "    \n",
    "    \n",
    "def fgam(gammas,U,inc1,inc2,E10,E20):\n",
    "    gamma1, gamma2 = gammas\n",
    "    uprop = U + gamma1*inc1 + gamma2*inc2\n",
    "    E1 = E_1(uprop,dx)\n",
    "    E2 = E_2(uprop,dx)\n",
    "    return np.array([E1-E10,E2-E20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "8b45ee29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4872 5504\n"
     ]
    }
   ],
   "source": [
    "# Initial condition\n",
    "U = np.zeros((2*m))\n",
    "for i, xc in enumerate(x):\n",
    "    U[2*i] = 1./np.cosh(xc)\n",
    "    \n",
    "# Integrate with RK4\n",
    "\n",
    "dt = 0.05*dx*dx\n",
    "T = 1.\n",
    "t = 0.\n",
    "UU = []\n",
    "UU.append(U)\n",
    "q = 18.\n",
    "\n",
    "E10 = E_1(U,dx)\n",
    "E20 = E_2(U,dx)\n",
    "\n",
    "ret = np.zeros_like(U)\n",
    "F = np.zeros_like(U)\n",
    "errs = 0\n",
    "steps = 0\n",
    "\n",
    "while t<T:\n",
    "    f1 = udot(UU[-1],dx,q,ret,F)\n",
    "    y2 = UU[-1] + 0.5*dt * f1\n",
    "    f2 = udot(y2,dx,q,ret,F)\n",
    "    y3 = UU[-1] + 0.5*dt * f2\n",
    "    f3 = udot(y3,dx,q,ret,F)\n",
    "    y4 = UU[-1] + 0.5*dt * f3\n",
    "    f4 = udot(y4,dx,q,ret,F)\n",
    "    inc1 = dt/6 * (f1+2*f2+2*f3+f4)\n",
    "    inc2 = dt/4*(f1+f2+f3+f4) # Only 2nd-order accurate\n",
    "    \n",
    "    # Solve for gamma1, gamma2 that give conservation\n",
    "    gammas, info, ier, mesg = fsolve(fgam,[1.,0.],args=(UU[-1],inc1,inc2,E10,E20),full_output=True)\n",
    "    if ier != 1: errs += 1\n",
    "    #print(t,mesg)\n",
    "    #print(fgam([1.,0.],UU[-1],inc1,inc2,E10,E20))\n",
    "    gamma1, gamma2 = gammas\n",
    "    #print(sum(gammas))\n",
    "    Unew = UU[-1] + gamma1*inc1 + gamma2*inc2\n",
    "    \n",
    "    UU.append(Unew)\n",
    "    t+=(gamma1+gamma2)*dt\n",
    "    steps += 1\n",
    "    \n",
    "print(errs, steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "522da83e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc084a99a00>]"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEDCAYAAAD9ZJllAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo/klEQVR4nO3df7xVdZ3v8ddbfokiyBEC5JCCQqWUFEegKU0SAcsJnabCMTlT/pg0m5nq1uA4RaM5l/HWNHnvjF1kKJybmne05DoaHSiqaQblkKiQGmAmHBHQg1Jp+IPP/WN/N65z2Pv82uv8fj8fj/XYa332d631/e6zz/7s9V1rf5ciAjMzszwc0d0VMDOzvsNJxczMcuOkYmZmuXFSMTOz3DipmJlZbpxUzMwsN04qrZD0IUlbJB2UVNNCuRWS9kja3Cz+PyQ9JulhSd+VdGyKD5b0TUmPSHpI0lmZdaan+DZJN0pSildJqpO0NT2OTHGlctvSft6R2VZtKr9VUm1n7KO9r4mZ9V1OKq3bDPwR8JNWyn0LmF8iXgdMjYi3Ab8Erk7xywAi4q3AOcBXJRX/Hjel5yenqbjdxcDaiJgMrE3LAOdmyl6e1kdSFbAEmAnMAJYUk0Re++jga2JmfZSTSisi4tGIeLwN5X4CNJaI/yAiXk2L64HqNH8K8MNUZg/wPFAjaRwwPCLWR+GXqbcA56d1FgAr0/zKZvFbomA9cGzazjygLiIaI2IfhQQ3P+d9IOlzkjakI5i/be01MbO+y0mla30cuC/NPwR8QNJASROB6cAEYDywM7POzhQDGBMRu9L8M8CYND8e2FFinZbiuexD0lwKRy8zgGnAdElnln8JzKwvG9jdFegJJK0BxpZ46pqIuDunfVwDvAp8O4VWAG8B6oFfA/8JvNbW7UVESOrUMXbauI+5aXowLQ+jkGRa6y40sz7ISQWIiDmduX1JfwqcB5yduptIXWKfzpT5TwrnXPbxehcZab4hze+WNC4idqWupz0p3kDhKKf5Og3AWc3i61I8r30I+O8R8b9bfSHMrM9z91cnkzQf+DzwgYh4MRM/StLRaf4c4NWI+EXqetovaVa6ImsRUDxaWgUUr+CqbRZflK7QmgW8kLazGpgraWQ6QT8XWN0J+/i4pGGpLeMlvSGHl87MeqOI8NTCBFxA4fzBAWA3hQ9lgOOBezPlbgN2Aa+k8pek+DYK5yI2pekbKX4i8DjwKLAGOCGzrRoKV51tB/4XoBQ/jsIVWVvTOlUpLuCfUvlHgJrMtj6e6rAN+Fgn7eMvUuwR4L+Ak1p6TTx58tR3p+IHiZmZWcXc/WVmZrnp9yfqR40aFSeeeGJ3V8PMrFfZuHHjsxExunm83yeVE088kfr6+u6uhplZryLp16Xi7v4yM7PcOKmYmVlunFTMzCw3TipmZpYbJxUzM8uNk4qZmeXGScXMzHLjpGJm1kvs3PciP3p8T+sFu5GTiplZLzHvaz/hY9/c0N3VaJGTSif4/Suv8fOn9nV3Ncysj/ndy22+j1+36ffDtHTUP9T9kod3Ps+6x/e2WG7kUYPY9+IrnPe2cbzw0iv85/bn+HDNhMPKSaXXLxNus3LbbXW9Duw5CF55NUru9/tbnmHs8CN5W/WIQtmA4vjY2YGypabLpdz585286+TjGDdiaJvq1fi7l/nhY3v40PTqJnXLtrGtr9PtG16/o/LC0yekdhQqXKx3NJl/PViuWS2NFP69TU9zWvUIJo46mu9tepoZE6s4fsSRAPzgF7sZfuQgZk2qYv0TjTyz//cAfOC04zmiWXtefu0g9z7yDO9/6zgGDRDP7P89659o5Pxpxx8qs/bRPQwdPIB3nTyqbS8GsPc3B/iPbc9ywdsLd6Nu/N3L/PiXr/9PnDbhWCaNOpqfbXuWA68e5L1vfv1WOw/8qpEXXnqFc04Zc9h226qj/x93P/Q0J48exqnjh3dsAwGvRXD3pqeZe8oYhg0ZyF0PNjBzYhWDBx7BT7c+26T4H5x0HM+88Ht27HuRP3zb8WU22nZnf3UdTzz7O0YMHcTzL75yKP6RmgmtvpdffPk1Vj30NB98RzVLP/hWBg3I99ii4qHvJVUB36Fwf5AngQ9HxGFf0yXVAn+TFr8cEStTfDrwLWAocC/wFxERkr4EXAYU36F/HRH3pnWuBi6hcPvdP4+I1Sk+H/g6MABYHhFLW6t/TU1NdGTsr4XL/ov1TzS2ez2AUcMGo2Z/+dJ/hsr+Nh3903Z0rxHB4IGl36C79x8AYOzwwgeiVPhAOPx1iMNiWa8djEMfnuOPbVtSaXj+JQCGDhrAyKMGHZbMoo0tfu0gPPvbA01iY4YPQSiTqNKjDk9YhTaXblupJh+MYEdjoe7HjziSp18otPuE444C4NfPvXhouThfVCxTlH3+jVVH8VTji4fqf+SgAU3KvLGq6botKW5nxNBBjBg66NBy1oSqoYfakd12sWx79pfV1r9bKcX6VI9s23uolJ37Xjo0f/TgAW0+iphQ1fF9FuvdkjHDh7T4fPF/EeCx6+Yf+vu3l6SNEVFzWDyHpHID0BgRSyUtBkZGxF81K1NF4V7sNRQ+szYC0yNin6QHgD8H7qeQVG6MiPtSUvltRHyl2bZOoXDzpxkUbpS1BpiSnv4lcA6FG0JtAC6MiF+0VP+OJhWA069fw97fHGi9IHDq8cPZ8vR+ALb87TyOHtK/DhJPXPzvADy59P0VbWf3/t8z8+/WMmb4EO7/67bdBfqCf/4ZDz71PHde8QdMP2Fkh/f93G8PMP3Law4tf/vSme36Vt9eL778Kqd8cTUAP/38bM644UeMP3YoP1v8XqDpa1qcB/j6wmksmDa+ybb+cc0v+cc1W/nz957MZ+a+iTNu+CE7Gl9i3X87ixNHHX3Y9trq3K//lEd37eeeT72bqeNH8NHl9/Mf25p+S8/WL7vtvN4THZHHvm+9/yn++ruPsPD0CXxm7hRmXL+WUcMGUz3yKDbteL7p/o47iidT0q5kn7teeIlXXwvGjjiSV147eOj9UbTy4zN4z5TDBg5u4uq7HuG2B57i+gumctHMEzpcl3JJJY/jngXAyjS/Eji/RJl5QF1ENKajmDpgfroH+vCIWB+F7HZLmfWb7+/2iDgQEb+icEfDGWnaFhFPRMTLwO2prFmS7w3pOvv+dh3pgrTukf1bdebbYtyIoUyoOopBA47gqMGHfzHtCTddzCOpjInCvcoBngFKdZCOp3BL3aKdKTY+zTePF10l6WFJK9I91lvbVqn4YSRdLqleUv3evS2fEzErOqyrrlM/Ppp2iXX03Fi5jRY/BCvdbvPVO/s16Uly+5v0MW1KKpLWSNpcYmpyJJCONvJ6V90EnARMo3Cf86/mtF0iYllE1EREzejRLR8qWl9S2afAYR+gXfj5mdu+0oaaX1zQ4c1VWp8+qL/nmjYllYiYExFTS0x3A7tTNxbpsdQvcxqA7CVP1SnWkOabx4mI3RHxWkQcBG6m0L3V2rZKxTvN1z8yrc1lO+Vbp3Wprv67+X3SO/X3RJtH99cqoDbN1wJ3lyizGpgraWTqxpoLrE7dZvslzVKhb2FRcf1iokouADZn9rdQ0hBJE4HJwAMUTsxPljRR0mBgYSrbaf6ghZO0n5v3ps7ctXWD5uc4OvvDo1POqXRy91d/MmpY4Sqrtl6F2F/kkVSWAudI2grMSctIqpG0HCAiGoHrKHzwbwCuTTGAK4HlFE64bwfuS/EbJD0i6WFgNvDptK0twB3AL4DvA59MRzSvAldRSGCPAneksp3qnk+9+7DYuVPH8snZJ3f2rq3dKkwDzT5BO/ukaKccqbj7Kzdz3vIGll08nSvOOqm7q9KjVHxda0Q8B5xdIl4PXJpZXgGsKFNuaon4xS3s83rg+hLxeylcltxlpo4f0ZW7s27U/EO+849Uyu/bup8k5p46trur0eN4mBbrVXrAFZOv6+xLijOZJLd2u/vLOln/+gVeD+HfH1hb9MZ3SX99bx939GDOnTqWS949kTt/3sBDzX78CHDbZbPYse/wEQf6GieVTnCE+yp6qMr+Ls2PFnrl71Ry1pMOHLvTEUeImz46HYC3Vo/gtgeeOqzMO086jndyXFdXrcu5+ysHP/j0mU2Wr3n/W7qpJtaVOv0X9SmTDDxCuf9OxTrPkIEdG0urr/CRSg6mjDmmyfLxvsTQcrL43De3OpZTd2p+AHXq+OFNxv5aMK38iLyXvHsiu15ofYBE612cVLpIR4ZZt8ONPmYI504dy6VnTOrA2pUOntpsa13wpf8T7ylcrrozr774Ft58f/P+tzCsnQOdNn8JPjf3TcyaeBwf+9YGBg84gq9+6LSy637hvFPata+8+f+wczip5OzCGW/s7ir0aQMyfdfd7WBv7Epqoc4dS9RNDRxwxKF75gw7ciADc75XR15WXfUuRh/T8hDx1jFOKjmb2tGb/lgXqOyr6eAe+gHZnXrrl/23VR/b3VXos/xfkpPZb2q53/vL5x/2+07rZY4cNIBF73z9/hO98DjFfT7W6ZxUcjK2hVvbXrvgVE6bcGzXVcY6zZh058peqzd22fUBXdXV1hP+uk4qXWDRO0/s7ioY0Ntu0mV9R085D9gVnFTM2qHpIJK9MKvk3P01YuggAAYOeH27Q9I9z9/uo/NDiiMa9wc+UW/9SL4fqL35SOWowYUP/kpHf7jxwrfz/x56mjdlfqs1bMhA7vnUu5k0+uhDsU1fPIeDvfj1srZzUukGPlfaN3TlZ+SxRw0G4EM1r9/T7p2TjmPcsYVzPJ+b9yZu3/AUOxpfKjly9h9Pr+bW+5/iQ9ML6//Ln57O3ZsaqB5Z2Q91Rx8zhI+/e+Jh8eZ1KNa/v1n72ff0u6sGnVTMOqgrj1SGDRnI1uvPZeARr38jue3yWYfmPzn7ZD45+2QiosnoxkXVI4/igWvmHFoef+xQrjzL9/zpbCeNHtbdVehyTiqd6IeffQ9DB/fvcYAsP4Pa8I23VEKxvkvqed2wFR2XSaqSVCdpa3ocWaZcbSqzVVJtJj493d1xm6Qb0y2FkfQlSQ2SNqXpfSl+jqSNaZ2Nkt6b2dY6SY9n1nlDJW3Lw6TRwxjXwqXG1vtk/4E7e5Ris9b0xK8QlXb2LQbWRsRkYG1abkJSFbAEmAnMAJZkks9NwGUU7jM/GZifWfVrETEtTcW7OT4L/GFEvBWoBf612e4uyqyzp8K2dZr+es8JM+v7Kk0qC4CVaX4lcH6JMvOAuohojIh9QB0wX9I4YHhErI/CdZq3lFn/kIh4MCKeTotbgKGS+s+1etbtmlxQ7AMV62Yd6+7s3DdupUllTETsSvPPAGNKlBkP7Mgs70yx8Wm+ebzoKkkPS1pRplvtg8DPI+JAJvbN1PX1BbXwaku6XFK9pPq9e/eWb51ZC3rlgJLWp53cAy4MaDWpSFojaXOJaUG2XDrayOu/7CbgJGAasAv4arM6nQr8PfBnmfBFqVvsjDRdXG7jEbEsImoiomb06J57rwozs/aYUHVUG0p1bvd7q1d/RcSccs9J2i1pXETsSt1Zpc5jNABnZZargXUpXt0s3pD2uTuzj5uBezLL1cB3gUURsT1Tz+K6v5F0K4XzN7e01r7u4At0zPq27//lGRw9uPMvru2JHyWVdn+tonDCnPR4d4kyq4G5kkambqy5wOrUbbZf0qzUVbWouH5KUEUXAJtT/Fjg34HFEfGzYgFJAyWNSvODgPOK65i9zmN/Wdd489jhbTxqqExP/IJaaVJZCpwjaSswJy0jqUbScoCIaASuAzak6doUA7gSWA5sA7YD96X4Demy4YeB2cCnU/wq4GTgi80uHR4CrE7lN1E44rm5wraZHcaJxKxlFR2fRcRzwNkl4vXApZnlFcCKMuUOu9FIRJQ8HxIRXwa+XKY6/WcYUOugnMf+8u9UrJsVfp7Qs96H/WtQmh6iBx6x9mnnTh0LFIYmqVQ2kfioxbpdD/ww8TAt1udddsYk/mTmCQwbku/b3UnF7HA+UrE+T1LuCQV6WqeD9UfDj+x5xwVOKmbtkD06ecu4Y8oXNOsC/+fSmd1dhcM4qXQDjyTbN5x6/OH3LTHrSm/sgsuW28tJxczMcuOkkpOh6b7cbbnnhZlZHnriiOc97yxPL/XZuVM45siB/NHbx7de2Hotn5y3nqRjPemd+y52UsnJ0UMG8ulzprSpbM/7bmFmlg/31ZiZ9Sud+7XWScWsPfyLR7MWOamYmfVSPfHXCU4q3aAnvhHMzPLgpGLWDr9/9WB3V8HskJ54SbGTilk7LP/pE91dBbMezUnFrB0O+jy99SA9sSu9oqQiqUpSnaSt6XFkmXK1qcxWSbWZ+PR0h8dtkm5MtxVG0pckNWTu7vi+FD9R0kuZ+Dda25aZWV/VEz/kKj1SWQysjYjJwNq03ISkKmAJMBOYASzJJJ+bgMuAyWman1n1axExLU33ZuLbM/FPZOItbatHcb4zs76q0qSyAFiZ5lcC55coMw+oi4jGiNgH1AHzJY0DhkfE+ogI4JYy67cqz22ZmfUWPfELaqVJZUxE7ErzzwBjSpQZD+zILO9MsfFpvnm86CpJD0ta0axbbaKkByX9WNIZmX20tK0mJF0uqV5S/d69e1tqn5mZtUOrSUXSGkmbS0wLsuXSEUJepzFvAk4CpgG7gK+m+C7gjRHxduAzwK2Shrd34xGxLCJqIqJm9OjROVXZzKxr9bzjlDYMKBkRc8o9J2m3pHERsSt1Qe0pUawBOCuzXA2sS/HqZvGGtM/dmX3cDNyT4geAA2l+o6TtwJSWtmVmZlmdewljpd1fq4Di1Vy1wN0lyqwG5koambqx5gKrU7fZfkmz0pVai4rrpwRVdAGwOcVHSxqQ5idROCH/REvbMjPrq3rgKZWKh75fCtwh6RLg18CHASTVAJ+IiEsjolHSdcCGtM61EdGY5q8EvgUMBe5LE8ANkqZRSKlPAn+W4mcC10p6BTiY9tHatszM+qSOnajv3ExUUVKJiOeAs0vE64FLM8srgBVlyk0tEb+4zP7uBO4s81zJbfUkf/2+N/N39z7W3dUwM+s0vklXF7r8zJO4/MyTursaZmadxsO0mJlZbpxUzMwsN04qZmaWGycVMzPLjZOKmZnlxknFzMxy46RiZma5cVIxa4czp3gAUrOWOKmYtcPxI47s7iqY9WhOKmZm/UrPHqXYzMzsECcVM7N+pXNHKXZSMTOz3DipmLVDdG53tFmv56RiZma5qSipSKqSVCdpa3ocWaZcbSqzVVJtJj5d0iOStkm6Md0KGElfktQgaVOa3pfiF2VimyQdTHeIRNI6SY9nnntDJW0zM7P2q/RIZTGwNiImA2vTchOSqoAlwExgBrAkk3xuAi6jcK/5ycD8zKpfi4hpaboXICK+XYwBFwO/iohNmXUuyqyzp8K2mZlZO1WaVBYAK9P8SuD8EmXmAXUR0RgR+4A6YL6kccDwiFgfEQHcUmb9ci4Ebu9oxc3MLH+VJpUxEbErzT8DjClRZjywI7O8M8XGp/nm8aKrJD0saUWZbrWPALc1i30zdX19odiVVoqkyyXVS6rfu3dvuWJmZtZOrSYVSWskbS4xLciWS0cbeV0bcxNwEjAN2AV8tVmdZgIvRsTmTPiiiHgrcEaaLi638YhYFhE1EVEzerTHcjIzy8vA1gpExJxyz0naLWlcROxK3VmlzmM0AGdllquBdSle3SzekPa5O7OPm4F7mm1zIc2OUiKiuO5vJN1K4fzNLS21zczM8lVp99cqoHg1Vy1wd4kyq4G5kkambqy5wOrUbbZf0qzUVbWouH5KUEUXAIeOSCQdAXyYzPkUSQMljUrzg4DzsuuYmVnXaPVIpRVLgTskXQL8msKHPZJqgE9ExKUR0SjpOmBDWufaiGhM81cC3wKGAvelCeCGdKlwAE8Cf5bZ55nAjoh4IhMbAqxOCWUAsAa4ucK2mR0mOnkwPrPerqKkEhHPAWeXiNcDl2aWVwArypSbWiLe0vmQdcCsZrHfAdPbUXUzs37KoxSbmVkr3jz2mO6uAuCkYtYu6uQRXs066vt/eWYbS3qUYrMew+dUzFpW6Yl6szY54bijqDmhqrurYWadzEnFusSPPze7u6tgZl3A3V9mZpYbJxWzdvBNusxa5qRiZma5cVIxM7Pc+ES9WTuUv6GCWff46edn8+LLr3V3NQ5xUjFrB59TsZ5mQtVR3V2FJtz9ZWZmuXFSMTOz3DipmJlZbpxUzMwsN04qZmaWm4qTiqQqSXWStqbHkWXK1aYyWyXVZuLTJT0iaZukG9OthYvPfUrSY5K2SLohE786lX9c0rxMfH6KbZO0uNK2mZlZ++RxpLIYWBsRk4G1abkJSVXAEmAmMANYkkk+NwGXAZPTND+tMxtYAJwWEacCX0nxU4CFwKmp7D9LGiBpAPBPwLnAKcCFqayZmXWRPJLKAmBlml8JnF+izDygLiIaI2IfUAfMlzQOGB4R6yMigFsy618BLI2IAwARsSezv9sj4kBE/ArYRiFRzQC2RcQTEfEycHsqa2ZmXSSPpDImInal+WeAMSXKjAd2ZJZ3ptj4NN88DjAFOEPS/ZJ+LOn0NmyrVPwwki6XVC+pfu/eva21z+wQ//bRrGVt+kW9pDXA2BJPXZNdiIiQlNf/3UCgCpgFnA7cIWlSHhuOiGXAMoCamhp/TpiZ5aRNSSUi5pR7TtJuSeMiYlfqztpTolgDcFZmuRpYl+LVzeINaX4ncFfqFntA0kFgVHp+Qpl1ysXNzKwL5NH9tQooXs1VC9xdosxqYK6kkekE/Vxgdeo22y9pVrrqa1Fm/e8BswEkTQEGA8+m/S2UNETSRAon9x8ANgCTJU2UNJjCyfxVObTPzMzaKI+kshQ4R9JWYE5aRlKNpOUAEdEIXEfhg38DcG2KAVwJLKdwwn07cF+KrwAmSdpM4aR7bRRsAe4AfgF8H/hkRLwWEa8CV1FIYI8Cd6SyZrk5/cSSV8ybWVLxKMUR8Rxwdol4PXBpZnkFhURRqtzUEvGXgY+W2ef1wPUl4vcC97aj+mbt8uGaCfzVnY90dzXMeiz/ot6sHeQbqpi1yEnFzMxy46RiZma5cVIxM7PcOKmYmVlunFTMzCw3TipmZpYbJxUzM8uNk4qZmeXGScXMzHLjpGJmZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlhsnFTMzy01FSUVSlaQ6SVvTY8nb4kmqTWW2SqrNxKdLekTSNkk3KnOzCkmfkvSYpC2SbkixcyRtTOtslPTeTPl1kh6XtClNb6ikbWZm1n6VHqksBtZGxGRgbVpuQlIVsASYCcwAlmSSz03AZRTuMz8ZmJ/WmQ0sAE6LiFOBr6TyzwJ/GBFvBWqBf222u4siYlqa9lTYNjMza6dKk8oCYGWaXwmcX6LMPKAuIhojYh9QB8yXNA4YHhHrIyKAWzLrXwEsjYgDAMUEEREPRsTTqcwWYKikIRW2wczMclJpUhkTEbvS/DPAmBJlxgM7Mss7U2x8mm8eB5gCnCHpfkk/lnR6ie1+EPh5MfEk30xdX1/IdqU1J+lySfWS6vfu3dtiA83MrO0GtlZA0hpgbImnrskuRERIihzrVQXMAk4H7pA0KR3RIOlU4O+BuZl1LoqIBknHAHcCF1M4+jlMRCwDlgHU1NTkVWczs36v1aQSEXPKPSdpt6RxEbErdWeVOo/RAJyVWa4G1qV4dbN4Q5rfCdyVksgDkg4Co4C9kqqB7wKLImJ7pp4N6fE3km6lcP6mZFIxM7POUWn31yoKJ8xJj3eXKLMamCtpZDpBPxdYnbrN9kualbqqFmXW/x4wG0DSFGAw8KykY4F/BxZHxM+KO5A0UNKoND8IOA/YXGHbzMysnSpNKkuBcyRtBeakZSTVSFoOEBGNwHXAhjRdm2IAVwLLgW3AduC+FF8BTJK0GbgdqE1HLVcBJwNfbHbp8BBgtaSHgU0UjnhurrBtZmbWTq12f7UkIp4Dzi4RrwcuzSyvoJAoSpWbWiL+MvDREvEvA18uU53pba64mZl1ioqSill/dOcV7wTKXlxo1q85qZi10/QTqrq7CmY9lsf+MjOz3DipmJlZbpxUzMwsN04qZmaWGycVMzPLjZOKmZnlxknFzMxy46RiZma5cVIxM7PcOKmYmVlunFTMzCw3TipmZpYbJxUzM8uNk4qZmeWmoqQiqUpSnaSt6XFkmXK1qcxWSbWZ+HRJj0jaJunGdFvh4nOfkvSYpC2SbkixEyW9lLnr4zfasi0zM+salR6pLAbWRsRkYG1abkJSFbAEmAnMAJZkks9NwGXA5DTNT+vMBhYAp0XEqcBXMpvcHhHT0vSJTLzktszMrOtUmlQWACvT/Erg/BJl5gF1EdEYEfuAOmC+pHHA8IhYn+4/f0tm/SuApRFxACAi9rRUiVa2ZWZmXaTSpDImInal+WeAMSXKjAd2ZJZ3ptj4NN88DjAFOEPS/ZJ+LOn0TLmJkh5M8TMy+yi3rcNIulxSvaT6vXv3ttJEMzNrq1ZvJyxpDTC2xFPXZBciIiRFjvWqAmYBpwN3SJoE7ALeGBHPSZoOfE/Sqe3deEQsA5YB1NTU5FVnM7N+r9WkEhFzyj0nabekcRGxK3VBleqmagDOyixXA+tSvLpZvCHN7wTuSl1ZD0g6CIyKiL1AsUtso6TtFI5qWtqWmZl1kUq7v1YBxau5aoG7S5RZDcyVNDKdoJ8LrE7dZvslzUpXai3KrP89YDaApCnAYOBZSaMlDUjxSRROyD/RyrbMzKyLVJpUlgLnSNoKzEnLSKqRtBwgIhqB64ANabo2xQCuBJYD24DtwH0pvgKYJGkzcDtQm45azgQelrQJ+DfgE23YlpmZdZFWu79aEhHPAWeXiNcDl2aWV1BIFKXKTS0Rfxn4aIn4ncCdZepScltmZtZ1/It6MzPLjZOKmZnlxknFzMxy46RiZma5cVIxM7PcOKmYmVlunFTMzCw3TipmZpYbJxUzM8uNk4qZmeXGScXMzHLjpGJmZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlpuKkoqkKkl1kramx5FlytWmMlsl1Wbi0yU9ImmbpBvTrYCLz31K0mOStki6IcUukrQpMx2UNC09t07S45nn3lBJ28zMrP0qPVJZDKyNiMnA2rTchKQqYAkwE5gBLMkkn5uAyyjca34yMD+tMxtYAJwWEacCXwGIiG9HxLSImAZcDPwqIjZldndR8fmI2FNh28zMrJ0qTSoLgJVpfiVwfoky84C6iGiMiH1AHTBf0jhgeESsT/efvyWz/hXA0og4AFAmQVxI4f71ZmbWQ1SaVMZExK40/wwwpkSZ8cCOzPLOFBuf5pvHAaYAZ0i6X9KPJZ1eYrsfAW5rFvtm6vr6QrYrrTlJl0uql1S/d+/eso0zM7P2GdhaAUlrgLElnromuxARISlyrFcVMAs4HbhD0qR0RIOkmcCLEbE5s85FEdEg6RjgTgrdY7eU2nhELAOWAdTU1ORVZzOzfq/VpBIRc8o9J2m3pHERsSt1Z5XqpmoAzsosVwPrUry6Wbwhze8E7kpJ5AFJB4FRQPGwYiHNjlIioiE9/kbSrRTO35RMKmZm1jkq7f5aBRSv5qoF7i5RZjUwV9LIdIJ+LrA6dZvtlzQrdVUtyqz/PWA2gKQpwGDg2bR8BPBhMudTJA2UNCrNDwLOA7JHMWZm1gUqTSpLgXMkbQXmpGUk1UhaDhARjcB1wIY0XZtiAFcCy4FtwHbgvhRfAUyStJlC8qgtdn0BZwI7IuKJTD2GAKslPQxsonDEc3OFbTMzs3ZqtfurJRHxHHB2iXg9cGlmeQWFRFGq3NQS8ZeBj5bZ5zoK51qysd8B09tXezMzy5t/UW9mZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlhsnFTMzy42TipmZ5cZJxczMcuOkYmZmuXFSMTOz3DipmJn1I0MGFj72Bx5R9pZTFalo7C8zM+tdPjt3CkMGHsEFb69uvXAHOKmYmfUjxxw5iKvf95ZO2767v8zMLDdOKmZmlhsnFTMzy03FSUVSlaQ6SVvT48gy5WpTma2SajPx6ZIekbRN0o3p1sJI+o6kTWl6UtKmzDpXp/KPS5qXic9PsW2SFlfaNjMza588jlQWA2sjYjKwNi03IakKWALMBGYASzLJ5ybgMmBymuYDRMRHImJaREwD7gTuSts6BVgInJrK/rOkAZIGAP8EnAucAlyYypqZWRfJI6ksAFam+ZXA+SXKzAPqIqIxIvYBdcB8SeOA4RGxPt2D/pbm66cjlw8Dt2X2d3tEHIiIX1G4v/2MNG2LiCfS7YhvT2XNzKyL5JFUxkTErjT/DDCmRJnxwI7M8s4UG5/mm8ezzgB2R8TWNmyrVPwwki6XVC+pfu/eveXaZWZm7dSm36lIWgOMLfHUNdmFiAhJkUfFMi7k9aOUXETEMmAZQE1NTd71NTPrt9qUVCJiTrnnJO2WNC4idqXurD0lijUAZ2WWq4F1KV7dLN6Q2fZA4I+A6c22NaHMOuXiZW3cuPFZSb9urVwZo4BnO7huT9eX2wZ9u31uW+/U29p2QqlgHr+oXwXUAkvT490lyqwG/i5zcn4ucHVENEraL2kWcD+wCPifmfXmAI9FRLaLbBVwq6R/AI6ncHL/AUDAZEkTKSSThcCftFb5iBjd5pY2I6k+Imo6un5P1pfbBn27fW5b79RX2pZHUlkK3CHpEuDXFE6qI6kG+EREXJqSx3XAhrTOtRHRmOavBL4FDAXuS1PRQpp1fUXEFkl3AL8AXgU+GRGvpX1eRSGBDQBWRMSWHNpnZmZtpMJFV9YRfeWbRSl9uW3Qt9vntvVOfaVt/kV9ZZZ1dwU6UV9uG/Tt9rltvVOfaJuPVMzMLDc+UjEzs9w4qZiZWW6cVDqgtw5cKWmFpD2SNmdiJQcEVcGNqY0PS3pHZp2Sg4N2J0kTJP1I0i8kbZH0Fyne69sn6UhJD0h6KLXtb1N8oqT7Uxu+I2lwig9Jy9vS8ydmtlVyMNbulsbve1DSPWm5L7XtSRUGzd0kqT7Fev37sqyI8NSOicLlytuBScBg4CHglO6uVxvrfibwDmBzJnYDsDjNLwb+Ps2/j8Ll3QJmAfeneBXwRHocmeZH9oC2jQPekeaPAX5JYWDRXt++VMdhaX4Qhd90zQLuABam+DeAK9L8lcA30vxC4Dtp/pT0fh0CTEzv4wHd/bdLdfsMcCtwT1ruS217EhjVLNbr35flJh+ptF+vHbgyIn4CNDYLlxsQdAFwSxSsB45VYcSEkoODdnrlWxERuyLi52n+N8CjFMZ+6/XtS3X8bVoclKYA3gv8W4o3b1uxzf8GnC1JlB+MtVtJqgbeDyxPy6KPtK0Fvf59WY6TSvu1eeDKXqLcgKAVD9zZXVKXyNspfKPvE+1L3UObKAyDVEfhm/jzEfFqKpKt56E2pOdfAI6jh7YN+Efg88DBtHwcfadtUPgC8ANJGyVdnmJ94n1ZSh6/qLc+IqJTBgTtUpKGUbj/zl9GxP7Cl9iC3ty+KIwaMU3SscB3gTd3b43yIek8YE9EbJR0VjdXp7O8OyIaJL0BqJP0WPbJ3vy+LMVHKu3X0oCWvdHudHiNmg4IWq6dPbb9kgZRSCjfjoi7UrjPtA8gIp4HfgS8k0LXSPGLYbaeh9qQnh8BPEfPbNu7gA9IepJCV/J7ga/TN9oGQEQ0pMc9FL4QzKCPvS+znFTabwNp4Mp0RcpCCoNc9lbFAUGh6YCgq4BF6WqUWcAL6XB9NTBX0sh0xcrcFOtWqV/9X4BHI+IfMk/1+vZJGp2OUJA0FDiHwjmjHwF/nIo1b1uxzX8M/DAKZ3tXAQvTFVQTeX0w1m4TEVdHRHVEnEjhf+mHEXERfaBtAJKOlnRMcZ7C+2kzfeB9WVZ3XynQGycKV2j8kkK/9jXdXZ921Ps2YBfwCoU+2Uso9EevBbYCa4CqVFYUbs+8HXgEqMls5+MUToRuAz7W3e1KdXo3hb7rh4FNaXpfX2gf8DbgwdS2zcAXU3wShQ/ObcD/BYak+JFpeVt6flJmW9ekNj8OnNvdbWvWzrN4/eqvPtG21I6H0rSl+HnRF96X5SYP02JmZrlx95eZmeXGScXMzHLjpGJmZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlpv/DyYsmbcWszZuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skip = 1\n",
    "Umag = np.zeros((len(UU),m))\n",
    "E1 = []\n",
    "E2 = []\n",
    "for j, U in enumerate(UU[::skip]):\n",
    "    V = U[::2]\n",
    "    W = U[1::2]\n",
    "    Umag[j,:] = np.sqrt(V**2+W**2)\n",
    "    E1.append(dx*np.sum(V**2+W**2))\n",
    "    E2.append(np.sum((np.diff(V)**2+np.diff(W)**2)/dx - 0.5*q*dx*(V[1:]**2+W[1:]**2)**2))\n",
    "\n",
    "plt.plot(E2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "a5263a23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc084affaf0>]"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEDCAYAAAAsr19QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApAklEQVR4nO3deXxddZ3/8dcn+9J0S9KFlpK2tEDZClQEFEU2WRxQXKaM68jPoo7joI4jjAvjwCju4rgWdURnWERBkH1VkKUlLd0X6L63ado0S5v1fn5/3JP0Jk2am/be3HNu3s/H4z5yzvecc8/3m9z7yfd8v9/zPebuiIhI9ORkOgMiInJkFMBFRCJKAVxEJKIUwEVEIkoBXEQkohTARUQiatADuJn92sx2mdmyFL1fh5ktCl4PDeC4E83sZTNrMbN/Pcx+F5rZQjNbZmZ3mllekD7KzB4wsyVmNt/MTgnSi4L1xWa23My+nvBek81snpmtMbN7zawgSC8M1tcE26sSjrkpSF9tZu9MSL8sSFtjZjem4xwiEnLuPqgv4G3AmcCyFL1fYxL7bOglbQzwJuC/gH/t47gcYDMwPVj/T+C6YPk7wM3B8onAM8GyAcOC5XxgHnBOsP57YHaw/HPgU8Hyp4GfB8uzgXuD5RnAYqAQmAysBXKD11pgClAQ7DMjlecY7M+FXnrpNfDXoNfA3f15YE9implNNbPHzWyBmb1gZicOQj52ufurQNthdisHWt399WD9KeC9wfIM4NngvVYBVWY21uMag33yg5ebmQEXAn8Itt0JvDtYvjpYJ9h+UbD/1cA97t7i7uuBNcDZwWuNu69z91bgHuDqFJ9DREIuLG3gc4F/dvezgH8FfjqAY4vMrNrMXjGzd6c4X7uBPDObFay/Dzg2WF4MXANgZmcDxwETg/VcM1sE7AKecvd5xP8Z1Ll7e3D8FmBCsDyBeE2fYPu+YP+u9B7H9JWeynOISMjlZToDZjYMOA+4L14hBOKX85jZNcSbLXra6u6dbbXHuftWM5sCPGtmS919rZn9BHhLsM8xQUAFuM/d/yuZvLm7m9ls4AdmVgg8CXQEm28Dbg/edynwWuc2d+8AZprZSOCBoH18RzLnFBFJVsYDOPGrgDp3n9lzg7vfD9x/uIPdfWvwc52Z/QU4A1jr7v/UuY+Zbejt/ZPh7i8D5wfvcykwPUivB/4xSDdgPbCux7F1ZvYccBnwPWCkmeUFNeCJwNZg163Ea/Zbgk7SEUBtQnqnxGN6S69N8TlEJMQy3oQSBML1ZvZ+iAdDMzs9mWODkSCdtfUK4jXuFanMn5mNCX4WAl8i3jGImY3sHOEB/D/geXevN7PKoOaNmRUDlwCr3N2B54g3wwB8FHgwWH4oWCfY/myw/0PA7GAEyWRgGjAfeBWYFow4KSDeKflQis8hImE32L2mwN3AduKdh1uA64iPfniceLvyCuBrSb7XecSbLxYHP6/rY78NvaSNC85fD9QFy8ODbY8CxwTL3wFWAquBGxKOPxd4PUi/HxgVpJ9GvDllCbAssSzER43MJ95ReB9QGKQXBetrgu1TEo75MvGRIauByxPSrwjOvxb4cjrOoZdeeoX7Ze6aTlZEJIoy3oQiIiJHZlA7MSsqKryqqmowTykiEnkLFizY7e6VPdMHNYBXVVVRXV09mKcUEYk8M9vYW7qaUEREIkoBXEQkohTARUQiSgFcRCSiFMBFRCJKAVxEJKIUwEVEIkoBXEQkjR5bup3axpa0vLcCuIhImuxpauVT/7eQ6+5Mzw2MCuAiImnS3hEDYGvdgbS8vwK4iEhEKYCLiERUUgHczD5nZsvNbJmZ3W1mRWb2GTNbY2YePA1HREQGUb8B3MwmAJ8FZrn7KUAu8Ud4vQhcDPQ6S5aIiKRXstPJ5gHFZtYGlADb3P01gIQnyYuIyCDqtwbu8ae+fxfYRPxZlvvc/clkT2Bmc8ys2syqa2pqjjynIiLSTTJNKKOAq4k/ePgYoNTMPpTsCdx9rrvPcvdZlZWHPFBCRESOUDKdmBcD6929xt3biD+B/bz0ZktERPqTTADfBJxjZiUWb/C+CFiZ3myJiEh/kmkDnwf8AVgILA2OmWtmnzWzLcBEYImZ/TKtORURkW6SGoXi7jcDN/dI/lHwEhGRDNCdmCIiEaUALiISUQrgEmq7G1t4cc3uTGdDJJQUwCXU/v4XL/PBX87LdDZEQkkBXEJtbU1TprMgEloK4CIiEaUALiISUQrgIiIRpQAuIhJRCuAiIhGlAC6R4O6ZzoIIre0x1u8Oz8goBXARkSR95U9Lecd3/8KeptZMZwVQABcRSdpLa2sBaGppz3BO4hTAJRLUgiJyKAVwEZEBCkuFQgFcRCRJZpnOQXcK4CIiSQpLzbtTUgHczD5nZsvNbJmZ3W1mRWY22czmmdkaM7vXzArSnVkZukL2vZEhLiw18X4DuJlNAD4LzHL3U4BcYDbwLeAH7n48sBe4Lp0ZFREJi7DUxJNtQskDis0sDygBtgMXEn/YMcCdwLtTnjsRkRAJS827UzJPpd8KfBfYRDxw7wMWAHXu3jkYcgswobfjzWyOmVWbWXVNTU1qci1Dju7EFDlUMk0oo4CrgcnAMUApcFmyJ3D3ue4+y91nVVZWHnFGRUSku2SaUC4G1rt7jbu3AfcDbwFGBk0qABOBrWnKo4iI9CKZAL4JOMfMSszMgIuAFcBzwPuCfT4KPJieLIqISG+SaQOfR7yzciGwNDhmLvAl4PNmtgYoB36VxnzKEKcWcAkTD8knMq//XcDdbwZu7pG8Djg75TkSEZGk6E5MEZEBMsIxnlABXCJBowhFDqUALiIyQGFpA1cAl9DaVd+c6SyIdNPZdLK7UU/kETms1o5YprMg0k1nzfu9P3spFHcHK4BLJITlklUkTBTARUSSFJbRJ50UwEVEjkAIWlAUwCUawvBlEQkbBXAJLQvb5MsiIaMALiKSpL1NB4cPhuGiUAFcQisMw7REEjW0tPe/0yBSABcROQJhqGAogEtoqQ1c5PAUwEVEIkoBXETkCGS+AUUBXCIiBM2NIqGTzFPpTzCzRQmvejO7wcxON7OXzWypmf3ZzIYPRoZFRCQumWdirnb3me4+EzgL2A88APwSuNHdTw3Wv5jOjMrQoy5MCbMwXBUOtAnlImCtu28EpgPPB+lPAe9NZcZEEmk2QpFDDTSAzwbuDpaXA1cHy+8Hju3tADObY2bVZlZdU1NzZLmUIUkhW+Twkg7gZlYAXAXcFyR9HPi0mS0AyoBeH1Hh7nPdfZa7z6qsrDza/IqIhEIYrgrzBrDv5cBCd98J4O6rgEsBzGw6cGXqsyciIn0ZSBPKtRxsPsHMxgQ/c4CvAD9PbdZkqEvsxAxDh5FI2CQVwM2sFLgEuD8h+Vozex1YBWwD/if12RMRCacwVCqSakJx9yagvEfa7cDt6ciUiIj0T3diSiSEoLIjEjoK4CIiEaUALqGl2WRFDk8BXEIrDJ1EIn0Jw+dTAVwiIQxPPxEZqHR/arM2gNc2tvDhX81jT1OvN4iKiERe1gbw37y0gRfe2M3/vrIx01kRkSyUzK306e7GydoALtGX2ImpBhSRQymAi4hElAK4iMgRCEO/ugK4iEhEKYBLaFlCF1AYajsiYaMALqEVhgnzRfoShk+nAriISEQpgEs0hKG6I5IgDHcHK4CLiESUAriElqX9PjaRaOs3gJvZCWa2KOFVb2Y3mNlMM3slSKs2s7MHI8MiImGQ+QaUJB6p5u6rgZkAZpYLbAUeAO4Avu7uj5nZFcC3gQvSllMZ0jQiReRQA21CuQhY6+4bif8DGh6kjyD+YGMRERkkST3UOMFs4O5g+QbgCTP7LvF/BOf1doCZzQHmAEyaNOnIcikiEjIhGISSfA3czAqAq4D7gqRPAZ9z92OBzwG/6u04d5/r7rPcfVZlZeXR5leS8B8PLWfWrU9nOhspFYYvi0jYDKQJ5XJgobvvDNY/CtwfLN8HqBMzJH7z0gZ2N7ZkOhsikmYDCeDXcrD5BOJt3m8Pli8E3khVpkREQi8EV4VJtYGbWSlwCXB9QvIngNvNLA9oJmjnHizPrNzJmZNGMaq0YDBPKyISGkkFcHdvAsp7pP0NOCsdmerP3qZWrruzmjdVjeK+T/badypZJgSVHZHQieSdmG0dMQA21O7PcE4knUw3YkqIheHehEgGcBERUQCXiAjDzG8iicLwkVQAl9AKwxdEJMwUwEVEIkoBXEJLnZgSZmG4QFQAl0gIw5dFJGwUwEVEIkoBXETkCIRhZJQCuERCCL4rIqGjAC6hpT5MScY/3/0aVTc+kulsZIQCuIhE2p8XZ+ZhYGG4KFQAl9AKwxdEJMyyPoCr7TQ7hGHiIJFEycSWdH9qsz6AS3SpDVzCJBYLXyUi6wO47uYTkVR46Aja2tMdfrI+gEuWCF/lR4aY+ua2buthaNbr94k8ZnYCcG9C0hTga8C5wAlB2kigzt1npjh/MoRl/ushEm79BnB3Xw3MBDCzXGAr8IC7/7BzHzP7HrAvPVmUoUqtXyKHN9AmlIuAte6+sTPBzAz4AN2fWC9DVFNLO5ff/gLLth79/3NTB4aEWQguEQcawGdzaKA+H9jp7m/0doCZzTGzajOrrqmpOZI8SoRUb9zLyu31fOvxVSl93xB8V0RCJ+kAbmYFwFXAfT02Xcthat/uPtfdZ7n7rMrKyiPLpQxJYegkEgmzftvAE1wOLHT3nZ0JZpYHXAOcleqMiYiESc8bd8JQvRhIE0pvNe2LgVXuviV1WRI5lO6olf6s3tHA3fM3ZTobgyqpGriZlQKXANf32NRbm7iIyKB75w+fB+Dasyel5f3D2KeeVAB39yagvJf0j6U6QyIiURCGq0LdiSnSi92NLezb39b/jjJkhCFg96QALpEw2CNSZt36NGfc8mTX+p0vbeDy218AYGNtE999YjVVNz7CL19Y1+vxv/rbet7Y2dC1vrG2KRSP4JLUCcMoKQVwkT4kTj5380PLWbm9HoAP/nIeP35uDQB3zeu90+yWh1dw1Y9fBOC1TXt5+3f+wv++srHXfSW1ttUdSMv7hrENXAFcwivzFZxe7W/tSGq/A23x/dbvbgJg4aa6oz73I0u209yW3PmHqvNuezbTWRg0CuASCVFqfUhXU8kr62r5p7sW8o1HV3ZLv/rHf+O637zatf7imt1c95tXQzl/dZQdMg48BL/egdzIIyIZVH8g3qm6ra65W/riLd3nnZnz22qaWjvY39bBsEJ9xbNZ1tTAv/qnZUP2ydQSLpmumYWgYpiVeraBh+H3nDUB/HfqIMpq+vtK2Pw06MjOpKwJ4JLdfvHXtZnOQtLCUDOT1Ot5ZfXkip297ziIFMAlEtQfJ5mW7OijwaQALpJiumEnO62raey2nsyfOd2fBAVwEZGIUgCX0Mrmemws5rR3xDKdDRmAI/k8pvvmTQVwkcCeplZqG1uO+n2S+aJ/8Q9LOP7Ljx31uWRoUwAXCZx5y1OcdevTg3KuPy48+meg7G1q7XbvwyNLth/1e6bDR349Pyvu0Ti0Np35a0QFcJEUG6w+zLU9OtX+6a6Ffe5750sbqLrxEQ5kYCTF869n58PMw9BXHckAXr1xLwA1DUd/uSsSPQOPHD8PxtHv3d+a6swMGSGI14foN4Cb2QlmtijhVW9mNwTb/tnMVpnZcjP7dtpzG7jzpQ19bvvW46uob9ZE/JI+/XVMhWGe6KFiqA/Z7HemG3dfDcwEMLNcYCvwgJm9A7gaON3dW8xsTDozmqyf/WUt+w60UV5akOmsiKTUwVAVwompB9kr62opzs/tmrK3p9rGFv64cAufOH8KlqKJvHu+Sxj+dQy0CeUiYK27bwQ+Bdzm7i0A7r4r1Zk7Un1Nsi8Dt2hzHat3NPS/o3Tpq1KY6tpif++WzbXT2XNf4eqfvNhnAP/c7xfzjUdXsaTHTI1HI4y/zYEG8MSn0E8HzjezeWb2VzN7U28HmNkcM6s2s+qamuzszMhm7/7Ji11P+xaJioagGbU9y+dgSDqAm1kBcBVwX5CUB4wGzgG+CPzeerlWcfe57j7L3WdVVlamIMsyVGRbBTJVl/Jd79fP9iz79UkvBlIDvxxY6O6dU3BtAe73uPlADKhIdQZFpHcK0JkVhiaqgQTwaznYfALwJ+AdAGY2HSgAdqcsZyJyVEIQXzJmqJQ9qQBuZqXAJcD9Ccm/BqaY2TLgHuCjPgj/klrbY+w7cHCY4OPLwnn3mQxdQyV4DLaFm/YO+N6PMD5JPpWSCuDu3uTu5e6+LyGt1d0/5O6nuPuZ7j4oj4L+5P8uYFXCqIjP3rNoME4r0qtMxuosj02HuOanL3Hlj14Y0DHp/Gcahv/TkbsT89lV3UcrtrZrRrdstaG2KdNZOCLpupFnwOOQwxBhUmxXkjXwwah5h+EfaOQCuAwds+e+kuks9Otov8SD3RGWbTE9k0E0DL9LBXCRFAtLG/hQuKW/rxIOxt8gDH9nBXCRDDqaIHAktc8wXPYPpmzvxOx3LhSRKNtY28S2umbaOmKcP60i5TfT9KbvWuGhW46mEtf/rfRH8ebSrzCMA1cAl6zU2h7j1kdW8NuXNx6ybUpFKcOK8mhqaWffgTZmHTea4cUHvwpzn1/btfzE8h20dxz8ou7Y10xTa3vX+tF+heNBIMuriWk01H9zCuCSdRqa2/jwr+azaHMdHzuvihnjh/PnJduo3rCXsqI82mPOyJICxpQVMW99LQs27SU3oWb+jUdXdS1f/7sF3d77nG8+0219/e4mTv7a4xxXXkp+rlGUn9tt/o3P3LWwa0KlPy/ZzrlTyykryu/avqF2P2OGF9Le4Ywozqe1PUaHO0V5OeTlxls4O2JObs7AQ1Xm64eZl8pK8h8WHP1TlFItUgG8TpPRSxK+/fhqFm2u4/bZM7l65gQAPvCmY5M6tiPmNLa00xFzciweoGPubN/XzLKt9bR3xKgsK2TJln3saWqlfFgBS7bso7aphbycHCrKCrtdWj+c8JizjpjzpT8u7Xa+i7//1z7zUpCbQ2vw4ONTJ4xg+74DADy9cifX/66615tanly+g/3BU3fe2NlA+bBCjh8zLKmyy8CE4R9kpAL4W7/1XKazICG3Y18zd83fxIfPOa4reA9Ebo4xovhgDfmMSQfnlX/XaccM+P3cnZb2GHk5RkNzO02t7dTtbyM3x2hu62Dx5jqaWjsoKchlb1MruTk5rN/dSFVFKQ3N7ezY18wjS7dTnJ/L5IpSdjfGKzFPLN/Z6/nmJFwx/H0wDLOsMI+Glnizz1U/fpGKYQW4Q3FBLvsOtFGUn0vFsAL+7rRjKCrIZeveA5wzZTSjSwsYN6KIwrzcAZc7LNLZ5ZHMW6c7yEcqgDe2tPeavmaX5qs+nMWb6/h99Waa22JMHFVMSUEuJQW5lBXlM7IknzOOHcWIkvz+3ygJme7Yuf+1LXTEnE+cPyWj+ehkFm9WARhVWsCo0gImjjq4/YxJo/o48qCfHGabu9PU2sHuhhZ2N7awd38bLe0dvL6jgdGlBby6YS+VZYU0tbRz34ItzDpuFO3B1cW+A200t3Wwfnf8hqkX3uh9KqOzq0YzcXQxw4vyqW1qpbQglxHF+Vx68liOH1PG8KK8lHUOd8Scl9buZntdM8OK8iguyKUwN4dpY8u69rn14RVdyz/7y9re3oZFm+sAqG1spaG5revhDzvrm3lkyQ4Wb6ljd2MLdfvbyM81SgvzKC3Io7Qwj8qyAqrKS1mydR+PLNnO67deTkHeoQP2VAMfgNhh5vW9+Puar7ovDy3exufvXURujjGyJJ+d9b3fyVacn0txQS7DCvMYN6KIscOLGF7UvWMvLyeHy08dx/gRxYOV/QGbt24PJ4wtY1J5SaazMijMjGGFeQwrzKOqovTghtPiPz72lsldSd95/+l9vk9Leweb9xygtT3Gy+tq2VZ3gJfW1rKzvpk1NY1s2rOf+uY2DrR1dLUr/+L5dQAMK8xj/IgicnOsqwnq2FElnDZxBFMTmm9eXLOb48pLKCnI40BbB+OHF5ETtO1v3rOfl9fV8uu/re82VUZvfvfKwY7peev3HHbfT/y2GjMozMuhuS0W/M7ghLFlVJYVMrmilLaOGPtbO2hqaWfL3v0s3LSXPU0Hm2vP/eYz3D3nnEPfPAQRPDIB/OGlmrRqoJZv28cN97zGrKrRzP3wWYwsKWDfgTbMoKUtPinY9n0HWLGtnu37mtkfXN7XNrWydEtd1+U6HOzYu+2xVVwyYyxjhhfiHv9imBmFeTmUFuYyb93hv1DpFIs5r23ay5VH0NQx1BXm5Xa1lc84Zvhh921tj7GnqZVFm/fy3Koa1u1upLggj/wcoy3m5BqsrWnitZf30pYwgueDv5x3yHtNGl1Cc1tH1y3yUypK+fZ7T2PGMcOpaWyhsbmdbXUH2LhnP+dNLefsyaMZU1aEu9PaESMWg9aOWFezV3NbBzUNLSzfto/7qrdw/rQK6g60saeplQkjiykryueik8YwdnjRYcu4b38bV/zoBbbWHaC2qZVLf9BLJTEEQ2AiE8Ab9KDiAfv6n1cwsqSAOz48q6uJpKt9twgqy+IdXOdP6/9BG3uaWqlpaOHOlzfw9Iqd8UdZObR0xHD3bl9UgPoDg//3WlvTSH1zO2dOGjno5x5KCvJyGDeiiMtGjOeyU8b3uV9HzNlWd4D65jZWbm9geFEeO+ub6Yg5r6zbwxu7GphcMYxRJflMH1vGW6dVMH1sWVIjbuKVhnjTVDEH2+iL8nM5dnQJx44uOWze+jOiJJ/zp1Vwz6ub+95JNXBJp/nr9/CVK09KSfv26NICRpcW8I33nMo33nPqIdtjMedAWwcvra3lE7+tprZp8EcMLdi4F4Czjuu/XVnSLzfHOHZ0vCnr5GNGdNuW2LQTdt+85lTueGEd62q6T64WgvitW+mz3d+dPjjNCTk58Y6gS2aM5eNvmcyWvQdYHHQkDZaFm/YyqiSfyYltwSIp8PlLph+SFoIWFAXwbHb8mGH9tvWlww2XTKMoP4e7528a1PMu2LiXMyeNGpTb5WVoeefJ4w5JG0gNPF2DsyITwEMw7UDknDZxRP87pcHwonyuOGU8jy7dTkt7x6Ccs25/K2trmjhTzSeSIu85I34fwTlTysnPPbpQubuxhedfr0lFtrrpN1dmdoKZLUp41ZvZDWb2H2a2NSH9ipTnTo7KzGNHZuzcV808hvrmdv66OvUf2t68FjTXnKEOTEmRN08pZ8NtV/bZJDfQex4qhhWmIlvd9BvA3X21u89095nAWcB+4IFg8w86t7n7oynPnRyVGeMPPxwsnd5yfAXlpQU8uHjboJxv8eY6zOC0iSMH5XwiAzWlMvV9MwMdhXIRsNbdNw5mO+NfVu/ioUWDEwiiom5/K/e8upm29hg76puJuTNhZPcbbKZUZm4OjPzcHK48bTz3vrqZnzy3BjMwjByDHLP4usXXjXgnqJnFl4P0nh5ftqPXu27d4YdPv8HwovgNLSKDbcW2ep5ddej0Bjvqm7uWO+/ITaWBftpnA3cnrH/GzD4CVANfcPe9PQ8wsznAHIBJkyYdUSafXLGT+Rsyd4NIGD22bAe3PbbqsPuMLi047PZ0+8CsY7ln/ma+88TqlLzfZ+957bDPQL38KMb9igzUyJKD36/bn3m9z/lpIH1Xw0kHcDMrAK4CbgqSfgbcQrwz9hbge8DHex7n7nOBuQCzZs06oq7I/3r3KXz1yhnUNLSwpqaBPyzYwqNLdxz2mGzv9GwPZqmb9+8XUV5acEiPeN4RTD+aaqdMGMHKWy6jI+bE3HGPP+Yr5hxc9/h64k8nvv32p9/odiNFa3uML1wynU9eMPWQcxl0Tb8qkg6rbrmME7/6eNd6ZdnBNu3mthinThjB/Z8+75Dj0vnZHEgN/HJgobvvBOj8CWBmdwAPpzhvXcyM4oJcJpWXMKm8hAtPHNu17aW1u/mHO+YxfewwRhTn8+qG+EXAj59bA0T3yeb96ZwapiA3J9SBKzfHjmgua4Db3nsa/3bZiZx5y1OcNH44K7fXM6Ik/6hHBIgciaL8XDbcdiUQn1LijhfWsauhmTFlRbTHYhTk5Qz6Z3MgZ7uWhOYTM0u8Xn0PsCxVmRqI86ZW8LvrzubRz57PfZ88j9W3Xsb3P3A6xUF70/0Lt/L3v3iZBxdtpa2j78vvqIkFlxg5WT7meXRpAVfPPIaV2+sB1MYtofC+sybQEXP+vDg+R1Nbu5OfO/jfxaQCuJmVApcA9yckf9vMlprZEuAdwOfSkL+knD+tsqsWWpiXyzVnTmTlLZcx/8sX8aXLTmT7vmb+5Z5FXPPTl3h4ybbDtqNGRWcN3IZAZfTMhClXFcAlDI4fU8bUylJeeCM+TLYtFsvIlWFS3wZ3bwLKe6R9OC05SqExZUV86oKpXP+2KTy8dDvffHQln7nrNYYV5nHVzGN4+/RK3jatkuKC6E1Y3zm9brbXwAFOTbghaViRAriEw3lTK/jjwi20dcRo6whxAI+6nBzjqtOP4fJTxvHcql38vnoz98zfxF3zNlFakMtFJ43lUxdM5aQMjpseqINNKBnOyCA4cdzByfxVA5ewOHdqOb97ZSPLt9XT3pGZJpQh9W3Iz83h0pPHcenJ4zjQ2sFrm/fywMKtPLh4G48s3c5lJ4/jslPG8Y4Tx1BakBvqOTU6m1CGQg28pODgx1QBXMLi1AnxK8Pl2/bR2hHLyGCCIfttKC7I5bypFZw3tYIvvvMEvvX4ah5ftp1HggdHFOTmcPGMMVxx6viuSeTDpLMGPgTidzdqQpGwmDiqmOFFeSzbWk9bR4wCBfDMGDO8iO994HRue++pPLViJ2t2NbJqRz2PLt3Bo0t3UJSfw1eunEFVeSnjRhSF4infPkRGoXR631kT+cOCLQwvSs2zO0WOlplxyoQRrNi2j/YOz8i9FwrgCfJzc7ji1IOjI2sbW3hk6Xa+9uByvvKng6MkJ44qxgyqyks5u2o0lWWFvOPE/h/TlEpDqQkF4LZrTuULl05Py+3IIkfq+DHDeGDhVgrzc8jv5cHH6aYAfhjlwwr5yLlVXHTSWKo37KH+QBt3vryRNbsaAdi858AhT/J+U9UoGprbMTNGFudz/duncMEJY1Ket6HUiQnxO9nC/DBlGZqqyktpaGmnoQXyVQMPpwkji5kwMz438IfPrQLiTRg761uoO9DKwo113P7M6+xqaKGhub3bU7VfXlcLxB/LVJyfy9Uzj0lJ52jXOPAhUgMXCaOqipKuZQ0jjBAzY9yIIsaNKOLEccP5hzd3n6jLPf7g1mvveAWAm+5fCsAN9y6irDCPt02v5F2njee84ysOPmh4ANx9yNS+RcKqqvzgFLFqQskiZsa5U8tZ940rqG1q5fWdDXzwl/MAaGhp55Gl8REvJ4wt47F/OZ+cAUbjmPsRzzEiIqlx7OgScix+RawmlCyUk2NUlhVSWVbIhtuuxN3Z3djKE8t38MOn32D1zgb+56UNXPfWgT2luyOm5hORTMvPzaGyrJCd9S0ZaUIZAjNphItZPKB/6JzjeOnGCwG45eEV7EqY+D0ZakIRCYdxweizTNzIowCeQQV5Ofz8Q2cCcPY3nqGmoSXpY2PuQ2YIoUiYdQ4fDu1shJI+l50ynjdVxWfbu/531UkfF/OhMwZcJMzGjYgH8Ew8REYBPATu++R5vHnyaBZuqmPqvz9KbWP/NfGY+5C7jV4kjDpr4DVJfG9TTQE8JP772jMA6Ig5Z936NGtrGg+7v6sGLhIKk0bHx4IfaO0Y9HMrgIfEmOFFbLjtSs6fVgHAVf/9N5pa2vvcP6ZOTJFQeOfJ47j+bVP43CXTB/3cCuAhc8dHZvGu08bT1NrByTc/wV9W7+p1P3ViioRDQV4ON11xEqNLC/rfOcX6DeBmdoKZLUp41ZvZDQnbv2BmbmYVac3pEFGUn8vts89gamX8Dq+P/c+rzPltddfsg51irnHgIkNdvwHc3Ve7+0x3nwmcBewHHgAws2OBS4FN6czkUJObYzzzhQu44tRxADy5YifPrupeE9c4cBEZaBPKRcBad98YrP8A+DcgAwNost9PP3gW9845B4Dr7qymue1gJ0lHTLfSiwx1Aw3gs4G7AczsamCruy8+3AFmNsfMqs2suqam5gizOXS9eUp51w0CNz+4vCtd48BFJOkAbmYFwFXAfWZWAvw78LX+jnP3ue4+y91nVVZWHnlOh7CFX70EgHurN/Off14BaBy4iAysBn45sNDddwJTgcnAYjPbAEwEFprZuNRnUcqK8jl9YvwBqg8u2oq7axy4iAwogF9L0Hzi7kvdfYy7V7l7FbAFONPdd6QhjwI8+Jm3cuu7T6G2qZU3djVqHLiIJBfAzawUuAS4P73ZkcO5+KSxADyzcpfawEUkufnA3b0JKD/M9qpUZUj6Nm5EETPGD+fZVTsZM7xIbeAiQ5zuxIyYC08cw6sb9rJiW71q4CJDnAJ4xFx52ngA1u9uoq0jluHciEgmKYBHzEnjh/P+syYCMGFUcYZzIyKZpGdiRtAX33kCI0vy+cyF0zKdFRHJIAXwCBozvIgvXzkj09kQkQxTE4qISEQpgIuIRJQCuIhIRCmAi4hElAK4iEhEKYCLiESUAriISEQpgIuIRJT1fNp5Wk9mVgNs7HfH3lUAu1OYnTDJ5rJBdpdPZYuuKJXvOHc/5JFmgxrAj4aZVbv7rEznIx2yuWyQ3eVT2aIrG8qnJhQRkYhSABcRiagoBfC5mc5AGmVz2SC7y6eyRVfkyxeZNnAREekuSjVwERFJoAAuIhJRkQjgZnaZma02szVmdmOm85MMM/u1me0ys2UJaaPN7CkzeyP4OSpINzP7UVC+JWZ2ZsIxHw32f8PMPpqJsvRkZsea2XNmtsLMlpvZvwTpkS+fmRWZ2XwzWxyU7etB+mQzmxeU4V4zKwjSC4P1NcH2qoT3uilIX21m78xQkQ5hZrlm9pqZPRysZ1PZNpjZUjNbZGbVQVrkP5d9cvdQv4BcYC0wBSgAFgMzMp2vJPL9NuBMYFlC2reBG4PlG4FvBctXAI8BBpwDzAvSRwPrgp+jguVRISjbeODMYLkMeB2YkQ3lC/I4LFjOB+YFef49MDtI/znwqWD508DPg+XZwL3B8ozgs1oITA4+w7mZ/tsFefs8cBfwcLCeTWXbAFT0SIv857LP8mY6A0n8Qc4FnkhYvwm4KdP5SjLvVT0C+GpgfLA8HlgdLP8CuLbnfsC1wC8S0rvtF5YX8CBwSbaVDygBFgJvJn7HXl7PzyTwBHBusJwX7Gc9P6eJ+2W4TBOBZ4ALgYeDvGZF2YK89BbAs+pzmfiKQhPKBGBzwvqWIC2Kxrr79mB5BzA2WO6rjKEve3BZfQbxmmpWlC9oYlgE7AKeIl7DrHP39mCXxHx2lSHYvg8oJ6RlA34I/BsQC9bLyZ6yATjwpJktMLM5QVpWfC57o4caZ4i7u5lFegynmQ0D/gjc4O71Zta1Lcrlc/cOYKaZjQQeAE7MbI5Sw8zeBexy9wVmdkGGs5Mub3X3rWY2BnjKzFYlbozy57I3UaiBbwWOTVifGKRF0U4zGw8Q/NwVpPdVxtCW3czyiQfv/3P3+4PkrCkfgLvXAc8Rb1YYaWadFZ7EfHaVIdg+AqglnGV7C3CVmW0A7iHejHI72VE2ANx9a/BzF/F/vmeTZZ/LRFEI4K8C04Ke8gLinSkPZThPR+ohoLNH+6PE24470z8S9IqfA+wLLvmeAC41s1FBz/mlQVpGWbyq/Stgpbt/P2FT5MtnZpVBzRszKybetr+SeCB/X7Bbz7J1lvl9wLMebzh9CJgdjOSYDEwD5g9KIfrg7je5+0R3ryL+PXrW3T9IFpQNwMxKzaysc5n452kZWfC57FOmG+GT7Ji4gvhIh7XAlzOdnyTzfDewHWgj3oZ2HfH2w2eAN4CngdHBvgb8JCjfUmBWwvt8HFgTvP4x0+UK8vRW4m2NS4BFweuKbCgfcBrwWlC2ZcDXgvQpxIPUGuA+oDBILwrW1wTbpyS815eDMq8GLs902XqU8wIOjkLJirIF5VgcvJZ3xops+Fz29dKt9CIiERWFJhQREemFAriISEQpgIuIRJQCuIhIRCmAi4hElAK4iEhEKYCLiETU/wdWBrylI+vhhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(E1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "67304f0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fc0a35beb50>"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmyElEQVR4nO2df6wkV3mmn7f6Xs/ENl7bMWs5trMYNGRloqxjLONVEsQuG7CtVQyriLX/CE4WxcnGSEHKamUSaWETIWWzIZHQZo2GZYQdsSYODsGKnCXGi4IirYGBOP4BcTwGI89osAUmGDAe39v17R91qrq6bvWP2923u+70+0ilrjp1zulT03fe/vo73/mOIgJjjDHrQbbqARhjjFkeFn1jjFkjLPrGGLNGWPSNMWaNsOgbY8wasbHqAUziDB2Ig5y16mGYVSKhoesdJ0N1R56rfB2cR1knE1FWT+chIIM8g9goDm3k9Hp51W2/nxHbGdoGbUOWAzkoioO8iI4bPg+IgCC91s6pvZa0RtjF0Et16mg8A3yXb38zIl7edq/zon+Qs3id3rjqYZhFo8k/MpWpqls/pzpXrW5WXFevQr1ecTPLYGMDeul1c4PY7A2OjR75GRn5Zla8boj+pugfFNsHxUtnw0vnwakL+mxe8CLnnv1C9b7/+L0z2frmQQ58s8cZ34YzvgcbLwa9F4PeVpBtB9lLOdlWTvZSjrb7aGtwsLUN29vQz4vXPCf6/aLzPH0h5DlEEPngy6YS9zwgivJonk8i8sl1zL7k0/Hxr4+613nRN2ahqOXXQaK09ENAsvQjY+gHRQT0Ixu6HvRd1K+3r/e72/EYsxdY9M2+R9l8U1M1DR8Idin2TQEP0c8HQh3pRqhRPxu0K2tH1uqQmhpl2bC1b8wMWPSNgYFvv6S09FUT9ABC5HnjWyJG1G/0b5vedAGLvjElaYJ3SLgbAh45w5Z+3fBu/DIY9FF+Yxizeiz6Zr1p86mnCJ8hyz0RocqlU15X50NiX/p1GmJvH75ZMRZ9Y2rsmHRts/bbZmZb6rX2Z8yKmTgDJulSSZ+R9GVJj0n69VR+vqT7JT2RXs9L5ZL0AUnHJD0s6cpaXzen+k9IunnvHsuYXSBBxlC0TTMCp7wRuchjcEQ+XGko8qcevZNhK990gmnCHraB34iIy4FrgFslXQ7cBjwQEYeAB9I1wHXAoXTcAtwOxZcE8B7gdcDVwHvKLwpjukg1udvU6qgddSqRt7ib7jJR9CPiZER8KZ1/F/gKcDFwA3BHqnYH8JZ0fgNwZxQ8CJwr6SLgzcD9EfFcRHwbuB+4dpEPY05fplpstAiq1brldXMgA79+lJE7Q+1rr9Uk7t6ztH8fs+/ZlU9f0iuAnwQ+B1wYESfTrW8AF6bzi4Gna82Op7JR5cZ0g7rfvib8w5O5amQ60M46TbeQDX/TIaYWfUlnA/cA74qI51WzYCIiJC3M1JB0C4VriIOcuahujRlmghXejLdXlSJnuJ2aq3InibzdP2aFTLWUUdImheB/NCL+LBU/k9w2pNdnU/kJ4NJa80tS2ajyHUTE4Yi4KiKu2uTAtM9i9hNT5H1Zpsui9MM3Uy/sWGgVELmqg1GCr1pKBpbr53feHTOOaaJ3BHwY+EpE/EHt1r1AGYFzM/DJWvnbUxTPNcB3khvoU8CbJJ2XJnDflMqMmZ9pRHVSnSqRW/uqWgBKoS+PvKVSs3024X0XMXZjpmQa985PAb8APCLpoVT2m8DvAndLegfwdeBt6d59wPXAMeAF4JcAIuI5Sb8DfCHV++2IeG4RD2FME02RxRMKC3yk/70Ze19lPm4G44+o3+iv8PtPTsdQjj3oTxy/MbtlouhHxN8weipqR87jiAjg1hF9HQGO7GaAxiyboUnc5r2GdR+xw/tTMI1v35gV4BW5xsCwD5/h82ha+/VvgpqV35ZKue7jN6YLeLtE032ak461icqZ0yq3+MgHAq3hhGs7KtWOpjlfL65NDk/z/lMNu/68zQlbT86aKbClb0yd2ircNuEvtj1sXFcXjXZVf144ZbqDRd+cPkxr9e/InT+w7MuQzWimWSbF6eci+oP2ylUJ/3DOHhFZDOfwaXnfqZ/Lm6eYBWHRN/uLyBm7v27aH3dqxvjwh3bPKm/kgephmnmtYaNd63aJu/HqZEr9j/mlYJeO2SX26ZvVMUGwtBvxnoW6pd1IrtYq2DHmYES7ITdPy/keMPHfzV8Ua40tfXP6Mq241q39xt64rRO5eeO62Vd9njdr3JtqzJ4DMHuHLX2zb5g5LcMuff0jrfVIE7e1Q3VLf9SvhN347mfAGTbNbrDom06y566diQNovJaMce2MbbdkVv7vZzqLRd/sfxbpI5910nXRfbT2ayE382OfvjFMyIdfFqfUymokWWtLKj46H78xq8WWvjFtKHlvGoLdFHg1IneqzbQs8qaj2NI3a0891/1Yqzwar83yHf028/D4m8CsHlv65vRggTl4BvcYkXtnzHVbm2nfbxyzPp8xDWzpm9OTWaJXJvj0gUGYZu26tdpQ/p1djqMcu9Ppmz3A5oM5fZjHfTJNTP2ocPhxYfKNlb4zYbeQWSAWfbM/GbUgKdPUu2ZNZKS1r/HO/wVptJSN/sXiBVlmRqbZI/eIpGclPVor+xNJD6XjqXIbRUmvkPSD2r0P1tq8VtIjko5J+kDae9eY2RmXiGyR1DMjjJrEDZYXsbOs5zanJdP49D8C/A/gzrIgIv59eS7p/cB3avWfjIgrWvq5Hfhl4HMU++heC/zlrkdsTC3TZuQ56vUW/x6z+uP3eEVu1FMsO3GamYGJln5EfBZo3cA8WetvA+4a14eki4BzIuLBtIfuncBbdj1aY7rAqBQMxuwD5nV+/gzwTEQ8USu7TNLfSvprST+Tyi4GjtfqHE9lrUi6RdJRSUe3ODXnEI0xxpTMG7J5E8NW/kngRyPiW5JeC/y5pNfsttOIOAwcBjhH59ueMuMZ5eOecdpoUYuoYpY0ydLo57Ev3yyAmUVf0gbw74DXlmURcQoK0zwivijpSeDVwAngklrzS1KZMcaYJTKPe+ffAH8fEZXbRtLLJfXS+SuBQ8BXI+Ik8Lyka9I8wNuBT87x3madWFQI5gIYl3tn5XTo38l0l2lCNu8C/h/wY5KOS3pHunUjOydwXw88nEI4Pw78akSUk8C/Bvwv4BjwJI7cMcaYpTPRvRMRN40o/8WWsnuAe0bUPwr8+C7HZ8xi6PqykK6Pz5w2+PegOf3ZL4K6X8Zp9jUWfWOMWSOcZdOcnuxnq7kae5dmic3pgi19c3qQLyclwaR0+nvGkp7PnP5Y9I0xZo2w6BvTQAta+bqofoxZJBZ9Y4xZIyz6Zr1ZtjVu69+sGIu+MdMwacJ2HwcLmfXCIZvGwO6iI+sCP207G/imI9jSN8aYNcKib9ae1iibeS3zlvaO5jFdwKJvjDFrhEXfmHE0jfPmhG3z2sa86TgWfWPq1F0wswp4vZ1dOqZjOHrHmJJJ+lxa9fXtskKTwzWt+6ZD2NI3py9ZVhy7pSnSs4ZlziL2s4zXmF0wzXaJRyQ9K+nRWtl7JZ2Q9FA6rq/de7ekY5Iel/TmWvm1qeyYpNsW/yjGTGARe8hO8unP1KeF3iyPaf7aPgJc21L+hxFxRTruA5B0OcXeua9Jbf6npF7aLP2PgOuAy4GbUl1jVscYf/u4Dc9Dg2OW9vbzm1UyzR65n5X0iin7uwH4WEScAr4m6Rhwdbp3LCK+CiDpY6nul3c/ZGOMMbMyz+/Kd0p6OLl/zktlFwNP1+ocT2WjyluRdIuko5KObnFqjiGatWIOC7q0zIcs9LZwzaYJrxgbttna727xLwOzQGYV/duBVwFXACeB9y9qQAARcTgiroqIqzY5sMiuzToy5daJUwn/JOYR/P28xaPZN8wUshkRz5Tnkj4E/EW6PAFcWqt6SSpjTLkxe0+WBLUprEG7sI8Q6hBDln0wYi53VJ87fj1oML6+LXqz98xk6Uu6qHb5VqCM7LkXuFHSAUmXAYeAzwNfAA5JukzSGRSTvffOPmxj9oiaMO/w2mjwWp/IbU7oVpejvlCMWSETLX1JdwFvAC6QdBx4D/AGSVdQ/Ek/BfwKQEQ8JuluignabeDWiOinft4JfAroAUci4rFFP4wxM9Pwm8/qix/Zzn550xGmid65qaX4w2Pqvw94X0v5fcB9uxqdMcskGs6aUQI+VKCdxRZ802GchsGsJyOEWHkS77bbqh3Nsh39p8CefHfvb8xeY9E3ZpwAN27tmMhtpt7xoizTcSz6xiSa7pmR/vwJkZVT92PMCrDom/3HHuSq2eGPr7+W5+XirLqKi0HYZYxozx4JvzKIUf4jY9pxpifTffZQ2Nq2MKx0fZJff5I/f9lbJvoLwEyBRd+YklaR3l0XrfXt3jEdwu4ds76MsPIH91tuNrJrqsrHszPUc0d/Y97XmGVh0TcGGlsctlzXGRGnv6s+jFkRFn2z9oxy4bT69dvi9OuM8OWPex9jlol9+mb/sZcTlpNcL22plefpbx48cWtmwJa+MYl63hw1UzIw8OXvSNTZEsGjiMXk0jdmwdjSN2YUs4q1Rd50GFv6xoxhR7x9c/csaXRdYzqILX2zPmQT8ifUGSfe4xZmTdO+yW7GZcycWPTN6UU+m6ldrpQtInbq58P1mv77Hf78urU/5Nef8SfAjM9jzCjs3jGmJB+5+eEwuzXMLdymQ9jSN6aFeX3z9u2brjJR9CUdkfSspEdrZf9d0t9LeljSJySdm8pfIekHkh5KxwdrbV4r6RFJxyR9QGoGvhmzAiJGx9K37YRV+XzSUVf3UTtnTfNexiyJaSz9jwDXNsruB348In4C+Afg3bV7T0bEFen41Vr57cAvU2yWfqilT2NWTzMCZ9SK3FGTuc28O9Z40zEmin5EfBZ4rlH2VxGxnS4fBC4Z14eki4BzIuLBiAjgTuAtM43YmGUwLp2Cojpa71vsTYdZhE//PwB/Wbu+TNLfSvprST+Tyi4GjtfqHE9lrUi6RdJRSUe3OLWAIRozgZxW98xQWL7KrJqDQzVrPyXhHFDvzxkTTEeYK3pH0m8B28BHU9FJ4Ecj4luSXgv8uaTX7LbfiDgMHAY4R+fbZjJLYyh1QusmKsMWfrRZ+4FTMJjOMrPoS/pF4N8Cb0wuGyLiFBSmeUR8UdKTwKuBEwy7gC5JZcZ0DtWt8vq8LYzOstms1+zHmI4wk3tH0rXAfwZ+LiJeqJW/XFIvnb+SYsL2qxFxEnhe0jUpauftwCfnHr0xC6S+MKs1OKeWcK086uWDPqi5dsLWvukUEy19SXcBbwAukHQceA9FtM4B4P4UeflgitR5PfDbkrYovJi/GhHlJPCvUUQC/RDFHEB9HsCYlbHDIm/ucTs2985wOxi9x64tf9MFJop+RNzUUvzhEXXvAe4Zce8o8OO7Gp0xS2JkmoQdxc2onWaehl32b8yScRoGY2qMjdEv60xaVjhpj1xjVohF35iStHHKkGum6eoRw5Z+PVyzpd2gX2O6gUXfmBbK3bMUqkXupEVZ2UDEpRiEbUa93fLHbMw0WPSNqVNG8ARFWE5LGgY1rne2L/sIr8w1ncOib9abMqKmZUHVkMumLFOQZYMwHCmKW7UQz7Z+KvF3BI9ZMRZ9Y2o+91KglacQy6iFWgqUwUavJvoZA79+3tJ2aGLYZr9ZPRZ9s77kQC+dl66YvPDjq5/Eu18TcUBZzuZGv+pCyeov6wy1S/1V7p76+xqzIiz6xiQUASFUCn8p5DW9znrBZs29k/WGfyXU2ygJvmP0TZew6BuT3DCRxDqSYFcCnqcJXUGm4IyN7appVqbWrAQ/qqNy75RfHNZ+0wEs+sbUGLLy+5D1KUQ7GfdZlnOgVxP9mnuHvKhfuXhKa9+YDmHRN2uNooy+iXQu6AfKCgtfSchLK72XBQdrot8rY/YjCX4+sPjpR/UepU/frh6zaiz65vQlr8zzdB3j88rWM2vmoH4UVntflfD3ejlnbmxVTXq9fCD4/bJNDOYCpnHr1H8N5J7lNXuLRd+sPZWFn87Jk38/h2y7iOQpI256Wc7BXk30y0nd0h20naz80rc/FA5qK9+sHou+MRQCH1lAXiy5Ld06eT/I+lFN5vYUnLUx2MKzp6gmcbNk5ZdunmqbxAinVTadwaJvTKLaRCWJdfQh66eY/T4QopflnNm09GNQJ+tTzAkkK995eEzXsOgbA1WGzVLwRZABUVruabFVL8s5u/di1ayX5UWUTr/065e/DKKx2bqV33QDi74xJTFw1VS5draDbCvIknsnUwyJfpbcO1kO2Vag7WTl92srcS34pkNMtUeupCOSnpX0aK3sfEn3S3oivZ6XyiXpA5KOSXpY0pW1Njen+k9Iunnxj2NMC5NEd0funST+20GWjsrSV86Z2UvV0VPN0i/rbg8Ef1e5d/zlYJbAtBujfwS4tlF2G/BARBwCHkjXANdRbIh+CLgFuB2KLwmK/XVfB1wNvKf8ojCmC1T5daoVuUG2nVeiT4hMwcuyH1RHYemrJvr5YFVuDPdrTBeYSvQj4rPAc43iG4A70vkdwFtq5XdGwYPAuZIuAt4M3B8Rz0XEt4H72flFYsxyKa3ruu89Uqz9dl4cW0G2DeSwkeWc03uxOjayvFiJuw3aqrXpx7Brx7590xHm8elfGBEn0/k3gAvT+cXA07V6x1PZqPIdSLqF4lcCBzlzjiEaM5p6fH5VVuXXH+TPKSz9IuwyU3CmBiGbmdJCrP7AtaMIIq3oHfW+xqyKhUzkRkRIiwtMi4jDwGGAc3S+/4eY5RABUiX42s4RkG310qKrotpZ2amhZkqWfraVk22ltMsbGZFR9WdMV5jWp9/GM8ltQ3p9NpWfAC6t1bsklY0qN2Y8mufPdASjhLiM0+8H5Dna6pO91CfbLhZuZQQvy16qjoxUvh1kL/XRVr9IpVCP3tnN+8/DXvw7mdOOef5K7gXKCJybgU/Wyt+eoniuAb6T3ECfAt4k6bw0gfumVGbMasnZkS5BUVj6bPXRVl5Z+pmCg8qro3TvFD79VH87r/qo9+nNU0wXmMq9I+ku4A3ABZKOU0Th/C5wt6R3AF8H3paq3wdcDxwDXgB+CSAinpP0O8AXUr3fjojm5LAxi6F01eyqTa1tJEt/u7Des+1BKoWzGt0Woh9oq4+2+4Vbp9oYfZdjKN/fmD1iKtGPiJtG3HpjS90Abh3RzxHgyNSjM2ZVpAgetrbRS72hidwDNTdKfSJXL20X9XsZYeE2HcUrco2pUblhcqCfw9YW2ugN5dQ/qN5wozK18tY2bG3B5kbVhxTEbn9xGLOHWPSNaaHa+GQ7CXl/sNhqg16jLsXE7dZWUb/akMWY7mHRN6aNJPqxtYU2NqqYfYBeM0qm3BZxu1/Ud74d02Es+sY0qKJuovDps7k9tLF51gx6K3fb2i58+vWtEZuLv4xZNQ7sNWYUeS2VwiTDvZ5R05uhmw5j0TdmEip20xpfh92HiBqzAuzeMWYUmaDXg16vCLdPmp43V1kpheOnumQWf9NdLPrGNAip0HcJnbEJGz0iG1j7/WgR/Uyw0SvqJ4vfoZqmi1j0jWlDKo6NDdjchJ6qxbXb9IeqFla+inpb24O2xnQQi74xNQorv1hQpV5WLLTa3EgZM4s6L8aw6BeWPlVdelll5dvaN13Dom9MM5whK47IMrSxQWz2yDdUCDtwquHeiYzi/mavqJ9lVR87+nXSNbNiLPpmPWlY4FGP0CndM8l6j80eea8Q/TzEi7WFV3mk8l4S/c2Nol3ZR+oz0PDmKRKT40CNWTwWfbO2jHK9hISkFLVTiHlsUFn6L8awCR8Zxf3NHrHZgywDaXz/C30SY6bHom9MnbTTVfQEuUA98jMK9w5ZkIf4fmxW1fNI5RsiP6OHXuwVrqGehvozpitY9M3pybQTqPV65Xm50CrLiJ7INzLyXhGlkyNeyM+omuQUUT15D/KNjGyzV7hxmu6iUvincetUY/KXhVk8Fn1jYOB7TxuaR0ba51bkmxqy9J/PD1bNCku/mMjNN0VsZpBHiutX6m/QvzGrZuY0DJJ+TNJDteN5Se+S9F5JJ2rl19favFvSMUmPS3rzYh7BmDkpwyvTBGz12suIjYzYzIge1UTuC/mB6igncqNHUW8jI3rZcD/Z8PsYs0pmtvQj4nHgCgBJPYpNzj9BsT3iH0bE79frS7ocuBF4DfAjwKclvTqiGfRszILZpdjWV9/mG1kRjtmjsvS/u8PSj+QGKlxBWYrLDFHl4N+LcRozC4tKuPZG4MmI+PqYOjcAH4uIUxHxNYo9dK9e0PsbMx8pxDKSPz+ydGyo8Ncnn34/z/h+fqA6+nk28On3IDYGbev9WdBNV1iU6N8I3FW7fqekhyUdkXReKrsYeLpW53gq24GkWyQdlXR0i1MLGqIxU1BOwIok+MlXnzbLKidyyyNPPwmiR+H77xVfFGUfFnvTNeYWfUlnAD8H/Gkquh14FYXr5yTw/t32GRGHI+KqiLhqkwPzDtGYkUSKp4+MKpNmSERvYOlHWpgVWWHpf2/7YHX086y6F0nwI0tt0hdIpDQN5XsZs0oWEb1zHfCliHgGoHwFkPQh4C/S5Qng0lq7S1KZMauhxeSJtBI3smS19wr3Tmnp9/NsKGSznxedROne6RXtilQ9QeRFLp8d7+uZLLMiFuHeuYmaa0fSRbV7bwUeTef3AjdKOiDpMuAQ8PkFvL8xi6GWJC0yQU3wi0icoJ+L728fqI5+XszWlhE8eY8iI2dWs+pt3ZsOMZelL+ks4GeBX6kV/56kKyhWljxV3ouIxyTdDXwZ2AZudeSOWSr1zU2GFmUx8OWXrhoNXDZl3h0E/cj4QX+wIrcfWZVls/xlEFlUfVV766oRylN//0xOxGaWxlyiHxHfB364UfYLY+q/D3jfPO9pzEJp7nFeRtoksc83kmsn1evn4sW66Oeq+oke5BuF+PdyIFLYZt29Y9eOWTFekWvWjrbJ1BA7Jl7LHbHKiVqAPM94Ke9V7fLSp1/+QsgG7h7S90ckK7/5rk68ZlaBRd+cvmQtU1Zt/nXVrPUqgofKrVNF9lBa+oP/NqWlX7mDkptH/UiCDyqt+3HvXR9z3z8FzN5h0TfrS2uytUG4JjUhp2bpn9oe/LcpLf1qLqD+C6GX8u4Qw4nWPLFrVohF35hENCZ3qy0Sa8UBbNfy6Q8FY6rRrtbvjrBNY1aERd+YGlEXeQ27ewAiRMRA0cvzoXqN9rbrTZew6Jv1Rhq26KuJ3Fr+nbrlHrDVH0zkVgZ85QqqzwMUfv36LwDZtWNWjEXfmF1QWPrD18bsJyz6xjRJ2TEri712KyItyKpdV+cwHO3jzc9NB1lUlk1jTg+yFst9WmO+rV5bf8asEFv6xiSqydq2SVkoVtiGBqtwSe6dujHfNvlr3TcdwqJv1pPGhGpzle5gi0OGBJxG9A6N6J36XrvN/ofewRO6ZkVY9I2p0Qy5rMrK+zFIp1xeD7WttXPIpukiFn1jWtMjlPl4GuUxLPTNedpQi1to3PsYs2Q8kWtMyShRHiXi09ax2JsOYUvfrC9t2xeqYa03bo9akVtvv6OPenvJC7TMSrHoGzOGIRdPitTJh5z8DE3mOlLHdB27d4yB4f8JaSOVWQW8Wpg1tDvWPIMzZnHM/aco6SlJj0h6SNLRVHa+pPslPZFez0vlkvQBScckPSzpynnf35hlU7p4nILB7EcWZX/8q4i4IiKuSte3AQ9ExCHggXQNcB3FhuiHgFuA2xf0/sbMzQ73jBoHtcVY6aiEv6Vua5/GrJi9+tF5A3BHOr8DeEut/M4oeBA4V9JFezQGY2aiLtTtIZs1S78lZLPZhzFdYhGiH8BfSfqipFtS2YURcTKdfwO4MJ1fDDxda3s8lQ0h6RZJRyUd3eLUAoZozGhGiXMV2bMjgmdwDFEJfnuH/hIwXWAR0Ts/HREnJP1T4H5Jf1+/GREhaVepBiPiMHAY4Byd7zSFphMoYOf25mJ3f93GrJa5Lf2IOJFenwU+AVwNPFO6bdLrs6n6CeDSWvNLUpkx3aBp3bdY52MncnekYrB5b7rFXKIv6SxJLyvPgTcBjwL3AjenajcDn0zn9wJvT1E81wDfqbmBjFkZIbUIdnmvXpGBw77h02/LvVO+jnL5GLNs5nXvXAh8Iq0w3AD+d0T8H0lfAO6W9A7g68DbUv37gOuBY8ALwC/N+f7GGGN2wVyiHxFfBf5FS/m3gDe2lAdw6zzvacwy2GG1VzfGJ1xr/XVgTIdwGgaz3uzIjTPIhR9D7hkGAt8U+hjOs1OFbWYt4u8vA7NivDjcrB/z/NU3VXwek97/+8wKsKVvTBtj3DTj0i+MdAsZ0xEs+saUNDcxHyXck+Lym+2yum/ImNXiH5jG1JmwYUpzFe6OVbm76MuYVWBL35h6jH5iaBOVRpz+8Mbo9X4GR+sErmP1TQewpW+MMWuELX1jdovd82YfY0vfrC3TZMPcrb4PrdualL3TmBVgS9+YNkblw2/udt5SqW1DdGO6gi19YxKjNj6Zdv3VqLZOyWC6hC19s5ZM7WJphu6XWyXWr8fUH/f+/i4wq8CWvjFNpFYXjeq5d2ri3yb84RBN01Fs6RvDLlwwZT79+vUi+zdmj7Glb0ydRVvntvZNx7Clb0yJhs+HVuWWhIha7gXVTXg12rX1a8yKsaVv1ptRKRgWgFMxmC4ys+hLulTSZyR9WdJjkn49lb9X0glJD6Xj+lqbd0s6JulxSW9exAMYszTKPXF3HBZys3+Yx72zDfxGRHwpbY7+RUn3p3t/GBG/X68s6XLgRuA1wI8An5b06ojozzEGY/aEaPkFMDXVRujO12C6x8yWfkScjIgvpfPvAl8BLh7T5AbgYxFxKiK+RrE5+tWzvr8xe8FUG2OVK65Grcad1N6YFbIQn76kVwA/CXwuFb1T0sOSjkg6L5VdDDxda3acEV8Skm6RdFTS0S1OLWKIxkzPOOFu2R+3tV5LP8Z0gblFX9LZwD3AuyLieeB24FXAFcBJ4P277TMiDkfEVRFx1SYH5h2iMZMpF2RNUzUGxzR4oZbpEnOJvqRNCsH/aET8GUBEPBMR/YjIgQ8xcOGcAC6tNb8klRnTPdq8N6NEvrmblhOumQ4zT/SOgA8DX4mIP6iVX1Sr9lbg0XR+L3CjpAOSLgMOAZ+f9f2N2RWlpd3cB3diu8Z1W/TOuPqTKMfjXwJmScwTvfNTwC8Aj0h6KJX9JnCTpCso/js8BfwKQEQ8Julu4MsUkT+3OnLHdI4FZsW0xW+6yMyiHxF/Q/uf9H1j2rwPeN+s72nMUhkZudO4nqadMR3BK3KNaWFea9+hmqarOPeOWXuq3PotfvUhF000Xmlcj3INpYVa3ibRdAFb+mZ9GfXXP0abm2GaY8M2R/Xj/3VmhdjSN4aBhR6lII+bhJ0Un19rW/Znd4/pChZ9s9bEbq3uZpjmLtPrROZ5XrNa/EPTmAZDVnk9Xf4IgR8qr2+qZXU3HcSWvjGJptiPFO0Jat6Mz7f4my5hS9+YOosWaAu+6Ri29I2ZRNPFE43rlnrGdBVb+mZ9mDZOftT+ttMkXJtlb1zH75slYkvfmIz2HPpzhGy25ta3iWU6gP8MjTFmjbDom/Vk1r/8Re2H6P95ZkX4T8+YkgWmVS5xemXTNSz6xhizRlj0jWlQ5eFh95Z/Pd2+F2WZLuLoHbN+TBMi2bqBSnGomXunLZpnGsF3qKZZAUu39CVdK+lxScck3bbs9zfGmHVmqaIvqQf8EXAdcDnFfrqXL3MMxhizzizb0r8aOBYRX42Il4CPATcseQzGGLO2LNunfzHwdO36OPC6ZiVJtwC3pMtTn46PP7qEsS2LC4BvrnoQC2S+5xm1ujWfucdF4M+o2/h5JvPPRt3o5ERuRBwGDgNIOhoRV614SAvDz9N9Trdn8vN0m2U/z7LdOyeAS2vXl6QyY4wxS2DZov8F4JCkyySdAdwI3LvkMRhjzNqyVPdORGxLeifwKaAHHImIxyY0O7z3I1sqfp7uc7o9k5+n2yz1eRSxy52djTHG7FuchsEYY9YIi74xxqwRnRX90yFdg6SnJD0i6SFJR1PZ+ZLul/REej1v1eMch6Qjkp6V9GitrPUZVPCB9Jk9LOnK1Y28nRHP815JJ9Ln9JCk62v33p2e53FJb17NqEcj6VJJn5H0ZUmPSfr1VL4vP6Mxz7OfP6ODkj4v6e/SM/3XVH6ZpM+lsf9JCm5B0oF0fSzdf8VCBxQRnTsoJnmfBF4JnAH8HXD5qsc1w3M8BVzQKPs94LZ0fhvw31Y9zgnP8HrgSuDRSc8AXA/8JUW6sWuAz616/FM+z3uB/9RS9/L0t3cAuCz9TfZW/QyNMV4EXJnOXwb8Qxr3vvyMxjzPfv6MBJydzjeBz6V/+7uBG1P5B4H/mM5/DfhgOr8R+JNFjqerlv7pnK7hBuCOdH4H8JbVDWUyEfFZ4LlG8ahnuAG4MwoeBM6VdNFSBjolI55nFDcAH4uIUxHxNeAYxd9mZ4iIkxHxpXT+XeArFCvf9+VnNOZ5RrEfPqOIiO+ly810BPCvgY+n8uZnVH52HwfeKC0uJWtXRb8tXcO4D76rBPBXkr6YUksAXBgRJ9P5N4ALVzO0uRj1DPv5c3tncnccqbnc9tXzJDfAT1JYkvv+M2o8D+zjz0hST9JDwLPA/RS/SP4xIrZTlfq4q2dK978D/PCixtJV0T9d+OmIuJIiq+itkl5fvxnF77d9HTN7OjwDcDvwKuAK4CTw/pWOZgYknQ3cA7wrIp6v39uPn1HL8+zrzygi+hFxBUUWgquBf76qsXRV9E+LdA0RcSK9Pgt8guLDfqb8OZ1en13dCGdm1DPsy88tIp5J/ylz4EMM3AP74nkkbVII5Ecj4s9S8b79jNqeZ79/RiUR8Y/AZ4B/SeFaKxfI1sddPVO6/0+Aby1qDF0V/X2frkHSWZJeVp4DbwIepXiOm1O1m4FPrmaEczHqGe4F3p4iRK4BvlNzMXSWhk/7rRSfExTPc2OKprgMOAR8ftnjG0fy9X4Y+EpE/EHt1r78jEY9zz7/jF4u6dx0/kPAz1LMVXwG+PlUrfkZlZ/dzwP/N/1aWwyrntkeM+N9PcXM/ZPAb616PDOM/5UUUQV/BzxWPgOFb+4B4Ang08D5qx7rhOe4i+Ln9BaF3/Edo56BIkrhj9Jn9ghw1arHP+Xz/HEa78PpP9xFtfq/lZ7nceC6VY+/5Xl+msJ18zDwUDqu36+f0Zjn2c+f0U8Af5vG/ijwX1L5Kym+oI4BfwocSOUH0/WxdP+VixyP0zAYY8wa0VX3jjHGmD3Aom+MMWuERd8YY9YIi74xxqwRFn1jjFkjLPrGGLNGWPSNMWaN+P9QS/896i8bXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for j, U in enumerate(UU[::skip]):\n",
    "    V = U[::2]\n",
    "    W = U[1::2]\n",
    "    Umag[j,:] = np.sqrt(V**2+W**2)\n",
    "    \n",
    "plt.pcolor(Umag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "812fd7e7",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}