Método predicción-correción .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Método predictor-correcctor: Runge-Kutta 4 + Adams-Bashforth-Moulton 4 pasos\n",
    "## Para resolver la ecuación $y'=f(t,y)$ sujeta a $y(t_0)=y_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un método *predictor-corrector* (o de *predicción corrección*) es la combinación de una técnica explícita con una implícita. Aquí aplicaremos el **método RK4**:<br/><br/>\n",
    "    \\begin{gather*}\n",
    "        y_{n+1}=y_n +\\frac{h}{6}(k_1+2k_2+2k_3+k_4) \\\\\n",
    "        \\text{donde:}\\\\\n",
    "        k_1=f(x_n,y_n),\\\\\n",
    "        k_2=f(x_n+\\frac{1}{2}h,y_n+\\frac{1}{2}hk_1),\\\\\n",
    "        k_3=f(x_n+\\frac{1}{2}h,y_n+\\frac{1}{2}hk_2),\\\\\n",
    "        k_4=f(x_n+h,y_n+hk_3);\\\\\n",
    "    \\end{gather*}<br/><br/>\n",
    "para inicializar el problema y obtener $\\{y_i\\}_{i=0}^3$. Después, usaremos la fórmula explícita de Adams-Bashforth (*parte predictora*) de 4 pasos:<br/><br/>\n",
    "    \\begin{gather}\n",
    "        y_{i+1}=y_i+\\frac{h}{24}\\left[ 55 f (t_i,y_i) -59 f (t_{i-1},y_{i-1})\n",
    "            +37f(t_{i-2},y_{i-2}) -9f(t_{i-3},y_{i-3})\\right]\n",
    "    \\end{gather}<br/><br/>\n",
    "combinada con la fórmula implícita de Adams-Moulton (*parte correctora*) de 4 pasos:<br/><br/>\n",
    "    \\begin{gather}\n",
    "        y_{i+1}=y_i+\\frac{h}{720}\\left[ 251 f(t_{i+1},y_{i+1})+646f(t_i,y_i)-\n",
    "        264f(t_{i-1},y_{i-1})+106f(t_{i-2},y_{i-2})-19f(t_{i-3},y_{i-3})\\right]\n",
    "    \\end{gather}<br/><br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Definir el método numérico"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definir el método RK4 (inicialización)\n",
    "#f0=función, (x0,y0)=condición inicial, xf=punto final, h0=tamaño de incremento\n",
    "def rk4(f0,x0,y0,xf,h0):\n",
    "    n0=int((xf-x0)/h0)\n",
    "    X=[x0+i*h0 for i in range(n0+1)]\n",
    "    Y=[y0]\n",
    "    for x in X:\n",
    "        k1=f0(x,Y[-1])\n",
    "        k2=f0(x+h0/2,Y[-1]+h0*k1/2)\n",
    "        k3=f0(x+h0/2,Y[-1]+h0*k2/2)\n",
    "        k4=f0(x+h0,Y[-1]+h0*k3)\n",
    "        Y.append(Y[-1]+(h0/6)*(k1+2*k2+2*k3+k4))\n",
    "    return X,Y[:-1]\n",
    "\n",
    "# Definir la predicción (xj=(ti,yi) son las duplas (t_{j-i},y_{j-i}))\n",
    "def pred(f0,x0,x1,x2,x3,h0):\n",
    "    return x0[1]+(h0/24)*(55*f0(x0[0],x0[1])-59*f0(x1[0],x1[1])\n",
    "                          +37*f0(x2[0],x2[1])-9*f0(x3[0],x3[1]))\n",
    "\n",
    "# Definir la corrección\n",
    "def corr(f0,x0,x1,x2,x3,x4,h0):\n",
    "    return x1[1]+(h0/720)*(251*f0(x0[0],x0[1])\n",
    "                           +646*f0(x1[0],x1[1])-264*f0(x2[0],x2[1])\n",
    "                           +106*f0(x3[0],x3[1])-19*f0(x4[0],x4[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Combinar las rutinas\n",
    "def pred_corr(f0,t0,y0,tf,h0):\n",
    "    n0=int((tf-t0)/h0)\n",
    "    T=[t0+i*h0 for i in range(n0+1)]\n",
    "    Y=rk4(f0,t0,y0,t0+3*h0,h0)[1]\n",
    "    for i in range(3,n0):\n",
    "        y=pred(f0,(T[i],Y[-1]),(T[i-1],Y[-2]),(T[i-2],Y[-3]),\n",
    "               (T[i-3],Y[-4]),h0)\n",
    "        y=corr(f0,(T[i+1],y),(T[i],Y[-1]),(T[i-1],Y[-2]),\n",
    "               (T[i-2],Y[-3]),(T[i-3],Y[-4]),h0)\n",
    "        Y.append(y)\n",
    "    return T,Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Definir la ecuación diferencial con las condiciones iniciales, intervalo de integración y tamaño de paso\n",
    "$y'(t)=ty^3\\sin(-ty),\\;y'(1)=-2,\\;t\\in[1,10]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "f=lambda t,y:t*y**3*math.sin(-y*t)\n",
    "g=lambda t,y:t*y**3*math.cos(-y*t)\n",
    "t_inicial,y_inicial,t_final=1,-2,10\n",
    "h=0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7faf54a63128>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAF1CAYAAAD4CWwHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8VPW9//HXN5NJJvtOAoQlsu+oARUVqFvVWndrrVat7VXba1u7WLvdXm9/tdLae7taXOvaxWql7iuKgiuLgCAgW4AQCAnZl0kyM9/fH2eAgFnJcmaS9/PxmMfMmTlzzmdGDG+++Zzv11hrERERERGRjsW4XYCIiIiISDRQcBYRERER6QIFZxERERGRLlBwFhERERHpAgVnEREREZEuUHAWEREREekCBWcRGRCMMdcaY5a5XUdfMcYUGWPO6OD1rxtjSo0xdcaYrF46pzXGjO2NY/WWzr4HEZG+pOAsIhHPGLPEGFNpjIl3u5bOhMNmfTjAlhtj/m6MSe/jc3qB/wPOstYmW2v3H8UxlhhjvtbDOmw4vMe2ei7WGLPPGNPriwYYY24zxjzW28cVEWmPgrOIRDRjzGjgVMAC57taTNfNsNYmA8cAGcBtfXy+XMAHrO/uG42jN/8uqALOabV9LlDZi8cXEXGNgrOIRLqrgfeAh4BrDjxpjMkyxjxjjKkxxnwAjGn9JmPM740xu8KvrzTGnNrqtduMMU8YYx4zxtQaYz4yxow3xvwoPDq6yxhzVqv9rzXGbAvvu90Yc2VXCrfW1gDPAJNbHesrxpgN4WNtM8bc0Oq1bGPMc8aYKmNMhTFm6RGhdqYxZq0xptoY87gxxmeMGQ9sCr9eZYx5PXysOcaY5eF9lxtj5rQ6zxJjzO3GmLeBBuBRnH+c/Ck8Uv6nVuc8wxizOTzif5cxxnTysR/F+W92wNXAI613MMYMC/+3qzDGbDHG/Eer1x4yxvyi1fZ8Y0zxkScxxpwN/Bi4PFzzmi4c+zZjzD+NMY+Ev//1xpjCTj6PiMhBCs4iEumuBv4avn3WGJMbfv4uwA8MBa4L31pbDswEMoG/AU8YY3ytXv88TsjLAD4EXsb5mTgc+DlwD4AxJgn4A3COtTYFmAOs7krhxpgM4EKc4H/APuA8IBX4CvBbY8xx4de+BxQDOTijyD/GGWk/4AvA2UABMB241lr7CTAl/Hq6tfY0Y0wm8Hy47iycNo7nj+h9/jJwPZACXAssBW4Kt3rc1Gq/84BZwIzw+T/bycf+NzDXGJMeblE5FXj6iH3+Hv6cw4BLgV8aY07v5LiHsda+BPwSeDxc84wuHvt84B9AOs4/alr/I0FEpEMKziISsYwxpwCjgH9aa1cCW4EvGWM8wCXAz6y19dbadcDDrd9rrX3MWrvfWhuw1v4vEA9MaLXLUmvty9baAPAETlhdYK1twQlWo1v1JoeAqcaYBGvtHmttZy0Rq4wxVUA5MJJwCA/X9by1dqt1vAm8ghMuAVpw/iEwylrbYq1daq1tHZz/YK0tsdZWAM/i/MOgLZ8DNltrHw1//r8DG3H+sXDAQ9ba9eHXWzr4LAustVXW2p3AGx2c8wB/uLbLgS/ihFP/gReNMSOAU4BbrbV+a+1q4H6cIN8jXTz2MmvtC9baIM4/nGa0cSgRkTYpOItIJLsGeMVaWx7e/lv4uRwgFtjVat8drd9ojPleuCWiOhxi04DsVruUtnrcCJSHw9SBbYBka209Tgi8EdhjjHneGDMxfI714TaButatIMBx1tp0nL7jhcDSA6PdxphzjDHvhVsJqnB6gA/UdSewBXgl3MbxwyO+j72tHjcAyW18Z+CMtu444rkdOKPpB+yia9o8ZwefHZzWjKtpo00jXFuFtba2g9qOVleOfeTn8ZlWFzOKiHREwVlEIpIxJgGnNWCeMWavMWYv8B2cEcJcIACMaPWWka3eeypwa/j9GeEQWw101p/bpvDI9Jk4o8EbgfvCz08JtwkkW2uXtvG+FpwRzwKcEet44F/Ab4DccF0vHKjLWltrrf2etfYYnNHh73a3hSGsBGekvrWRwO7W5R1ZbndO0MlnX4rzXeUCR04RWAJkGmNS2qmtHkhs9VpeR2V089giIj2i4CwikepCIIhzYd3M8G0STii7GngKuM0Yk2iMmUyrCwdx+nYDQBkQa4z5GU5PcbcZY3KNMeeHe52bgLpwXV15rwenj7kR2AbE4bSMlAEBY8w5QOuLEM8zxowNX4BXEz5Pl851hBeA8caYLxlnOrjLcb7H5zp4TynOLCA9Fm4v+Txw/hGtJlhrdwHvAHeEL26cDnwVp4cdnP7xc40xmcaYPODmTmoefeACyi4cW0SkRxScRSRSXQM8aK3daa3de+CGczHXlcBNOG0De3Fm3Hiw1XtfBl4EPsH5Vb2frrcmHCkG56K9EqACmAd8o5P3rDHG1OFMw3YNcJG19kALwbeAf4Zf+xJOD/AB44DXcML5u8CfrbVLultweB7n88J17wd+AJzXquWlLb8HLg3PnvGH7p6zjRrWd9ALfgUwGuc7XQT8t7X21fBrjwJrgCKc/u/HOzjNE+H7/caYVV04tohIj5gjBgNERERERKQNGnEWEREREekCBWcRERERkS5QcBYRERER6QIFZxERERGRLlBwFhERERHpgoheLSk7O9uOHj3a7TJEREREZABbuXJlubU2p7P9Ijo4jx49mhUrVrhdhoiIiIgMYMaYHV3ZT60aIiIiIiJdoOAsIiIiItIFCs4iIiIiIl0Q0T3ObWlpaaG4uBi/3+92KQOKz+cjPz8fr9frdikiIiIiESnqgnNxcTEpKSmMHj0aY4zb5QwI1lr2799PcXExBQUFbpcjIiIiEpGirlXD7/eTlZWl0NyLjDFkZWVpFF9ERESkA1EXnIFBH5off/xxioqKevWYg/07FREREelMVAZnt91+++1MmTKF6dOnM3PmTN5///0O958/f/5Rz0c9Z86cw7Yfe+wxdu7cSXcXhnniiSeYMmUKMTExmhtbRERE5ChEXY+z2959912ee+45Vq1aRXx8POXl5TQ3N/fZ+d55553Dtq+66qqjOs7UqVN56qmnuOGGG3qjLBEREZFBRyPO3bRnzx6ys7OJj48HIDs7m2HDhgGwePFijj32WKZNm8Z1111HU1PTp96fnJx88PGTTz7JtddeC0BpaSkXXXQRM2bMYMaMGQcD84H9rbXccsstTJ06lWnTpvH4448DsGTJEubPn8+ll17KxIkTufLKK7HWfuq8kyZNYsKECb33RYiIiIgMMgrO3XTWWWexa9cuxo8fzze+8Q3efPNNwLlo8dprr+Xxxx/no48+IhAIsHDhwi4f91vf+hbz5s1jzZo1rFq1iilTphz2+lNPPcXq1atZs2YNr732Grfccgt79uwB4MMPP+R3v/sdH3/8Mdu2bePtt9/uvQ8sIiIiIsAAaNX40VNrKa359Mju0cpNjeeOi6e3+3pycjIrV65k6dKlvPHGG1x++eUsWLCAY489loKCAsaPHw/ANddcw1133cXNN9/cpfO+/vrrPPLIIwB4PB7S0tIOe33ZsmVcccUVeDwecnNzmTdvHsuXLyc1NZXZs2eTn58PwMyZMykqKuKUU045mo8vIiIiIu2I+uDcUcjtKx6Ph/nz5zN//nymTZvGww8/zMyZM7v03tazV3Rn+re22i8OONA2cqC2QCDQ5eOKiIiIdCgUhIAfWvzO/YHbwe3GDl5rY7ulEYIth47femavi+6GhIz+/4xdFPXBub9t2rSJmJgYxo0bB8Dq1asZNWoUEydOpKioiC1btjB27FgeffRR5s2b96n35+bmsmHDBiZMmMCiRYtISUkB4PTTT2fhwoXcfPPNBINB6uvrSU1NPfi+uXPncs8993DNNddQUVHBW2+9xZ133snGjRv754OLiIhIZLD2UABtaXBCaUuDsx1obPv5tl4LNEIgPMFBR9PSmhiIjYfYBOfemwCxPufm9R16HOuD+BRIygm/ltDG+8LbHm/H54xQCs7dVFdXxze/+U2qqqqIjY1l7Nix3Hvvvfh8Ph588EEuu+wyAoEAs2bN4sYbb/zU+xcsWMB5553HiBEjmDp1KnV1dQD8/ve/5/rrr+eBBx7A4/GwcOFCTjrppIPvu+iii3j33XeZMWMGxhh+/etfk5eX1+XgvGjRIr75zW9SVlbG5z73OWbOnMnLL7/cO1+KiIiIOCOznwqsDeGQe2C7i68R/k3zgd84G3P44wPB1JsA3sRW24mHno9PhuQh4dcOPO87fH9PXFQGWLeYjloA3FZYWGiPnHN4w4YNTJo0yaWKBjZ9tyIiMuAFA9BSD80NTkhtrj/ivqGd18PPtzQeCrCtGRMemfVBXGKrcNoqyLYOtp29FuPp/+9mEDPGrLTWFna2n0acRUREJPIEA9Bc5wTX5jpoqgtvt37cTugNtDNpgLVOIPUmhsNtknMfl3TosTcRUoaG90n69L6xCRCjSckGKwVnERER6blAUzvhNhxwm+qgubbV43pnu+WIC+UPtCTEeCAu2Qmv8ckQl9LqcTKkDm8VaJOOCLc+tR9In1BwFhERGaxCQWiqhaaa8H0t+GsO3z7ycaD58H7bA2Ljjgi3Sc72gcdJQ1o9Hw6/8ckKuRJVFJxFRESizcHAW3tEwK1pFX5bv1Ydbl8wh4feGI8zC8Jht1TnlpwL2eOOeD7FmRVBZJBScBYREelvoZATcv3V4VtVq8dt3ELBw99vYg4Pu77UQ8E2KQcyx7R6LU2BV6SXKDhHoccff5wTTjiB0aNHu12KiMjgZK0zu0JjpXPrLPi2XuwBDgXfhHQn2La+peW32k53ArFHf12LRAL9n3gUbr/9dv72t7/h8XiIiYnhnnvu4YQTTmh3//nz5/Ob3/yGwsJOZzn5lDlz5vDOO+8c3H7sscfYs2cPl19+ebeOc8stt/Dss88SFxfHmDFjePDBB0lPT+92PSIiA4q1zkVqBwJwY0Wrx61vVRAKHN7X601wVjhLSHfuD4TdlLxDodeX7owGe7zufUYR6TUKzt307rvv8txzz7Fq1Sri4+MpLy+nubm5z87XOjQDXHXVVUd1nDPPPJM77riD2NhYbr31Vu644w5+9atf9UaJIiKRIdAEDfsPv7UOvg0VzsiwDR3+vrikcABudUsZCkMmh7cznSAcG+fO5xKRiKHg3E179uwhOzub+HinVyw7O/vga4sXL+b73//+wZUDFy5ceHC/A5KTkw+uFvjkk0/y3HPP8dBDD1FaWsqNN97Itm3bAFi4cCFz5sw5uL+1lh/84Ae8+OKLGGP46U9/yuWXX86SJUu47bbbyM7OZt26dRx//PE89thjmCOuUD7rrLMOPj7xxBN58skn++T7ERHpFaGQE3Jbh+D68iNCcRUHV1cDZ1Q3MQsSs8P3mZA2EvJmHArEvjS1PYjIUdNPj24666yz+PnPf8748eM544wzuPzyy5k3bx5+v59rr72WxYsXM378eK6++moWLlzIzTff3KXjfutb32LevHksWrSIYDB4MFwf8NRTT7F69WrWrFlDeXk5s2bNYu7cuQB8+OGHrF+/nmHDhnHyySfz9ttvc8opp7R7rr/85S/dbvUQEemRUNAZ8a3fB3X7nBBcX9YqCJcf0QdsnBaIxKzDbzkTDj32pWshChHpV9EfnJ/5FtTu7b3jpeTB+X9o9+Xk5GRWrlzJ0qVLeeONN7j88stZsGABxx57LAUFBYwfPx6Aa665hrvuuqvLwfn111/nkUceAcDj8ZCWlnbY68uWLeOKK67A4/GQm5vLvHnzWL58OampqcyePZv8/HwAZs6cSVFRUbvB+fbbbyc2NpYrr7yyS3WJiLQr0OyE3/pwEK7bFw7GZc7zTTXOftY6F8MlZkFyjjOfb1IODJkESeHR4YRM8Prc/TwiIp2I/uDcQcjtKx6Ph/nz5zN//nymTZvGww8/zMyZM7v03tYtFH6/v4M9D2ePnGi+ldbtIB6Ph0Ag0OZ+Dz/8MM899xyLFy/+VCuHiAjgLHNcXwZ1e6G21Lmv2xd+bp8zk8QBHq8TgJNyIDkchjNPCAfjbKctQj9rRKQToZClvjlArT9AbqoPT0zk/tyI/uDczzZt2kRMTAzjxo0DYPXq1YwaNYqJEydSVFTEli1bGDt2LI8++ijz5s371Ptzc3PZsGEDEyZMYNGiRaSkpABw+umnH2ztCAaD1NfXk5qaevB9c+fO5Z577uGaa66hoqKCt956izvvvJONGzd2qe6XXnqJX/3qV7z55pskJib2wjchIlEl2OIE39q94VC899DjujKwQcA4C2IkZUNynvMbuJQ8GDLl0EhxnH5+iMghLcEQtf4Adf4ANf4Wav0Bav0t1DUFDj6uPfg4QENTgCOHAg2QGB9Lii+WH50zkRRf5M5Co+DcTXV1dXzzm9+kqqqK2NhYxo4dy7333ovP5+PBBx/ksssuO3hx4I033vip9y9YsIDzzjuPESNGMHXq1IO9zL///e+5/vrreeCBB/B4PCxcuJCTTjrp4Psuuugi3n33XWbMmIExhl//+tfk5eV1OTjfdNNNNDU1ceaZZwLOBYJ33313L3wjIuKqUMjpD67ZDTUlzu1AKG4oPzSDRIzXCb8pQ50V4VKHwfDjnICclKML5kQGoWDIUutvoabRCb3VjS3UNLYcDMA14UBcG95uDoY+dYzYGEOyL5ZUn5cUXyzJ8bGkhB+PykokxXdoOzk+lqS4WGIieES5M6ajFgC3FRYW2hUrVhz23IYNG5g0aZJLFQ1s+m5FIkwo5PQMtw7F1cXOvb863AZhnOCbOuzQLSUvHIiznRFkERmQDrQ4OIHXCb9O8A0cFoIPvBYMHZ75Yowh1RdLaoLXuR143Crstg69Pu/A/XlijFlpre10wQ0NMYiIuMFa54K66p1QtSscjMMBubkeCF9Ql5QDqcOdQJw1BgrmOo/jU9U/LDIAtARDVDe2UNXQQnVjM1UNrcJuOPgeCMZHjvgaICk+lrQEL6kJzqhvaoKXtAQveam+w55L9cUS69EsND2l4Cwi0hdCIad/uCocjKt3Oo9rSw+1TyRmQfpIZ4nlIRNh7OnhUJzibu0i0m3+liA1jS1UNrRQ1dBMVWML1Q0tVIXD8IHthubAYRfoe2IM6Qle0hO9pCfGHQy+o7ISw4HX2U7xDewR32ih4CwicjSCAWeEuCociKt3OQG5vsx53RinZSJthBOOR86B6V90nlP7hEjE8rcEqag/EHabw+G35fDthpZPjf7GeWLISPKSlhDnhOBwGB6RmRgOxV7SE+LweWM0s1UUi8rgbK3VH7peFsm97iKuaaqDyiKo3O7cV2x3wrINQUysMzqcPgrSR8C4M51V6pKy1UIhEiGaAkGqGlqoqG+msqGZyvoWKhqaqaxvPvhcrf/wKVx93hjSE+PICAfdtEQvWUlxjMlJPhiIUxO8Gv0dpKIuOPt8Pvbv309WVpbCcy+x1rJ//358Pi0+IIOMtc7sE62DcWWRs9QzQFwyZIx2bkNnwOQLnH5jjRiL9LtAMERlQ0s4ADuht6K+JXzvPFfd2HLYVGdejyEjMY6MpDgyw/fD031MHZZKZpKznRIfqzwhXRZ1wTk/P5/i4mLKysrcLmVA8fl8B1cfFBlwGipg/1bYvxn2b3EetzQABlJyIaPACccTz3XuEzJcLlhk4AuFLNWNLeyvb6K8rpn9dc2tHjexv675YDuEtZZYTwzpCd6DgTczMY6clHgm5CWTkRhHZlIcqT5vVE91JpEv6oKz1+uloKDA7TJEJNI0N0DFtnAwDofj+jKnbSIhA7LGOrfJF0LmMRCf7HbFIgNOQ3OA/XXNlIeD75GhuHVbhAHSEr1kJ8eTlRRHVnI8IzITOXZEBlnJThBWO4REmqgLziIyyDVUQNnG8G2T01oRCoI3wZmuLWssjD4Fjr/WmbVCv4IV6ZGG5gD7apooq2uirNa57av1U1bbREV982FzAyfGx5IdDsFZyXFkJcUzJif54LbaIiTaKTiLSOSx1lke+kA4LtvozFoBkJAJOeMhZyIUzHMuztOqdyLdEgiGqKhvZl/toTDcOhg3NB8aGU6MjyUnOZ6cFOc2KiuRwtEZ5CTHk5kUp7mBZVDR3zYi4q76cihdB6UfOwG5psR5PjkXciY4AXn8Z51p3WL0F7RIR/wtQfbVNLG3xk9pjf+wYFzZ0Aw4/cKeGENWkhOEh6TGk5MczzE5SQfDcWKc4oFIW/R/hoj0j0AzlH8CpeudoFy+GUIBp50idwrkTg7PWjFM7RUiRwiFLOX1TZRWN1Fa4z8YjEtr/OyvayYUnlI0PtZDbmo8Q1J95Kb6GJ+bzMljs8hJjicjMU4Xzon0kIKziPQua6Gu1AnHe9c5QbmxAjxxkD0ecqfCjCsgexx4vG5XK+K6Wn9LOAQ3sbfaT2mtn9JqZ7spEAQgxhiykuPIS/UxJNVHXqqPGfnp5KbFk5UUj0eBWKRfKDiLyNGz1uk9LlkNe9bAvo8h2AzJeZA3FYYfB8d+GZKy3K5UxBWNzUFKqhvZU+WnpLqRkirncXldExanbSLZ5yUvNZ7c8CjxrNGZ5Kb4GJIar1klRCKMgrOIdI21zgwWew6E5I1Oq0VaPgybCZPOg3k/gNh4tysV6RctwRB7q/1OGK72HwzIe2v8BMLzDyfEeRialsDQNB/D0hOYNz6H4ekJZCfHq21CJAr1KDgbYzKBx4HRQBHwBWttZRv7FQG1QBAIWGsLe3JeEekHNSVQvAJ2r3BmtrAhZwaLoTNg6iXORXtqtZABrMbfQnFFI7urGimubKC40hkx9rcEMcYQG2PITfUxNN3H8PQECkdlMnS6M2ocF6sLWUUGop6OOP8QWGytXWCM+WF4+9Z29v2Mtba8h+cTkb7QXO+0W+xeAbtXQVMtpA6F/Fkw7QvhkKxfUMnAUt3YcjAQO7cGSqoaaQ44o8UpPi/5GQkMz0jgmJxkTh3njBYnxKl9QmSw6unfhBcA88OPHwaW0H5wFpFIEAo5S08XL3dGlKt2OouHDDsW8gvh+K+AL9XtKkV6rLqxhV0VDQdD8YER45ZgCGMMKb5Y8jMSyM9IZHxuMqdNHMLQNJ/6ikWkXT0NzrnW2j0A1to9xpgh7exngVeMMRa4x1p7b3sHNMZcD1wPMHLkyB6WJyIEmp2+5J3vwq4PoKXRmd0ivxBO+Q6kj9T0bxKVQiHL3ho/O/Y3sLOiPnzfQH2Ts3hHaoKXERmJ5GckMDEvlTMm5TI03Ud8rIKxiBydToOzMeY1IK+Nl37SjfOcbK0tCQfrV40xG621b7W1YzhU3wtQWFho29pHRDrgr4HiD2DHu7B3LZgYGDoTRp4IhV+F+GS3KxTpMn9LkF0VTiA+EIyLKxsJhkLEGENumo+RmYmMykzkc9OHMjIzkRSfeu9FpG90GpyttWe095oxptQYMzQ82jwU2NfOMUrC9/uMMYuA2UCbwVlEuqmxEorehqKlzqIi8ckw4gSYeC7M/6Eu4JOI19gcpGh/PdvLnVtReT37651V7uI8MYzITGBkVhJjhyRz+qQhDEtPwKtlnkXEBT1t1XgGuAZYEL5/+sgdjDFJQIy1tjb8+Czg5z08r8jg1VTntF1sf9NZptqXCqNPcUaTs8ep7UIiUkswRHFlI0Xl9Wwrr2d7eR0lVX6stfi8HkZnJ1GQncQJBZl8oXAE2clxGP1ZFpEI09PgvAD4pzHmq8BO4DIAY8ww4H5r7blALrAo/AMwFvibtfalHp5XZPBo8TutF9vfcma+8Ppg5ByY/kUYMhliNPImkcFaS2lNE9vK65zR47J6dlQ0EAiG8MQ4I8cF2UmMz03ms1NyGZaWoLmMRSSqGGsjt424sLDQrlixwu0yRPqXtc68yVsXQ9EyZ3vEbCiY58yhrGnhxGUtwRA79jewZV8dW8vq2LKvjsqGZgyQl+ajIDuJguxkCrKTGJGZoIvxRCTiGWNWdmWdEf0NLBIJ/NWw7U3Y8pozPVzOBBh7hjM1XFyi29XJINXYHGRr2aFwvLWsjsbmIJ4Yw8jMJMblJjO7IJMvzR5JRlKc2+WKiPQ5BWcRN1jrzHix+RVnijhvgjOifOp3IWO029XJIFPfFOCT0lo2l9axeV8t28sbCIRC+GI9HJPjXJR35uRcbswZQ1K8/toQkcFLPwFF+kugCbYvhU0vQMVWyJsG48+GOd+GWI3WSd8LBENsL69n495aNu2tZfO+WpoCIRLjPIwbksL43BQuKxzB6KwkLRktItIGBWeRvlRf7owqf/IytDRAwVw48RuQPdbtymQAs9ayp9rPpr214ZBcQ1VjC7ExhtFZSUzIS+GzU/L4z8+M1fLRIiLdoOAs0tuqdsLHzzjTxcWnwLiz4LzfQmKm25XJAORvCbJpby3rS2r4eE81xZWNAAxNS2BiXgozR6TzxVkj1IMsItILFJxFekNlEXz8tHOBX0oeTL4AZv8HxMa7XZkMIFUNzXxcUhMOyTVUNjQTHxvDhNwUJg9L44a5Y8jPSND8xyIifUTBWeRoVWxzwvL2pZA6DCZfCCd8Xf3K0mPWWnZXNR4MyRv31uBvCZGe6GXy0FSmDEvjkuPzydQosohIv1JwFumO6t2w7knYtgTSRsCUC+Gkm7SstRy1AyH5o+Jq1hRXs7m0lkDIMjwjgclDU5k/IYcb541RL7KISARQcBbpjL/a6Vne+Bz40mDaZc4FfgrLchT21foPhuQNe2poDoQYlp7AjPw0Pj9jKONzx+P1aEYLEZFIpOAs0pZAM2x5Fdb9y1nyevL5cMn9zsV+Il1U1dDMR7urWVtczbrd1TQ0B8lOjmfGiDROmziEb8wfg8+rkWQRkWih4CzS2p618OGjUL7ZmQ3js3dASq7bVUkUCIYsm/bWsmpnJR/urKK8rom0BC/T89OYNTqTa+eM1uIhIiJRTj/FRRqrnL7lDc85S10fdw3kTXW7KolwFfXNfBgOyRv21BCylvF5KRw7IoNbz5nAkBSf2yWKiEgvU3CWwcla2PE2fPiY08M89RK44h/gVdiRT2s9mrxqZyXldc1kJHo5dkQ6Z03J5dtnjFNfsojIIKDgLIPk4ChRAAAgAElEQVRLYxWs+buz7PWIE2H+DyFjtNtVSYRpCgT5qLia97dXsHpXFS3BEBNyUzhuVAY/PGeiRpNFRAYpBWcZHErXwwf3Qe0emHEFXPWUZsWQg+qaAqzaUckH2ytYX1KNJ8YwbXg6swoyuO7kAk0FJyIigIKzDGTBFtjwrDPCnDYCTrgBhkxyuyqJAJX1zby/vYIPtlewpayOpDgPx43M4MzJarsQEZH2KTjLwNNQASsfdBYpmXieM42cL83tqsRFtf4WlhdV8M6W/WwqrSU9MY7ZBZlcMXsEY3KSiYnREtUiItI5BWcZOCqL4L2FULEdCq+Dk78DMRo5HIz8LUFW7ajkna37+Wh3NQleD7MLMrm0MJ/xQ1IUlEVE5KgoOEv0270S3v0zGAMnfh2GH+92RdLPAsEQa4qreWdLOR/uqiLGGI4f5bRe3HzGOGLVeiEiIr1AwVmik7Ww7Q1nhDnzGDj9vzQ7xiCzu6qRtz4pY9mWcur8AWbkp3HKuByun3cM8bG6mE9ERHqfgrNEF2th86vw/t0w7Fi48G5IynK7KukHjc1B3tu+n7c+KWPLvjqGpvmYOz6H2y+cSnpinNvliYjIIKDgLNHBWmfu5Q/ug5EnwqUPQEKG21VJH7LWsmVfHUs2lfH+9v2A4cRjMrli9kjGDUnGGPUpi4hI/1JwlshmLWx8zgnMBXPhCw9rhowBrCkQ5IPtFSzesI8t++oYOySZ+RNy+PJJo/B51X4hIiLuUnCWyLVtCSz7HYw+Bb74N4hPdrsi6QPldU28sXEfSzaV0dAc4IRjsrjyhJGM1aiyiIhEGAVniTy7V8Kbd0LOBLj0L5CY6XZF0oustWzcW8viDaUsL6okLcHLaROH8IsLp5KRpF5lERGJXArOEjnKPoEldzi9y+f9FlKHul2R9JJQyLJqZyUvrdvLhr01TMpL5fRJudwwb4xW6RMRkaih4Czua6iAJQugYT+c9lPIGuN2RdILWoIh3t9WwUvr91BU3sBxI9O5+Lh8Jg1NUQuGiIhEJQVncU+wBZbfD5+8DPN+AKPmuF2R9JC/JciyzeW8tH4v+2qbOKEgk6+cXMCYHPWni4hI9FNwlv5nLWx+Bd75I8z8Elz1lJbGjmJNgSBvfVLO82tLqG5s4ZRxOdx8xjjyMxLdLk1ERKRXKThL/9q/FV79GQyZBFf8QzNlRKlAMMS72/bz7JoSSmuamDs+hx+dO4ncVJ/bpYmIiPQZBWfpHy1+WPZbKF0HZ/0CMgvcrki6KRSyrNxZybNrStheXs+Jx2Rx02fGMTJLI8siIjI4KDhL39uy2AnNs/8D5v8QdGFYVNmwp4ZFH+5mfUk1x4/K5MsnjmJcborbZYmIiPQ7BWfpO7Wl8MpPIXkIXPF3iFfYihZltU08vXo3SzaVMXZIMhcfN5wfnTNRs2GIiMigpuAsvc9aWPs4rPk7nPn/YOh0tyuSLvC3BFm8YR/PrNlNrCeGC2YM4y/XziIuVhduioiIgIKz9LaaEnjxVicsX/kkeLxuVyQdsNby4a4qnlxZzK6KBs6cnMuCi6drBT8REZE2KDhL77AWPnwM1v0Lzr7DmTVDIlZVQzOLPtzNqx+XMm14GtedPJqxQ9RKIyIi0hEFZ+m52lJ4/ruQPys8yqw/VpHIWssH2yt4fPkuavwBLj5uOA99ZbZaMURERLpICUd6ZtNL8O6f4JxfQe4Ut6uRNuyva+KpVbt5Y9M+jhuZwXfOHM+ITE0hJyIi0l0KznJ0mhuchUxiPHDlE+BNcLsiacVay6qdlTz23k4amgNcclw+j1w3m1iPRpdFRESOloKzdN+etfDyj+Hkb8O4M92uRlrxtwR5Zk0Jz6wuYcqwVL6r0WUREZFeo+AsXWctLL8fti2BSx+E5By3K5Kw4soGHntvJ+t2V3P+jGHcd3UhCXEet8sSEREZUBScpWua6pwLAHMmwhcehRj9yt9t1lre3bafv763E6/HcNWJo7j17AlapERERKSPKDhL5/ZtgBdugXm3QsGpblcz6LUEQ7zw0R6eWFHMjBFp/Ozzk8lN9bldloiIyICn4CwdW/sEfPQEXHI/pOS5Xc2gVtcU4B8f7GTxhn2cPTWPe68+nsQ4/S8sIiLSX/S3rrQtFITXbgMbgi/+TXMzu2hvtZ8H39nOxyU1XDF7JI997QQ8MWrHEBER6W9KQ/Jp/mr49zdg4udg5pfcrmbQ2lpWx91LtlLXFOArJxfww7Mnqn9ZRETERQrOcrj9W+HZb8Pp/w0jZrldzaC0cW8NC5dsxRNj+Pq8MYzL1VLYIiIikUDBWQ7ZtgSW/Q4uugfShrtdzaCztriKe97cRoovlu+cMZ7R2UlulyQiIiKt9Cg4G2MuA24DJgGzrbUr2tnvbOD3gAe431q7oCfnlT6w+u+w+WWnnzlOC2b0pxVFFdz71jby0nz8+HOTGJ6uVRhFREQiUU9HnNcBFwP3tLeDMcYD3AWcCRQDy40xz1hrP+7huaU3WAtLfwN1ZXDJA84S2tIvVu+q4k+vb2FUViK/uHAqQzSlnIiISETrUXC21m4AOrtgaTawxVq7LbzvP4ALAAVntwUD8ML3IfMYOOdXoAvP+sWGPTX88fXNZCfHc/tFUzUHs4iISJTojx7n4cCuVtvFwAn9cF7pSHM9LLoBJl8I0y51u5pBYWtZHX9cvJn4WA8/OmcSIzLVEiMiIhJNOg3OxpjXgLZWvviJtfbpLpyjrWFM28H5rgeuBxg5cmQXDi/d5q+GJ78KJ38LCua6Xc2AV1zZwB8Xb6EpEOSbp49jTE6y2yWJiIjIUeg0OFtrz+jhOYqBEa2284GSDs53L3AvQGFhYbsBW45SQwU8eR185ieabq6PVTe28Oc3tlC0v56bzxjPpKGpbpckIiIiPdAfrRrLgXHGmAJgN/BFQKtquKG2FJ76Gpx1Owyd7nY1A1ZzIMRf39/B4g37+Pr8Mfxo7CS3SxIREZFeENOTNxtjLjLGFAMnAc8bY14OPz/MGPMCgLU2ANwEvAxsAP5prV3fs7Kl26p2wb++Cuf+RqG5j1hrefGjPVzzlw9I9Xl55LrZnDw22+2yREREpJcYayO3G6KwsNCuWNHm1NDSHVU7YdHX4YI/OjNoSK9bvauK3776CbMLMrnu5AIS4jStn4iISLQwxqy01hZ2tp9WDhzoqnfDv78BF94FGaPdrmbAKa9r4jcvbwLgN5fNICcl3uWKREREpK8oOA9ktXudKec+/3uF5l4WCIZ45N0dvLFpH987awIzR6S7XZKIiIj0sR71OEsEqyuDf30NzvstZI1xu5oB5Z2t5Vz9lw9Ijo/l4a/MVmgWEREZJDTiPBA1VoYvBLwTsse5Xc2Asbfazx0vbiArKZ6FVx1PWoLX7ZJERESkHyk4DzTNDfCv/4Az/weGaBq03hAKWf76/g5e3bCPH587kYl5mo9ZRERkMFKrxkASDMC/b4Q5N8GwY92uZkDYuLeGax9aDsbw0LWzFJpFREQGMY04DxTWwvPfgckXwjHz3a4m6vlbgvxh8WZ2VjTw60umk5fmc7skERERcZmC80Cx+OeQOxWmXux2JVHv7S3l/PH1zVx3cgE/OHui2+WIiIhIhFBwHghW/AWwcMINblcS1eqaAtzxwgY8MYb7ri4kxaeL/0REROQQBedot20JFC2Di+93u5Ko9s7Wcv6weDPfOm0cc7RMtoiIiLRBwTmalW+GZb+FL/4dYnSd59FoaA7w65c20RwMaZRZREREOqTgHK0aKuC578DF90FcotvVRKXlRRX87yub+Pr8scwbn+N2OSIiIhLhFJyjUaDZWUr7rF9A6lC3q4k6LcEQv3vtE/bVNHHPlwu1kImIiIh0iX6/H41e/S+Y+SUYNtPtSqLOjv31XPfQcibkpXLnZTMUmkVERKTLNOIcbT56EmJiYcpFblcSVay1LPpwN0+vLuGXF01jRKbaW0RERKR7FJyjSenHsPZx52JA6bIafwu3Pb2eEZmJPHBNIbEe/aJFREREuk/BOVr4q+HFH8AlD4BH/9m6as2uKu54cQPfO2sCs0Znul2OiIiIRDElsGhgLTz7bTjtp5CS63Y1UcFay2Pv7+TtzeXcfdXxpCfGuV2SiIiIRDn9zjoavH835M+CkSe6XUlUqG8K8L0n1lDT2MKfrzxOoVlERER6hUacI92etc7KgJc/5nYlUWHLvlr+69/r+ebpY5kzRisAioiISO9RcI5kzfXw8o/h0gfBGLeriXjPrinhX6uK+d0XZ5Kb6nO7HBERERlgFJwj2cs/gVO+A8la1a4jwZDlVy9txFrLfVcX4tWsGSIiItIHFJwj1cfPQFwSjD3d7UoiWnVjC7c8sYazp+Zx8XH5bpcjIiIiA5iCcySq2wfL74Mrn3S7koi2tayOny5ax63nTGTmiHS3yxEREZEBTsE50lgLL94KZ/0CYuPdriZivbFpH39Ztl39zCIiItJvFJwjzbp/Qc4EGDrD7UoikrWWB5ZtZ+PeWu67uhCf1+N2SSIiIjJI6CqqSFK7F1Y9Aqd+z+1KIlIwZPmfZz+msTnInZdOV2gWERGRfqUR50hhrbOk9md/CR6v29VEnIbmAN9/Yg1nTs7lomN1EaCIiIj0PwXnSPHx05AzEfKmul1JxCmrbeK7/1zN1+eP0aImIiIi4hoF50jgr4Hl92sWjTZs2VfLT/+9jp9fMJXxuSlulyMiIiKDmIJzJHj9FzD3FvBqdojWVu2s5P9e+YQ/fPFYhmjmDBEREXGZgrPbdq+E5jo4Zp7blUSUpZvLePDtIhZedRwpPvV8i4iIiPsUnN0UDMDin8PF97tdSUR54aM9PLumhD9feZxmzhAREZGIoeDsphV/gSkXQXKO25VEjMeX7+SD7ZX84Ypj8Xo0W6KIiIhEDgVntzRUwKbn4aqn3K4kYtz71lZKqvzceel0YmKM2+WIiIiIHEbB2S1LFsDcH0CMWhGstfz2tc1gLf/9+ckYo9AsIiIikUfB2Q1lm6BhP4w+2e1KXGet5X9f+QSfN4abThvvdjkiIiIi7VITqRsW/xxO/5nbVbjOWsuvX95EQpyHm04b53Y5IiIiIh1ScO5vn7wCORMgY5TblbjKWsuClzaS6vPyn58Z63Y5IiIiIp1ScO5PoSC8+0c45btuV+Iqay13vLiRrKQ4vj5/jNvliIiIiHSJgnN/WvtPmHgexCe7XYmrfvPKJrKS4rh+rkKziIiIRA8F5/4SaIbVf4Xjr3W7Elfd9cYWvJ4Ybpin0CwiIiLRRcG5v6x6GGZcAbHxblfimofe3k51YwvfPl0XAoqIiEj0UXDuD831sOEZmH6525W45p8rdrGlrI4fnTNR8zSLiIhIVNI8zv3hg3th1tfAMzi/7ufWlvDB9gp+fcl0hWYRERGJWhpx7mv+atj2Jkw63+1KXLFsczkvfrSXBRdP0zLaIiIiEtUG5xBof1p+P8y+HgbhSOv6kmruW7qNe758PLEe/RtNREREopvSTF9qbnBGm8ef7XYl/W5XRQO/eG4Dv7t8Jj6vx+1yRERERHpMwbkvrXoYjrsaYgbX11xR38wPnlzLnZdNJyMpzu1yRERERHpFjxKdMeYyY8x6Y0zIGFPYwX5FxpiPjDGrjTErenLOqBFogo3Pw5SL3K6kXzU2B/nO46v5r/Mmk5+R6HY5IiIiIr2mpz3O64CLgXu6sO9nrLXlPTxf9Fjzd2f6uZjB06YQClm+/+Qabph7DJOHpbpdjoiIiEiv6tGIs7V2g7V2U28VM2CEQvDRk4Nu3uY7X9nE3HHZzBmb7XYpIiIiIr2uv5pvLfCKMWalMeb6jnY0xlxvjFlhjFlRVlbWT+X1sk9egrGnQ+zg6e99cmUxgWCIy2eNdLsUERERkT7RaauGMeY1IK+Nl35irX26i+c52VpbYowZArxqjNlorX2rrR2ttfcC9wIUFhbaLh4/sqx6GC662+0q+s3yogoWbyjlT186zu1SRERERPpMp8HZWntGT09irS0J3+8zxiwCZgNtBueot2ctpOVDQobblfSLXRUN/O61T7jny4V4tMCJiIiIDGB93qphjEkyxqQceAychXNR4cD0/j1wwo1uV9Evav0t3Pqvtfzqkukkx2stHRERERnYejod3UXGmGLgJOB5Y8zL4eeHGWNeCO+WCywzxqwBPgCet9a+1JPzRqzaUmiqhuxxblfS56y1/PCpj/jumeM17ZyIiIgMCj0aJrTWLgIWtfF8CXBu+PE2YEZPzhM1Vj4IhV91u4p+sfDNrcwenUnh6Ey3SxERERHpF4NrSbu+FAzAjrfhmPluV9Lnlm0uZ+u+eq4+aZTbpYiIiIj0GzWm9pYtr8LYM8AM7Avkdlc1svDNLdx/9SzMAP+sIiIiIq1pxLm3rP4rzLzS7Sr6lL8lyA//tZZfXjSNhLjBsyKiiIiICCg4947qYoj1QdLAXjHvf579mGvnjGZUVpLbpYiIiIj0OwXn3rDqUTjuarer6FNPrSomKymO0yflul2KiIiIiCsUnHsqGICd78DoU92upM9sK6vjmTUl3HzGwJ9mT0RERKQ9Cs49te0NGHPagL0osCkQ5GdPr+eXF00j1qM/LiIiIjJ4KQn11EdPwLQvuF1Fn7njhY1cM2c0w9IT3C5FRERExFUKzj3RVAuNVZA23O1K+sSrH5cCcOZk9TWLiIiIKDj3xMbnYdLn3a6iT+yt9vPwO0X86NyJbpciIiIiEhEUnHvi42dg8vluV9HrQiHLfz29jtvOn0x8rOZrFhEREQEF56NXWwpeH/jS3K6k1z363g7mjMli7JAUt0sRERERiRgKzkdr3b9g6iVuV9HrtpbVsXRzGdecNNrtUkREREQiioLz0dryKow90+0qelVLMMRtz6zn5xdMJSZmYE6vJyIiInK0FJyPRmURpA6D2Di3K+lVd72xhUuOy9fUcyIiIiJtUHA+Ghueg0kXuF1Fr/qouJptZfVcMHOY26WIiIiIRCQF56OxbQkcM8/tKnpNSzDEgpc28N+fn4wZoCsgioiIiPSUgnN31e6FhAyIjXe7kl5z95KtfKFwBFnJA+cziYiIiPQ2Befu2vgcTDrP7Sp6zebSWj7eU8P5M9SiISIiItIRBefu2rIYxp7hdhW9Ihiy/L/nN/AztWiIiIiIdErBuTv81RATC3FJblfSKx5+p4jPTsllaJpm0RARERHpjIJzd2xbAmM+43YVvWJXRQPLtpRzxayRbpciIiIiEhUUnLtj8ysDZtGTX76wgZ98bpIWOhERERHpIgXnrrIWqndD+gi3K+mxVz8uZUxOMmNykt0uRURERCRqKDh31d61kDfN7Sp6rKE5wAPLtvGfnxnrdikiIiIiUUXBuas2vwrjznK7ih770+tb+I9TjyEhzuN2KSIiIiJRRcG5q3a+ByNPdLuKHtmyr5YdFQ2cPinX7VJEREREoo6Cc1f4a8CbAB6v25UcNWstd7ywkR+dM9HtUkRERESikoJzV+x8F0ad7HYVPfLs2j0cNyqD/IxEt0sRERERiUoKzl2x/S0oONXtKo6avyXI397fwddOLXC7FBEREZGopeDcFfs2QM4kt6s4ave9tY1r54wmPlYXBIqIiIgcLQXnzjRUgC8VYqLzqyqt8fPhrio+OyXP7VJEREREolp0psH+tOOdqO5v/s3Lm/jumeMxRisEioiIiPSEgnNnipZCwVy3qzgqa3ZV4YkxTB2e5nYpIiIiIlFPwbkzZZsge7zbVXSbtZb/e/UTvntW9NUuIiIiEokUnDvSWAkJGRCFbQ4vfLSX2QWZDEnxuV2KiIiIyICg4NyR3Sshv9DtKrqtJRji0feKuO5kTT8nIiIi0lsUnDtSvALyZ7ldRbc9saKY82cMJyFO08+JiIiI9BYF547sWQN5092uolsam4M8u6aELxTmu12KiIiIyICi4NyeUAiCzeCNrh7hh98t4qoTRxHr0X9aERERkd6kdNWeiq2QNdbtKrqluqGFZZvLOWeqFjsRERER6W0Kzu0pXh51/c13v7WVG+YdQ0xM9M0CIiIiIhLpFJzbU7w8qmbUKK3x88neWk4Zm+12KSIiIiIDkoJzeyp3QPoot6vosj+/sYX/PG2sltYWERER6SMKzm1pbgBvQtQsfFJa46ek2s9xIzPcLkVERERkwFJwbsu+DZA7xe0qumzhkq3cOG+M22WIiIiIDGgKzm3ZuxbyprldRZfsq/FTXNnI8aM02iwiIiLSlxSc27J3bdQsfHLPW9u4cd4xbpchIiIiMuD1KDgbY+40xmw0xqw1xiwyxqS3s9/ZxphNxpgtxpgf9uSc/aJqF6SPdLuKTpXVNrFjfz2FozPdLkVERERkwOvpiPOrwFRr7XTgE+BHR+5gjPEAdwHnAJOBK4wxk3t43r4TCoKJiYoLA+99ayvXz1Vvs4iIiEh/6FFwtta+Yq0NhDffA/Lb2G02sMVau81a2wz8A7igJ+ftU/u3Qlbkh9Hyuia2ldUzu0CjzSIiIiL9oTd7nK8DXmzj+eHArlbbxeHnItLaj9ex1kb+UtsPLNvO105Vb7OIiIhIf+k0OBtjXjPGrGvjdkGrfX4CBIC/tnWINp6zHZzvemPMCmPMirKysq58hl5VPexUVqSc1u/n7Y5afwvrS2o48RiNNouIiIj0l9jOdrDWntHR68aYa4DzgNOttW0F4mJgRKvtfKCkg/PdC9wLUFhY2G7A7isZiXFUNjT392m75e8f7ORLs0dolUARERGRftTTWTXOBm4FzrfWNrSz23JgnDGmwBgTB3wReKYn5+1LGUmRHZybAkHe2FjGmZPz3C5FREREZFDpaY/zn4AU4FVjzGpjzN0AxphhxpgXAMIXD94EvAxsAP5prV3fw/P2mYxEL5UNLW6X0a6nV5dw/sxheGI02iwiIiLSnzpt1eiItW1fRWetLQHObbX9AvBCT87VXxK8Hhqbg26X0aZQyPLvD3fzl2tnuV2KiIiIyKCjlQOPEMl9w4s37uPUcTn4vB63SxEREREZdBSco8jf3t/Bl06I/BUNRURERAYiBec2ROKY89riKkZnJ5GW4HW7FBEREZFBScG5Db4I7HN+6J0irp0z2u0yRERERAYtBec2ZCR5qYigKen21fjxtwQZlZXkdikiIiIig5aCcxsyEuOorI+c4PzYezu46oRRbpchIiIiMqgpOLchPYJWD/S3BFm5s5KTxmS5XYqIiIjIoKbg3IbMpMhZBOWZNSWcP2NYRE+TJyIiIjIYKDi3IT0xjqoIGHG21vLM6hIumDnc7VJEREREBj0F5zZkJMZREQE9zsuLKpmen6YFT0REREQigIJzGzIT46iKgFaNv3+wUwueiIiIiEQIBec2pCd5XR9xLq9roqE5QH5Goqt1iIiIiIhDwbkNKfGx1DUFXK3hyZXFfKFwhKs1iIiIiMghCs5tcHsGi1DI8tYnZcyfMMTVOkRERETkEAXnCLRsSzlzxmThidEUdCIiIiKRQsE5Aj2+YpfaNEREREQijIJzO2KMIRAM9ft591b7iTGGIam+fj+3iIiIiLRPwbkdaQleavz9f4Hg48t3cblGm0VEREQijoJzO9ITvf2+emAoZHl/+37mjMnq1/OKiIiISOcUnNuRkeilsp8XQXlv235OKMgiRhcFioiIiEQcBed2pCXGUd3YvyPOT64q5pLjh/frOUVERESkaxSc25Ge4O3XZbdr/C3UNLZopUARERGRCKXg3I6MxLh+bdV4bs0ePj9jWL+dT0RERES6R8G5HemJXqr78eLAl9bv5bNT8vrtfCIiIiLSPQrO7UhL8FLV2D8jzlv21TIiIwGf19Mv5xMRERGR7lNwboczHV3/BOcnVhRzmeZuFhEREYloCs7tSI6Ppa6p7xdACQRDrCupZkZ+Wp+fS0RERESOnoJzO4zpn7mU39m6n5PHZvfb+URERETk6Cg4u+zp1SWcr9k0RERERCKegrOLGpuDVNQ3ae5mERERkSig4NyBGGMIBEN9dvzXNpRy+qTcPju+iIiIiPQeBecOpCV4qfH33QWCL3y0h89NG9pnxxcRERGR3qPg3IGMRC9VfbQIyoHjZiTF9cnxRURERKR3KTh3ID3R22fLbj//0R7O0WiziIiISNRQcO5AWmIc1Y19M+L82selnKn+ZhEREZGooeDcgfSEvlk9cHdVI5lJ8STEaYltERERkWih4NyBjMS4PmnVeH5tCefNUJuGiIiISDRRcO5AeqKX6j64OHDZlv2cPCa7148rIiIiIn1HwbkDaQleqhp7d8S5uLKBISnxxMXqqxcRERGJJkpvHchIiuv1HueX1u3l3Gl5vXpMEREREel7Cs4dSIrzUNfUuwugLN1czslj1aYhIiIiEm0UnDtgjOnV45VUNZKVHEd8rGbTEBEREYk2Cs6d6M3o/OK6vZw7VbNpiIiIiEQjBedO2F481tLNZZwyTm0aIiIiItFIwbkTMcYQCIZ6fJy91X4yEuPwedWmISIiIhKNFJw7kZbgpcbf8wsEX1y3h7OnajYNERERkWil4NyJjEQvVb2wCMrSzeXMG5/TCxWJiIiIiBsUnDuRnujt8bLbVQ3NJHg9atMQERERiWKxPXmzMeZO4PNAM7AV+Iq1tqqN/YqAWiAIBKy1hT05b39KS4yjurFnI85LNpUxb4JGm0VERESiWU9HnF8Fptr/397dxfZV13Ecf3/Xh21tB13ZA4OuwuRhIPLkQHBGGSNIEMEbEhI1xIfghVE0JipyYeKVF8ZoYqIh+EAiwRCEYBBxMGYIFzzIgAB2OIVRYGPtNjrWwcrafr3of8ZB/9vp2v7PWXm/bv7/nnO282nzy/6fnf7O+WWeDfwLuOkQx67JzHOPptIMB6ZqTO2K88Ob+rl05ZJpSiRJkqQyTKk4Z+a6zDxw59xjQPfUI1VL5/zWKU3VeHdkjN3v7B7hGfUAAAeuSURBVGdRx9xpTCVJkqRGm845zl8B/lpnXwLrIuKpiLhhGs854zrbWtg9hZsDn9yyiwtP7prGRJIkSSrDYec4R8RDwETPUbs5M++tHXMzMALcXuevWZ2ZWyNiCfBgRGzKzEfqnO8G4AaAnp6eAt/CzOpsa2HwnSO/4vxQ73auu6D870OSJElTc9jinJmXHWp/RFwPXAWszcwJF9rLzK211/6IuAe4EJiwOGfmLcAtAKtWrZrOhfuOSGdb6xHPcc5MNm8f4rSlHdOcSpIkSY02pakaEXEF8H3g6sx8u84x7RGx4MB74HLg+amct5HaW5sYGj6yBVD+3T/EKUs6iIhpTiVJkqRGm+oc518CCxiffvFMRPwaICJOiIj7a8csBR6NiGeBJ4C/ZOYDUzxvw0yl9D7Yu521Z/g0DUmSpNlgSs9xzsxT6mzfClxZe/8ScM5UzlO2I63OT7y8i699csW0ZpEkSVI5XDmwgCOZaL1zaJgF81pobfZHLEmSNBvY6gqYE8HI6Nik/syGFwdY42qBkiRJs4bFuYDOthZ2T/KRdBte7GfN6c5vliRJmi0szgUc19HKzr3FF0EZHhll7/AIC9tbZzCVJEmSGsniXMDijrns2DNc+PjHXtrFxSuOm8FEkiRJajSLcwGLF8xlYKh4cV7fu521ZyydwUSSJElqNItzAYs75jJQ8IpzZvKfgSE+vLh9hlNJkiSpkSzOBSyaxBXn3m17WHn8Ma4WKEmSNMtYnAuYzBXn9a4WKEmSNCtZnAs4dn4Lu98u9ji6J195kwtO6prhRJIkSWo0i3MBc+YEY3n49QP79+yjq62FliZ/rJIkSbONDa+gInOWN2zqZ81Kp2lIkiTNRhbngtpam9iz79DTNf7+4gCXnGZxliRJmo0szgWd2Dmfbbv31d2/b/8owyNjHNvW0sBUkiRJahSLc0HLjp3H1sF36u5/dPMOVp+yqIGJJEmS1EgW54KWdc5n62D9K85/e+ENLj/T1QIlSZJmK4tzQT1dbfTtenvCfSOjY7zx1j6Wd7U1OJUkSZIaxeJc0MmL2tmyY++E+57YsosLfXazJEnSrGZxLmheSxPDI6MT7lv3wnY+c9bxDU4kSZKkRrI4T8KcCEZGxw7alpls7t/DqUs6SkolSZKkRrA4T0LPcW1s2XnwPOeNfYOc3d1ZaIEUSZIkHb0szpNw7vJOnu5786Btdz75Ktd+rLukRJIkSWoUi/MknN+zkI19g//7emh4hIGhYVYsdpqGJEnSbGdxnoTuhfN5ZedeMhOAe595nc+ds6zkVJIkSWoEi/MkRARnd3eysW+QfftHue/ZbVz5UYuzJEnSB0Fz2QGONl/4eA833f0cba1NfP3TK5jb3FR2JEmSJDWAxXmSlne18eNrPsK7o2OsPP6YsuNIkiSpQSzOR8CbASVJkj54nOMsSZIkFWBxliRJkgqwOEuSJEkFWJwlSZKkAizOkiRJUgEWZ0mSJKkAi7MkSZJUgMVZkiRJKsDiLEmSJBVgcZYkSZIKsDhLkiRJBVicJUmSpAIszpIkSVIBkZllZ6grIgaAV8rOIRYBO8oOoUpybKgex4bqcWzoUMoaHx/KzMWHO6jSxVnVEBH/yMxVZedQ9Tg2VI9jQ/U4NnQoVR8fTtWQJEmSCrA4S5IkSQVYnFXELWUHUGU5NlSPY0P1ODZ0KJUeH85xliRJkgrwirMkSZJUgMVZE4qI5RGxISJ6I+KFiLix7EyqlohoioinI+K+srOoWiKiMyLuiohNtX9DLi47k6ohIr5T+0x5PiLuiIh5ZWdSOSLitxHRHxHP/9+2roh4MCI2114XlplxIhZn1TMCfDczzwAuAr4REWeWnEnVciPQW3YIVdIvgAcycyVwDo4TARFxIvAtYFVmngU0AdeVm0ol+j1wxXu2/QBYn5mnAutrX1eKxVkTysxtmbmx9n4P4x98J5abSlUREd3AZ4Fby86iaomIY4BPAb8ByMx3M3Ow3FSqkGZgfkQ0A23A1pLzqCSZ+Qiw6z2brwFuq72/Dfh8Q0MVYHHWYUXEScB5wOPlJlGF/Bz4HjBWdhBVzgpgAPhdbSrPrRHRXnYolS8zXwd+CvQB24Ddmbmu3FSqmKWZuQ3GL+ABS0rO8z4WZx1SRHQAfwK+nZlvlZ1H5YuIq4D+zHyq7CyqpGbgfOBXmXkesJcK/rpVjVebr3oNcDJwAtAeEV8sN5U0ORZn1RURLYyX5tsz8+6y86gyVgNXR8QW4I/ApRHxh3IjqUJeA17LzAO/obqL8SItXQa8nJkDmbkfuBv4RMmZVC3bI2IZQO21v+Q872Nx1oQiIhifo9ibmT8rO4+qIzNvyszuzDyJ8Rt7Hs5MrxoJgMx8A3g1Ik6vbVoL/LPESKqOPuCiiGirfcasxRtHdbA/A9fX3l8P3Ftilgk1lx1AlbUa+BLwXEQ8U9v2w8y8v8RMko4O3wRuj4hW4CXgyyXnUQVk5uMRcRewkfEnNz1NxVeJ08yJiDuAS4BFEfEa8CPgJ8CdEfFVxv+jdW15CSfmyoGSJElSAU7VkCRJkgqwOEuSJEkFWJwlSZKkAizOkiRJUgEWZ0mSJKkAi7MkSZJUgMVZkiRJKsDiLEmSJBXwXxsfUuxQ6lhLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Sacar la solución y graficar\n",
    "PC1=pred_corr(f,t_inicial,y_inicial,t_final,h)\n",
    "PC2=pred_corr(g,t_inicial,y_inicial,t_final,h)\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(PC1[0],PC1[1],lw=0.7,label='Solución 1')\n",
    "plt.plot(PC2[0],PC2[1],lw=0.7,label='Solución 2')\n",
    "\n",
    "plt.title('Adams-Bashforth-Moulton')\n",
    "plt.legend()\n",
    "#plt.savefig('grafica.pdf',format='pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}