{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Multi-step methods\n",
    "\n",
    "Our goal is to implement Adams-Bashforth methods for solving linear differential equations. We focus on the 2nd-order method, but include code for the 4th-order method. Adams-Bashforth methods are called linear multi-step methods because they use a linear combination of function evaluations at multiple timesteps (i.e. a weighted average of values of f at previous time steps). These methods are worth implementing because they have the accuracy of the R-K methods but require less computatation at each step. The drawback is that their regions of stability are quite small, so we may need to take small steps in the interest of stability, rather than accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import numpy.linalg\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def degradation_model( t, y, params ):\n",
    "    y1 = y[0]\n",
    "    y2 = y[1]\n",
    "    dydt = np.zeros( (2,) )\n",
    "    dydt[0] = -params[0] * y1\n",
    "    dydt[1] = -params[1] * y2\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forwardEuler( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler( fn, t, y0, params )\n",
    "    Given a function handle fn for dydt, find the value of\n",
    "    y at timepoints t, given y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[:,i-1],params);\n",
    "        y[:,i] = y[:,i-1] + h*yprime;\n",
    "    \n",
    "    return (t,y)\n",
    "\n",
    "def explicitTrapezoidal( fn, t, y0, params ):\n",
    "    \"\"\"[t,y] = explicitTrapezoidal( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[:,i-1],params);\n",
    "        yhatFE = y[:,i-1] + h*yprime;\n",
    "        yprimehatFE = fn(t[i],yhatFE,params);\n",
    "        y[:,i] = y[:,i-1] + h/2*(yprime+yprimehatFE);\n",
    "\n",
    "    return (t,y)\n",
    "\n",
    "def RK4( fn, t, y0, params ):\n",
    "    \"\"\"\n",
    "    [t,y] = RK4( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    Uses a 4th order explicit Runge-Kutta method\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        f1 = fn(t[i-1],y[:,i-1],params);\n",
    "        f2 = fn(t[i-1]+h/2,y[:,i-1]+0.5*h*f1,params);\n",
    "        f3 = fn(t[i-1]+h/2,y[:,i-1]+0.5*h*f2,params);\n",
    "        f4 = fn(t[i],y[:,i-1]+h*f3,params);\n",
    "        y[:,i] = y[:,i-1] + h/6*(f1+2*f2+2*f3+f4);\n",
    "    return (t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implement the method\n",
    "\n",
    "Use the explicit trapezoidal method to compute the value at the first time step after the initial condition. \n",
    "\n",
    "Then, use the Adams-Bashforth method to compute the remaining time steps. Write the code so that the ODE function is called as few times as possible.\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "f_b &=& f(t-\\Delta t,\\widehat{y}(t-\\Delta t))\\\\\n",
    "f_a &=& f(t,\\widehat{y}(t))\\\\\n",
    "\\widehat{y}(t+\\Delta t) &=& \\widehat{y}(t) + \\frac{\\Delta t}{2}\\left(3f_a-f_b\\right)\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adamsBashforth2( fn, t, y0, params ):\n",
    "    \"\"\"\n",
    "    [t,y] = adamsBashforth2( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    \"\"\"\n",
    "    h = t[1]-t[0];\n",
    "    if abs(np.mean(np.diff(np.diff(t)))) > 1e-8:\n",
    "        print( 'The time steps must be evenly spaced' );\n",
    "        return;\n",
    "\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run and plot result from AB2.\n",
    "\n",
    "Run the AB2 method on the degradation model with degradation rates (2,5) and initial conditions (1,1) with time beginning at 0 and ending at 2.\n",
    "Stephanie\n",
    "Then plot the true solution (solid lines) and the forward Euler's solution (dashed lines), explicit Trapezoidal's solution (dotted lines), and AB2's solution (dash-dotted lines).\n",
    "\n",
    "Then plot the difference over time between the true solution and each of the other methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "params = (2,5)\n",
    "y0 = np.array([1,1])\n",
    "dt = 0.1\n",
    "t = np.arange(0,2,dt)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine how the error depends on the time step size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,0.2,10)\n",
    "errFE = np.zeros( (dts.size,) )\n",
    "errET = np.zeros( (dts.size,) )\n",
    "errAB2 = np.zeros( (dts.size,) )\n",
    "params = (2,5)\n",
    "y0 = np.array((1,1))\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the error vs step-size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure()\n",
    "# plt.plot( dts, errFE, 'C0.-' );\n",
    "# plt.plot( dts, errET, 'C1.-' );\n",
    "# plt.plot( dts, errAB2, 'C2.-' );\n",
    "# plt.legend( ['FE Error','ET Error','AB2 Error']);\n",
    "# plt.xlabel( 'Size of Time Step');\n",
    "# plt.ylabel( 'Error');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implement AB4\n",
    "\n",
    "Use the RK4 to compute the value at the first three time steps after the initial condition. \n",
    "\n",
    "Then, use the Adams-Bashforth method to compute the remaining time steps. Write the code so that the ODE function is called as few times as possible.\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "f_d &=& f(t-3\\Delta t,\\widehat{y}(t-3\\Delta t))\\\\\n",
    "f_c &=& f(t-2\\Delta t,\\widehat{y}(t-2\\Delta t))\\\\\n",
    "f_b &=& f(t-\\Delta t,\\widehat{y}(t-\\Delta t))\\\\\n",
    "f_a &=& f(t,\\widehat{y}(t))\\\\\n",
    "\\widehat{y}(t+\\Delta t) &=& \\widehat{y}(t) + \\frac{\\Delta t}{24}\\left(55f_a-59f_b+37f_c-9f_d\\right)\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adamsBashforth4( fn, t, y0, params ):\n",
    "    \"\"\"\n",
    "    [t,y] = adamsBashforth4( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    Uses an explicit 4th order Adams-Bashforth method\n",
    "    \"\"\"\n",
    "    h = t[1]-t[0];\n",
    "    if abs(np.mean(np.diff(np.diff(t)))) > 1e-8:\n",
    "        print( 'The time steps must be evenly spaced' );\n",
    "        return;\n",
    "\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare AB4 output to RK4 output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = (2,5)\n",
    "y0 = np.array([1,1])\n",
    "dt = 0.01 # 0.1 is unstable for AB4\n",
    "t = np.arange(0,50,dt)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine how the error depends on the time step size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,1e-2,10)\n",
    "errRK4 = np.zeros( (dts.size,) )\n",
    "errAB4 = np.zeros( (dts.size,) )\n",
    "params = (2,5)\n",
    "y0 = np.array((1,1))\n",
    "\n",
    "# plt.figure(figsize=(8,12))\n",
    "# for i in range(dts.size):\n",
    "#     t = np.arange(0,50,dts[i])\n",
    "#     sol = scipy.integrate.solve_ivp(fun=degradation_model, t_span=(0,t[-1]), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "#     (t,y) = RK4( degradation_model, t, y0, params)\n",
    "#     plt.subplot( dts.size, 1, i+1 )\n",
    "#     plt.plot( t, sol.y.T )\n",
    "#     plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "#     plt.plot( t, y.T, '--')\n",
    "#     errRK4[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "#     (t,y) = adamsBashforth4( degradation_model, t, y0, params)\n",
    "#     plt.plot( t, y.T, ':')\n",
    "#     errAB4[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "# plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the error vs step-size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure()\n",
    "# plt.plot( dts, errRK4, 'C0.-' );\n",
    "# plt.plot( dts, errAB4, 'C1.-' );\n",
    "# plt.legend( ['RK4 Error','AB4 Error']);\n",
    "# plt.xlabel( 'Size of Time Step');\n",
    "# plt.ylabel( 'Error');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Code\n",
    "### Illustrate the trapezoidal rule and AB2 in terms of numerical quadrature\n",
    "\n",
    "Plot the method computing the area under a known curve (not an ODE)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAElCAYAAAChykXrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABT+klEQVR4nO3dd3gc1dXH8e9R75Ll3iu2cQfbgG1MJ4ABG0zvLQFeSEjovYeWAHHoJITeQ+8EMMZUg417x703ySouqvf9Y0awllVWtla7K/0+z6PHuztlz87sHp+ZuXeuOecQERERkcgQE+4AREREROQ3Ks5EREREIoiKMxEREZEIouJMREREJIKoOBMRERGJICrORERERCKIirMmwMxGmtn8GqY/a2Z/3cV1LzWzw3Y9ul1jZn81s41mtrah31tEQk95KzqZ2Rlm9r9wxxHtVJwFycwKA/7KzWxbwPMzwh1fTZxzXzvneu3KsmaWYGYPmNlK/7MuNbNx9RxiXWPqBFwJ9HHOtQlTDD3N7L9+os0zsxlmdoWZxYYjnsrMzJlZjxqm17hfw/Wfl9Qv5S3lrYD3v83MSvztsdnMvjOzYfX9Ps65l5xzv6vv9dbEzPb3P0+emeWY2bdmNtSfdq6ZfdOQ8dQHFWdBcs6lVfwBy4FjA157qWI+M4sLX5QhcT0wBNgHSAcOAn4OZ0BAJ2CTc259VRNDvQ/MrDswCVgB9HfOZQIn4W2n9F1YX1yl52Zmof5tRuJ+lXqmvBVR3++w5i3fa/53oSXwDfCWmVkVsUTEQWYwzCwD+AB4GMgG2gO3A0XhjGu3Oef0V8c/YClwmP/4IGAlcC2wFngBaIb3ZdkA5PqPOwQsPwG4B/gRyAfeBbIDpu8HfAdsBqYDB/mvDwMKA/62A0v9aYnAOGC1/zcOSAyMMWD9e+ElqgLgNeBV4K/VfNYPgL8EuS1qjQG4AdjoL3dGwHoSgfvx/gNZBzwBJFfxfocB24Byfxs8C3QBHHCBv/xEvAOPm4BlwHrgeSDTX0fF/OfhFVi5wMXAUGCGv90fqeEzvwh8WMP0HbZ3FdvpNuANfz35wO/978RdwLf+5+sB9AY+A3KA+cDJAet7FngU+NDfj5OA7v60if7n2+Jvo1Pqsl/xvsPlfhyFwDU1fS+D+U7rL/x/KG9Vty2aSt66DXgx4Hlff30t/HgeBz7CyxuHAe2AN/3vwxLgMn+5dv5nya60bzYC8cC5wDcB04YDPwF5/r/Dq9oPlWMEkvBy5Cb/s/0EtK7icw0BNlfzmffE+76V+dt9c237rbZ93mC/13AnjGj8Y+ckVwrc5+/wZKA5cAKQgnfU9l/gnYDlJwCrgH5Aqv8DqPhCtve/jKP8H+rh/vOWlWKIB74C7vGf3wH8ALTCOyr6Drgz8MvmP07A++Ff7q/jRKCE6pPcTf4X+BKgP2A1bIvaYigFHvS304F4SaCXP/0fwHt4Rz7pwPsVn62KmH79PP7zLnhJ5nl/eyYD5wO/AN2ANOAt4IVK8z+BlwB+h/cDfsePvT1eYjywmvdfC5xXw/djh/iq2E63+dv8OH8fJ/vfieV4CTMOyMRLwOf5zyuSXx9/Hc/634t9/OkvAa8GvJ8DetQQY9D7NZjvJTV8p/UXGX8ob1W3LZpK3rotYH8lAn8HlvvPn8Urnkb4+y8FmALc4m/7bsBi4Ah//vHAHwLW/XfgCf/xufjFmb9dcoGz8PLUaf7z5pX3QxUxXuRvzxQgFhgMZFTxuTLwvmvPAUcBzSpN/zWegNeq3W+17fMG+72GO2FE4x87J7liIKmG+QcBuQHPJwD3Bjzv468jFu9I9oVKy38KnFPptcfxjg5j/OeLgFEB04/gt6PTg/gtyR2Ad3RoAfN+R/VJLha4FO+MTpG/7Dmu6m1RWwylQGrA9NeBmwHzv/zdA6YNA5ZUE9Ovn8d/3gUvaXULeO0L4JKA573wknlcwPztA6ZvIuAME95/PH+p5v1LgCNr2N87xFfFdroNmFhp+gTgjoDnpwBfV5rnSeBW//GzwFMB00YB8wKe11acBb1f/ec1fi+p4Tsdyt+i/oL/Q3nrHFf1tmgqees2f39txivixgOD/WnPAs8HzLsvfuEW8Nr1wDP+498D4/3HhncgeYD//Fx+K87OAn6stJ7vgXMr74eAGCuKs/P9fTwgiO/2nv5nWOnvr/fwz7Kx85m8GvdbTfu8IX+vja2dQbhscM5tr3hiZil4lfmReJcKANLNLNY5V+Y/XxGw/DK8o8EWQGfgJDM7NmB6PPBlwPovwvsC7eucK/dfbuevJ3Cd7aqItR2wyvnfuIB5q+TH+yjwqJlVHNk9bWY/OufmVrHummLIdc5tqWJ6S/wjtYDmD4aXYOsicJtWFUsc0DrgtXUBj7dV8TytmvfZBLStY2yVrajltc7Avma2OeC1OLzLTxUCe3xtpfp4d1LH/VoRT43fS6r/TgduV4kcylvBxdBY8hbA6865M4OIozPQrlL+iQW+9h+/CTxsZm2BnniXa79mZ5U/D/7z9jXEWOEFoCPwqpll4V3ivNE5V1J5Rn+fngtgZr39ecfhnamrLJj9Vt0+bzDqEFA/XKXnV+Id8ezrnMvAO+oD7wtQoWPA4054R0cb8X4gLzjnsgL+Up1z94LXvRy4ExjjnMsPWMdqvB9U4DpXVxHrGqB9pUagnYL6kM5tc849indauk8Vs9QWQzMzS61i+ka8pNI34DNnOq/hal0E7oeqYimlfgqFz/Eu/1RnC96PH/i1cW3LSvNU/s5Ufm0F8FWl70Gac+7/djXo6lSzXyvHV+P30lfdd1oik/JWcDE0lrxVlzhW4J1JCtyf6c65UQDOuVzgf3hn+E/Ha1JRVU6r/HnA+0yr/Mc75Erg116szrkS59ztzrk+eO3WjgHOrvVDODcP7yxavyo+FwS336rb5w1GxVlopOPt/M1mlg3cWsU8Z5pZH/9o9Q7gDf9o70XgWDM7wsxizSzJzA4ysw5m1hHv9OrZzrkFldb3CnCTmbU0sxZ4bQVerOJ9v8f7sV9mZvFmNhav3VKVzOwv/vsnm1mcmZ3jf76pVcweTAy3+93cR+L92P7rH0X/G/iHmbXy37e9mR1RXVxBeAW43My6mlkacDdeT6XS3VhnhVuB4Wb2dzNr48fbw8xe9I/wFgBJZna0mcXjtX9JrON7fAD0NLOz/P0Ub2ZDzWzPIJdfh9dOpEpB7NfKy1f7vQyYp7rvtEQH5a3Gnbfq4kegwMyu9bdhrJn1M//2FL6X8YqlE/3HVfkIL4+d7u+HU/AK5A/86dOAU/19OsRfFwBmdrCZ9fcPbvPxDgTKqcTMepvZlRW5yP++nYbXjhC8XNbBzBIA6rDfdtrntWyzeqXiLDTG4TXu3Ij3BfmkinlewKvu1+I17rwMwDm3AhiD11NkA94RzNV4++pQvNPbb9hv9yqa7a/vr8BkvF47M/F6Ne10g0bnXDEwFu8UcA7ekc9bNXyWrcADfpwb8dpxnOCcW1zFvLXFsBbv6HU1XgP2i/2jHPDarPwC/GBm+Xhnp3bpHke+p/G28US8nkbbgT/txvp+5ZxbhNdGoQsw28zy8E7zTwYKnHN5eA2Rn8I7QtyC1xaiLu9RgNfg91S87bWW3xpvB+M24Dnz7md0chXTa9uv9+D9h7XZzK6q5XtZocrvtESNcShvNdq8VRd+wX0MXrvDJXjb8Cm8jkoV3gP2ANY656ZXs55N/nquxGsOcg1wjHOu4oz6zUB3vO17OzsWeW3werXnA3PxOpIENuuoUIDXRm6SmW3B++7O8t8TvLZ1s4G1ZlbxvrXtt5r2eYOwqs9ESiiZ2QS8Ro9PhTuWhmJmB+F95g61zCpRqCl+p5uapriPlbeankjZ5zpzJiIiIhJBVJyJiIiIRBBd1hQRERGJIDpzJiIiIhJBVJw1QWZ2j5n9Jdxx7C4ze8LMbq5hujOzHkGsp4s/b73flNnMbjOzqm4NUNW8b5rZUfUdg0gohSufmNmzZrZTz05/2rlm9k1Dx+S/d9C/+WhT23Y1s9ZmNtHMCszsgXp6zxq3p5n9aGZ96+O9IomKsybGzFri3ZvmSTM7I6Br+zYzKw94XhjuWGvjnLvYOXdnqN/HzJb626fQzNb6/ynU9UaTwbiPKm4jIBKpAvOJ/7zK/7z939BhDR1fffJ/96Xm3RU/ovkFTUlAPp9rZjXdOLu+XIh3240M59yVtc1cmX9vujrddghvAPM76vpekU7FWdNzLvCRf9fsl/y7zqfhDRi7uuJ55btc+zcCbMqO9bfJILxByK+v7zdwzv0IZPg3YxSJBufi55NwBxJK5t0t/gS8wcGrG/4o0rwWkMv/ArxoZq1rWWZ3dQbmVDNaQI1248rFe8DB5t8QvLFQcdb0HIV3M78a+UeJj5vZR/6N/Q427473U80s38xWmNltAfNXXBq80MxWm9kaM7sqYHqMmV1nZovMbJOZvW7eXcgxs0cCz9j5R6e3+dP2NLMJ/s1QZ5vZ6Eox/jXg+dX++642s/MrfZ5qY68L59xavAGdB/nr3elIr6azBGa2n5l953+e6f49dQJNAI7eldhEwiCofBLIzEaZ2Rz/0teqSnniGDOb5v8+vjOzAQHT9jKzn/3lXsO7CW4tb2WPmFmemc0zs0P9F08ysymVZrzCzN6tYV0n4A0YfgdwTqVlu5rZV35cn+GNNRo4/b/+Gfc88y759Q2Y9qyZPWZmH/u571sza2Nm48ws1497r4D5r/W3WYGZza/4TLVxzn2Kd7PW7v56mpnZB2a2wX+fDyxgtA/zzoAu9t9niZmdUekz3e8vt8T8phhm9qy/ba7xP8thZpbof5bV/t84M0v05z/IzFb6n2kt3ugIH+ON6Vnxf0HFeJYJZva8H89sCziAdd74sFPwBqxvNFScNT39gflBzns6cBfesCff4N3p/mwgC6+A+D8zO67SMgfj3TX6d8C1AUXKn4DjgAPxBpDNxRuYGOfcHwOO8Pb3p71r3tBH7+ON4dbKX8dLZrbTHbjN7EjgKuBw//0rF0fBxF4rP4EdhXd36bou2x74EO/SZbYf75vmXRqqMBcYWNd1i4RJXfJJhf8AFznn0vHGPxwPXvGFd4f8i4DmeJdK3/P/g08A3sG7Q3w23lA6tV2m2xdYhFcs3Qq8Zd4B4XtAV9txKLSzgOdrWNc5eMXDq0BvMxscMO1lvOKgBd74oedUWvZjvJzUCm/0gZcqTT8Zb4i3FkAR3lBVP/vP3wAeBPDz3h+Bof62OwJYWtMG8JczMzsaSADm+C/HAM/gnenqhDds1yP+/KnAQ8BR/vsMxxtmqcK+ePu8BfA34D9mZs65c/3P9jc/n38O3Ajsh3cwOxBvyK2bAtbVBm9/dsbLz5Wv4FSMZzkab9tn4e2/Ryp9zEaXN1WcNT1ZeEdQwXjXOfetc67cObfdOTfBOTfTfz4DL1kdWGmZ251zW5xzM/F+/Kf5r18M3OicW+mcK8IbXuhECziV7Rcp7wB/cs5NxftRpwH3OueKnXPj8cZkO42dnQw845yb5Zzb4q//V0HGXpN3zKwAb1ia9VQ97mBtzsS7BPSRH8dneMPGjAqYpwBvH4lEgyyCzycVSoA+ZpbhnMt1zv3sv34h8KRzbpJzrsw59xxesbKf/xcPjPMHxH4D+KmW91kfMP9reAXF0X7+eQ3/8qR/JqsLv433uAMz64R30Pmyc24d8AX+ANz+tKHAzc65IufcRLwDyl855552zhUE5L2BZhY4DNLbzrkp/hmgt4Htzrnn/SGUXsNrRgFQhjd8Wx8zi3fOLfWHkqvOyWa2GSjEK2juds5t9mPa5Jx70zm31R8q7i52zIflQD8zS3bOrXHOzQ6Ytsw5928/vueAtnjDc1XlDOAO59x659wGvCGazqr0Prf6266mS+Pf+HmzDK9Ar1yINbq8qeKs6cnFOxMWjBWBT8xsXzP70j8VnodXcLWoYZlleGfJwDsyetu/XLEZ70inDP9H7Z8lewMvAb7qL9MOWOEPVBu4zvZVxNquiveua+w1Oc4/ijwI6F3HZSt0Bk6q2Ab+dtgfL7lVSMe7fCISDSrnk1K8IqqyeLyiDLwzXqOAZf7lwGH+652BKyv9Pjri/bbbAasqtWXa4Tdeharmr8hHzwGnm5nhFQuv+8VTVc4C5jrnpvnPX/KXjffXl+sfEO4Ul3kDht9rXnOOfH470xWYP9YFPN5WxfM0AOfcL3htx24D1pvZqxWX/WzHZiGd/GVfd85lOedS8S5nnm1mF/nzp5jZk2a2zI9rIpBlZrH+ZzkFL0euMbMPzax3QExrKx4457b6D6vrINWOHfdT4D4A2OAXpbVZG/B4K5BkO7ZRa3R5U8VZ0zMD6BnkvJUbdb6MdwTW0TmXCTwBWKV5OgY87oQ3cCx4hdNRfrKo+Etyzq3ypz+MN8Bt4Cnv1UBHMwv8nnbCG0y8sjVVvHddY6+Vc+4rvIGf7/df2gKkVEw3r+NEy52XBLxt8EKlbZDqnLs3YJ49gSoHERaJQJXzyXKgk1/0AF4hgHdJbxmAc+4n59wY/7V3gNf9WVcAd1X6faQ4517B+323D1wvO//GK6tq/tV+DD8AxcBIvOYbVQ2oXeFsoJt57cbW4l1mbIFXYK4BmvmXAquK63S8AeEPwxs0vIv/ep1zjx/3y865/fEKWYfXw5uAy4BpzrnlVSy3FO/y6rH+S1fiDfS9r3MuAzggMC7n3KfOucPxDhznAf/elXjxtnfngOeB/yfAzv/H7Opd8Rtd3lRx1vR8RN0u5wVKB3Kcc9vNbB+8xFPZzf5RWV/gPLzT8uAVQ3eZWWfwLmGa2Rj/8UV+TGdUOks2Ce8o6Rozizev8fyxeG0PKnsdONfM+vj/GVS+7BhM7MEaBxxuZgOBBXhHcUf7R9I34V16qMqLwLFmdoR/RJ3kN4oNHGD3QLwkKhINKueTScB24Dr/+50K3It3+X6ZmSWYdwufTOdcCd4BWcVv/t/Axf5ZbjOzVP93lY7XDqsUuMzPBWPx2i/VpFXA/Cfh/Qf+UcD05/HaLpU456q8d5d/Vq+7/16D/L9+eAd7Zzvnlvmf7Xb/s+3PbwUQeHmnCNiEdxB3dy0xV8vMepnZIeY1qN+Od1atvJbFKpbtABwJVFyeTPeX3+y3w7s1YN7WZjbG33dFeJdFg3qfKrwC3OTn+xbALXh5sDrrgOaVLvvWyMySgMHAZ7sYY0RScdb0PA+MMrPkXVj2EuAOv+3VLfx2xBvoK7zG8l8A9zvn/ue//k+8M1f/85f/Aa9hKXhtyLoBqwNOzd/gnCvGS3RH4d075zG8hDiv8ps65z7GK5rG++8/fhdiD4rfduJ54BbnXJ6/7qfwzuhtAaq8T49zbgXeUfQNwAa8MwVX4/8OzWwoUOi8W2qIRIMd8ol/afBovMv/K4HFeJexTg64xHgWsNS/nHYxXrsknHOTgT/gFUy5eL/jc/1pxcBY/3kO3mW3t2qJbRJeQ/yNeG2qTnTObQqY/gJeoVVTsXAOXtvbmc65tRV/ePnsGL+wOR0vl+XgFTmBHQuexztjuAqvMf4PtcRck0S8Qncj3mW+VtR8S59T7Ld7Vv4EfIvX5gu8XJnsr+sH4JOA5WKAK/DOcOXgFd//t4sx/xWveJ0BzMTr6FDtvRz93P4KsNi/tN2uunkDHAtMcL91HmgUNLZmE2RmdwPrnXPj6nGdXYAlQLxzrrS+1tuUmNmbwH+ccx/VOrNIhAhFPmkIfkG5HtjbObcw3PHIrjGzScAFzrlZ4Y6lPqk4k3qh4kxEoomZXQEc45w7JNyxiFRW72MJioiIRDIzW4rX+P248EYiUjWdORMRERGJIOoQICIiIhJBGs1lzRYtWrguXbqEOwwRaUBTpkzZ6Jyr7r5yUUU5TKRpqSl/NZrirEuXLkyePDncYYhIAzKz2u4SHzWUw0Salpryly5rioiIiEQQFWciIiIiEUTFmYiIiEgEUXEmIiIiEkEaTYeAYGwpKuXJiYt48ftl5G4toVlKPGcO68xFB3QnNbFJbQoRiTLKXyJNR0jPnJnZUjObaWbTzGynbkjmecjMfjGzGWa2d8C0c8xsof93zu7GsqWolOMf+5Ynv1pMztYSHJCztYQnv1rM8Y99y5YijTgkIjuKlBym/CXStDTEZc2DnXODnHNDqph2FLCH/3ch8DiAmWUDtwL7AvsAt5pZs90J4smJi1i2aStFpeU7vF5UWs6yTVt5cuKi3Vm9iDReYc9hyl8iTUu4z4WPAZ533hhSP5hZlpm1BQ4CPnPO5QCY2WfAkcAru/pGL36/bKfEVqGotJynvl5C+6xkUhLiyEyOJyM5nkz/LyMpjrhYNc8TkZ00SA4LzF9d81YzbM0sJrfuzcKsDhSVwlNfL6FtZjKJcTFkJMWTmfJb/spMjicpPrY+PquINJBQF2cO+J+ZOeBJ59y/Kk1vD6wIeL7Sf62613dgZhfiHa3SqVOnGgPJ3VpS4/StxWVc++bMaqdnpybQJiOJtplJtMlMol1WMm0ykujcPIXuLdNolppQ4/pFJCpFRA4LzF99Ni3hjHmfcda8/7E5IZUprXoxuXVv7nkln/zE1CqXT46P/TV3tcmsyGPJdGiWTI+WabTLSiY2xqrfCiLSoEJdnO3vnFtlZq2Az8xsnnNuYn2t3E+U/wIYMmRIjSO4N0uJJ6eGAi0rJZ4PLxvJlqJS8raVkL+thDz/b/PWEjYUFrE2bzur87bz8/LcnYq97NQEurVIpXvLNLq3SmXPthn0a5epok0kukVEDgvMXx92G8E37Qey9/oFDFk3jyHr53Poyp/BDNuzLyVD9qNg4FA2tetGXlEZedtK2FRYzLr87azJ28YPizaxrqCIsvLf3i4hLoZuLVLp1tLLYT1bp9OvfSads1OIUdEm0uBCWpw551b5/643s7fx2l4EJrZVQMeA5x3811bhXRYIfH3C7sRy5rDOPPnV4iovbSbGxXD2sM60z0oOen3bS8pYvXkbSzdtYfGGLSzasIVFGwr5Yt56Xptc9Ot87bOS6d8+k37tM+jfIYu9OmWRkRS/Ox9FRBpIpOSwyvkrLzGNLzvuzZcd9yY5Fq7qUs7oouUUTpyIe+E/NHv+KVpkZ5O6/wjSRh5A6v4jiGv2W5O3snLHhoIiVuRuZdH6QhZv3MKi9YXMXVPAp7PX/Vq4pSfG0addBv3bZ9K/QyaDOmbRKTsFMxVsIqFkXlOJEKzYLBWIcc4V+I8/A+5wzn0SMM/RwB+BUXgNZx9yzu3jN6adAlT0fPoZGFzRfqMqQ4YMcTWNS1fR26lyo9rEuBg6N0/h7UtG1Ft39M1bi5m9Op9Zq/KYuSqPWavyWLppKwBm0KdtBkO7ZLNv12yGds2mRVpivbyvSFNjZlOqaahfH+uOmBxWl/xVmpPDlm+/pXDi12z55hvKcnPBjOQBA0g9YCRpBxxAUt++WEzV7WiLS8tZuL4gIH/lM3dN/q/v2zojkX26Nmefrtns0yWbPVql6eyayC6oKX+FsjjrBrztP40DXnbO3WVmFwM4554w7/DrEbyGsluB85xzk/3lzwdu8Je/yzn3TE3vV1txBgH3CfphOblbi2mWksCZ+3VqkPsE5W8vYebKPH5amsOPS3L4eXku20u8ZNe9ZSoH9GzJQb1asW/XbDXeFQlSiIuziMphu5K/XFkZ22fPpnDi1xROnMj2mTPBOWJrOKtWlZKycn5ZX8iUZbn8uCSHSUs2sS7fu0KQlRLP8O7NOahnKw7s1ZLWGUk1rktEPGEpzhpaMMVZJCkuLWfmKq9Y+27RJn5YvIni0nKS4mMY1q05B/duxcG9WtExOyXcoYpErFAWZw2tIXLY7p5Vq+CcY0XONiYt2cSkJTlMXLCB9QVesdanbQYH9fIONvfulKWe7iLVUHEWBbYVl/HD4k1MmL+eCQs2sMy/DNq3XQaj+rdlVP+2dG1RdU8skaZKxdmuq6+zauAVa3PXFDBhwXomzN/AlGW5lJU7mqcmcES/Nozq15b9umWrUBMJoOIsCi3ZuIXP56zj41lr+Hn5ZgB6t0nn6P5tGTWgLd1bpoU3QJEIoOKs/tTXWTXwmnF8vWAjH89aw/h569laXEazlHiO6NuGo/q3ZUT35irUpMlTcRblVm/exiez1vLRzDVMXpYLwMCOWZw4uAOjB7QjM0W9P6VpUnEWGq6sjO2zZnln1b7+erfOqm0vKWPC/A18PGsNn89Zx5biMlqkJXL8Xu04cXBHerVJb4BPJBJ5VJw1ImvztvP+9NW8+fNK5q0tICE2hsP7tObEwR0YuUcLHY1Kk6LirGFUeVYtJobk/v3rdFatolB76+eVjJ+3ntJyR//2mZywd3tGD2pPtu4LKU2IirNGyDnH7NX5vDFlJe9OW0Xu1hJaZyRy2j6dOH2fTrRSjylpAlScNbxaz6odcCCpI4bXelZtU2ER701fzRtTVjJ7dT7xscZR/dpy9rDODO7cTPdSk0ZPxVkjV1xazvh563nlx+V8tWADcTHGEf3acPZ+ndmna7aSnDRaKs7Crz7Oqs1dk89rP63gzZ9XUrC9lD3bZnD2sM6MGdSOlIRwDwEtEhoqzpqQpRu38NKkZbw+eSV520ro1Tqdc4Z3Yeze7XX/NGl0VJxFlt09q7a1uJR3pq7m+e+XMm9tAelJcZw0uCPnjeii2wpJo6PirAnaVlzG+9NX8/wPS5m1Kp8WaQmcO7wLZ+7XmawUteuQxkHFWWTb1bNqzjmmLMvl+e+X8fGsNZQ7OLp/Wy48oBv92meG6dOI1C8VZ02Yc45JS3J44qtFTJi/gZSEWE4d2okLRnat01iiIpFIxVn0qOmsWtrI/UkdeUCVZ9XW5m3n6W+X8PKk5RQWlTJyjxZcfGB3hndvriYbEtVUnAngtev498TFvDd9NQBjBrXnT4f0oItubitRSsVZ9Pr1rNpXE72zaps313hWLW9bCS9PWs7T3y5hQ0ER/dpncNkhe3B4n9Yq0iQqqTiTHazavI2nvl7My5OWU1ruGLtXe/50yB50aq42HRJdVJw1DnU5q7a9pIx3pq7i8a8WsWzTVvq3z+Qvh+3BIb1bqUiTqKLiTKq0vmA7T0xYzIuTllFe7jhxcAcuPbiHGt5K1FBx1jgFc1YtrveevDN9DQ+NX8iKnG0M7JDJXw7vyUE9W6pIk6ig4kxqtC5/O49PWMTLk5ZT7hyn7tORPx/ak5bpiVBUCN89BD89BVtzICUbhv4ehl8GiRpCSsJLxVnjV9tZtaQR+zM+rSv/+HE9qzZvY69OWVx/1J7s0zVb+UsimoozCcqavG08+uUvvPrjChLjYvjj/m25aMGFxGxeCqXbf5sxLgmadYXff64EJ2Gl4qzpqe6sWmL//iztNoB/FbVmUnwrjt4ziwfyryQxf5nyl0QkFWdSJ4s3FHLfJ/PYc96jXBz/PkmU7DxTXBKM+DMcfEPDByjiU3HWtFV3Vq04LZOc5gn0b7uIrLZbiEus9P+c8pdEABVnsktK7ulKfFEOABtmpZHaqpiUVsW/zZDSAq5ZFKboRFScyY52OKv22fuUFRmYIzm7hNS220lrV0RSsxLMUP6SsKspf2lcDKlWfFEuAGUlRv6yZDbOSafd0M1kdt3mzbBtUxijExHZUVx2NpnHHkvmscfibvk323PjKFydROGaRDbOSmfjrAxiE8tIa1tEatstpObm1joGqEg4qDiT6qVkw9ZNxMY7uhy+kZXfZLN6UjOKC+No0a8AUrJRnygRiUSWlk1yzCaSm5fQsn8Bpdtj2LI2kcI1iRSuTiRvaQqM2J+k/v1JO/AA0kYeQFLfPjWOASrSUPQtlOoN/b3XNgOITXB0OnATmV23snF2Oit+yObZbQfz2Zx1NJZL4yLSiATkL4C4pHIyu2yj/bDNdByTw6SD+zJl2FEUFhSw8eFHWHrSSSzYfySrr72WvA8+pDQ3N4zBS1OnM2dSveGXwZz3IHcJlG7HYqHtPpuJy3Bsmp5K6+1ruHLzVwwd2JU7j+tHOw0HJSKRolL+qlAWk0BZalu+b3Mun64wWnQ6lEvPNvZaN5/Y6TPYPP5L8t59D2JidFZNwkYdAqRm/n2C3I9PwbYcShMy2NB5NBs2dSf+38+wLasFVw8+l3UZrbjmyN6cuV9nYmN0sVMahjoESI1+vc/Zf3BbN/2av9b1OJXyuGTmbSzm8Sn5rMwvZUTHJM4flE5WghGzeDGx06cTO30GMYsXg3PEZGeTXsMYoCJ1pd6astu2bdvGnDlzdngtZsECEh/8B2bGC6Mu5eXtzdmrUxb3jh1ArzbpYYpUmhIVZxKsuXPnsnXr1p1eLyl3vD1vC2/OLSQp1jhnYDoHd0n+bZSB/HxiZ8wgdvoMYmfOxAoKvLNqAwaQdsBInVWTXabiTHZbVcUZgK1bR+L9DxCzfj2rL7ySqza3o7ColIsP7M4fD+lBYlxsGKKVpkLFmQSruuKswsr8Up6YksfcjSX0b5XARYMzaJtWqeVPeflvZ9WmTSdmyRKdVZNdpuJMdlt1xRkAW7aQOG4csXPmkvKHixjX/kDenraaXq3TeeDkgfRrn9mwwUqToeJMglVbcQZQ7hyfL97GCzMKKHNw1oA0juieQkx1Y3XqrJrsBhVnsttqLM4ASktJ+M9/iJv4NRnHHsOCs//Mte/NI2dLMZcdugf/d1B34mOVoKR+qTiTYAVTnFXYtK2MxyfnMXVtMQNaJXDJ0ExaptRyFaDirNq06cRO/+2sWsUYoDqrJpWpOJPdVmtxBuAcce+9R8Lr/yV58GAy/v4gt321ivemr2ZAh0wePHkgPVqpLZrUHxVnEqy6FGcAzjk+W7KN56YVEGNw3l7pHNw5oC1abXRWTWqh4kx2W1DFmS/2++9JfPJfxLdtS6d/Pcln+Ync9M5MthSXcc0RvTh/RFdi1KNT6oGKMwlWXYuzCmsLS3nkJ68t2tB2iVw8OIOspDq2pa3xrNpIUg8YSdqIEcRmZdU5PoleKs5kt9WlOIPfenLGxsTQ8ZGH2dJ7ANe/NZPP565jRI/mPHjyIFpnJNW+IpEaqDiTYO1qcQZeW7QPF27lpZkFJMfH8KehmezdNnHXg9FZNUHFmdSDuhZnsGNPznZ330XGscfy2k8ruP39OSQnxPL3Ewdw6J6tQxSxNAUqziRYu1OcVViRX8qDP2xmeV4px+yRwpn904mP3c2rADqr1mSpOJPdtivFGbBDT84Wl15Kiz9eyqINW/jTK1OZuyafc4d34bqjepMUr1tuSN2pOJNg1UdxBlBc5nh+RgEf/7KVrllxXL5fFu3T63GwncCzajNmYIWFOqvWSKk4k922y8UZVOrJeSxt7/orxRbLfZ/M45lvl7Jn2wwePm2QOgtInak4k2DVV3FW4afV23n0pzyKy+D3e1W6cW19KS8nZlHFaAU6q9bYqDiT3bZbxRl4PTnffY+E/3o9OTs88jBxzZoxft46rvrvDLYWl3LXcf05YXCH+gtaGj0VZxKs+i7OwLvlxkOT8pi1oZiRnZK4aHAGyXEhPKOls2qNiooz2W27XZz5fu3J2a4dnZ58goQuXVifv53LXp3KD4tzOG2fTtx6bB9d5pSgqDiTYIWiOAMoc463527htdmFtMuI4+phWXTIqMfLnNXRWbWop+JMdlt9FWewc0/OlKFDKS0r54HPFvD4hEX0b5/JY2fsTcfslHp5P2m8VJxJsEJVnFWYsa6If0zKo7jUccnQDEZ0TA7Ze1WptrNqBxxIUp89dVYtgqg4k91Wn8UZ7NyTM3P0aAA+m7OOK16fRowZ/zhlIIf0Vm9OqZ6KMwlWqIsz8C5zPvD9ZuZvKmFUjxTOHphOfDju6aizalEhrMWZmcUCk4FVzrljKk3rDDwNtARygDOdcyv9aX8DjgZigM+AP7saglViC636Ls6AKntymhnLNm3h/178mTlr8rn04O5ccXgvYnXTWqlCqIuzhspfoBwWag1RnAGUljtemFHABwu3skd2PFcNy6JFbUM/hVrFWbVp0737qumsWkSoKX81xJ74MzC3mmn3A8875wYAdwD3AJjZcGAEMADoBwwFDgx9qNKgUlMpuvZaSkeOZOOjj7L6mmspLy6mc/NU3rpkOKcM6cijXy7igud+Im9bSbijlaZJ+UvqJC7GOG9QBlful8WK/FKu+XwT8zYWhzeojAzK9t+f4j9eyrbHH2P7bbdRMmYMW/Pz2fjwIyw98UQW7j+S1ddeR96HH1K2eXN445XQFmdm1gHv6PGpambpA4z3H38JjPEfOyAJSAASgXhgXegilbCJi6P4ogspPukk8t9/n+XnnU9pbi5J8bHcd+IA/npcP75ZuJHjH/2WX9YXhjtaaUKUv2R3DO+YxL2HNic53rh1Qg6fLw79WbugxMRQvkcPSk48ge133sHWRx+h6OKLKerVi81ffMHqK69iwfARLDn1NDb++9/hjrbJCvWZs3HANUB5NdOnA2P9x8cD6WbW3Dn3PV6yW+P/feqc2+no1cwuNLPJZjZ5w4YN9R68NBAzSo8bQ9EfL2XbjBksPfU0ipcuBeDM/Trz8h/2I29bCcc/+i1fzlsf3lilKRlHCPMXKIc1dh0z4rjv0Ob0bZXA41Py+ffP+ZSWR1g778xMykZWfVat8Jtvwx1dkxWy4szMjgHWO+em1DDbVcCBZjYV77T/KqDMzHoAewIdgPbAIWY2svLCzrl/OeeGOOeGtGzZsv4/hDSosmHD2H799RTn5rDklFPZ6re/2adrNu/9aX86NU/h/Od+4vEJi2gsHVkkMjVE/gLlsKYgLSGGG/dvxuieKXyyaCt3TMwlv6i6ej/MKp1Va//E4+GOqMkK5ZmzEcBoM1sKvIqXoF4MnME5t9o5N9Y5txdwo//aZryj0B+cc4XOuULgY2BYCGOVCFHeqyfbb7uN0tRUlp13PnnvvQdA+6xk3rh4OMcMaMd9n8zjslensb2kLMzRSiOm/CX1JjbGOGdgBpftk8mCTcVc+/kmluVFfjtai48PdwhNVsiKM+fc9c65Ds65LsCpwHjn3JmB85hZCzOriOF6vJ5PAMvxjkjjzCwe76i0uka50si41q3ZfustlO3Rg9XXXMuGhx/BOUdyQiwPnTqIa47sxQczVnPav39gY2FRuMOVRkj5S0LhwM7J/PXg5pSWO24cn8PUtcpfUrUG7zdrZneY2Wj/6UHAfDNbALQG7vJffwNYBMzEa9cx3Tn3fkPHKmGUllZlT04z45KDevD4GYOZuyaf4x/7ll/WF4Q7WmkilL9kd/XIjufeQ5vTJi2Wu7/J5ZNFEdJRQCKKbkIrQQnJfc6CUc2YnADTV2zmgucmU1RaxhNnDmZEjxYNH5+ElW5CK8FqqPucBWtbaTn/+CGPKWuKOLZnCmcNSCe2vgdO300DBw4kLq4BhqJqosJ9nzORXVfRk/PSS3bqyTmwYxbvXDqctplJnPP0j7z+04rwxioiEqTkuBiuHZHFqB4pvL9gK/d/t5ntpRHaUUAanIoziQplw4d7PTlzduzJ2aFZCm/833CGdW/ONW/O4G+fzFNPThGJCrFmXLBXBucPSmfy6iJumZBD7nZ1dBIVZxJFvJ6ct1KakrJDT86MpHiePncop+3TiccmLOLK/06npExHoCISHY7eI5VrR2SxMr+MG8bnsLqgNNwhSZipOJOo4tq0Yfttt/7Wk/ORR3HOER8bw93H9+OKw3vy1s+r+P1zk9lSpAQnItFhSLskbj+oGdtLyrlx/CZ+yYn8W21I6Kg4k+gT2JPzkUd26Ml52aF7cM/Y/ny9cAOn//sHNulWGyISJfbITuCuQ5qTFBfDrRN0q42mTMWZRKeKMTlPPHGHMTkBTtunE0+cOZh5aws48YnvWZETOT20RERq0i49jrsPyaZNWiz3fJPLV8u2hTskCQMVZxK9zCg9/jivJ+f06Tv05Pxd3za89Pt9ydlSzNjHv2PO6vzwxioiEqRmybHccXA2e7ZI4KEf83hv/pZwhyQNTMWZRL2y4cPZfsMNO/XkHNIlm/9ePIy4GOOUf33P5KU5YY5URCQ4qfEx3DSyGcM7JPHcjAJenFmgnuhNiIozaRR26sn5vndD9p6t03nz/4bTMi2Rs/7zI18v3BDmSEVEghMfa1y+Xya/65bM2/O28NTUAspVoDUJKs6k0fi1J2ePHqy++ppfe3K2y0rmtYuG0aVFKhc8O5lPZ68Nd6giIkGJMePCvTM4rlcqnyzayiM/5VFWrgKtsVNxJo1LWhpF111L6cj9d+jJ2TI9kVf/sB9922dwyUs/8/bUleGOVEQkKGbGmf3TOL1fGl8t284DP2ympEwFWmOm4kwan7g4ii+66LeenOdfQGluLpkp8bx4wb7s2zWbK16fzos/LAt3pCIiQTEzTtgzjfMHpTNpVRH3fJur4Z4aMRVn0jgF9uScNu3XnpypiXE8fe5QDunVipvemcUTXy0Kd6QiIkE7eo9ULh2Swcx1xdw5MZctJSrQGiMVZ9KoVdWTMyk+lifOGsyxA9tx78fzeGT8wnCHKSIStEO6pnD5flkszClRgdZIqTiTRq+qnpzxsTH84+SBHDeoHff/bwEPfaECTUSix/COSVw5LIsluSXc8VUuW4pVoDUmKs6kSdipJ+ejjxIbYzxw8iDG7t2eBz9bwLjPF4Q7TBGRoO3bPomrhmexdHMJt0/MoVAFWqOh4kyajsCenA8/wuprr8NKS/j7iQM5cXAHxn2+kAf/N183ehSRqDG0XRLXjMhiWV4pt3+VQ4EKtEZBxZk0Lb/25DyB/PfeY/n5F+DyNvO3EwZwypCOPDT+F+5XgSYiUWRw2ySuHd6MFfl+gVakAi3aqTiTpseM0uOP36EnZ+mK5dwztj+n7dOJR79cxN8+VYEmItFj77aJXDuiGavyS7lVBVrUU3EmTZbXk/P6X3tybp/6M3cd14/T9+3E4xMW8U91EhCRKLJXm0Su278ZqwtKuWNijjoJRDEVZ9Kklffq5fXkTE5m2bnnUfDhB/x1TL9f26A9PkH3QROR6DGwdSJXD2/G8rxS/vp1Ltt0m42opOJMmrzKPTk3Pf4Y947tz+iB7bjvk3k8/c2ScIcoIhK0wW0TuWJYFr/klnDXNxpJIBqpOBOBnXpyrrv+eu4/bk+O7NuGOz6Yo6GeRCSq7Ns+ib/sm8n8jSXc++1mijQWZ1SJC3cAIhHD78lZ3ro1+W+8Scnq1fzjn/+kuKycm96ZRWJcDCcN6RjuKEVEgjKiYzIlZfDIT3nc/91mrhmRRXyMhTssCYLOnIkEqujJeYnXk3PV6afzz5EtGLlHC659cwbvTV8d7ghFRIJ2UJdkLhycwc9ri/jHD5spK9cZtGig4kykCmUjfuvJufqM03mobwxDumRzxWvT+HLe+nCHJyIStN91S+H8QelMWlXEY5PzKNdtgiKeijORagT25Fz7hwt4qOU6erdN5+IXp/DT0pxwhyciErSj90jllL5pTFi2neemF+g+jhFOxZlIDX7tydm9O7k3XM9jpVNpn5XE+c/+xOzVeeEOT0QkaCftmcrRe6TwwcKtvDl3S7jDkRqoOBOpTUVPzv33Z8uTj/Pkyg/JinOc8/SPLNmoBCci0cHMOHdgOgd2TuKV2YV8/IvyV6RScSYSjPh4ii/2xuQs/eRD/jXzeZK3b+HMpyaxJm9buKMTEQlKjBmXDMlkaLtEnppawNfLlb8ikYozkWAF9ORk9kye+P5Rktev5qz//EjOluJwRyciEpS4GOOK/bLo2zKeh3/MY8qa7eEOSSpRcSZSRxU9OS0/j39+/TCpC2Zx/rM/sa24LNyhiYgEJSHWuG5EM7pkxXH/95uZv0kHmJFExZnILqjoyUlqCnd/8wTNvh/Pn175mdIyDZMiItEhJT6GG/dvRnZSLPd8k8uqgtJwhyQ+FWciu6iiJ6fr0YNrJr9My7df4uZ3ZqqLuohEjcykWG4+oBlmxl8n5pK7XVcAIoGGbxLZHenpFF13LQlP/Yezv/mUz5/YyD+SLiex3Q+8Nv81NhdtJisxi1N6ncJ5/c4jJT4l3BGLiOygTVocN+7fjFsn5PDXb9ayb59JTMwbz5ZZW5S/wiTkZ87MLNbMpprZB1VM62xmX5jZDDObYGYdAqZ1MrP/mdlcM5tjZl1CHavILqnoyXnCWA5bMYVOD57PKz/+h9yiXByO3KJcnpn9DKd/dDpbS7aGO1qpA+UvaSp6ZMfzp/2S2JD1EB9v+ojCskLlrzBqiMuafwbmVjPtfuB559wA4A7gnoBpzwN/d87tCewDaMwciVxmlI4dy9dnDKXn2mJuf24brXN+u7xZVFbEyoKVPDPrmTAGKbtA+UuajDWx44lPysHZjm3PlL8aXkiLM/9I8mjgqWpm6QOM9x9/CYzxl+sDxDnnPgNwzhU651SyS8R7putC7jgtlrRtcNfzZfRasWOB9tr818IYndSF8pc0NV/kfEEZJVVOU/5qWKE+czYOuAaorgvbdGCs//h4IN3MmgM9gc1m9pZ/SeHvZhZbeWEzu9DMJpvZ5A0bNoQgfJG6KSwrZH5H48azYylIhuYFO3YO2Fy0OTyBya4YRwjzFyiHSWQpLCuscbryV8MJWXFmZscA651zU2qY7SrgQDObChwIrALK8DoqjPSnDwW6AedWXtg59y/n3BDn3JCWLVvW8ycQqbu02DQA1mUbV18Qy3d9dvyJZSVmhSEqqauGyF+gHCaRpSJ/VUf5q+GE8szZCGC0mS0FXgUOMbMXA2dwzq12zo11zu0F3Oi/thlYCUxzzi12zpUC7wB7hzBWkXpxaPahxFs8AKVxtsO0hNhETul1SjjCkrpT/pImJzB/VRYfk6D81YBCVpw55653znVwznUBTgXGO+fODJzHzFqYWUUM1wNP+49/ArLMrOJQ8hBgTqhiFakvR7U4ilYJrXZKcK48joTylpzb99zwBCZ1ovwlTVF1+YvyOMqKsjmmy2nhCawJavCb0JrZHWY22n96EDDfzBYArYG7AJxzZXiXBL4ws5mAAf9u6FhF6iopNombu93MqBajyIjLwDCaJTZjYPrxrJn7B174fm24Q5TdoPwljVlg/kqPTf81f43tfhbFy//IZS/N0TB1DcQay93MhwwZ4iZPnhzuMBqtbdu2MWeODv7rIjMzkx49egBQXu7482vTeH/6ah47Y29G9W8b5ugaBzOb4pwbEu446oNyWGjNnTuXrVvVabYuBg4cSFycd6/6z+as48IXJvO7Pq15/IzBxMRYLUtLbWrKXxq+SaQBxMQYfz9xAIM7N+Py16YxdXluuEMSEQna4X1ac9PRffh09jru/WReuMNp9IIqzsyslZkdb2aXmtn5ZrZPQFsLEQlCUnws/zprMK0yErnwhSmsydsW7pBERIJ2/ogunLVfZ/41cTFvTFkZ7nAatRoLLDM72Mw+BT4EjgLa4t148SZgppndbmYZoQ9TpHFonpbIf84ZytaiUi58forab4hI1DAzbjm2DyN6NOeGt2YyZVlOuENqtGo7+zUK+INzbqhz7kLn3E3Ouaucc6OBgcBU4PCQRynSiPRsnc5Dp+3FrNV5XPXGdBpLu08RafziY2N49PS9aZeVxEUvTGHVZl0BCIUaizPn3NXOueVm1rWKyR2dc+84594MUWwijdahe7bmuiN78+GMNTw8/pdwhyMiErSslASeOmcoRaXl/P65yWwpKq19IamTYNuNVVWAvVGfgYg0NRce0I2xe7fnwc8W8PHMNeEOR0QkaD1apfHwaXsxf20+V74+nfJyXQGoT7W1OettZicAmWY2NuDvXCCpQSIUaaTMjLuP789enbK44vXpzFqVF+6QRESCdlCvVtwwak8+mb2WcZ8vCHc4jUptZ856AccAWcCxAX97A38IaWQiTUBSfCxPnjWYrJR4Lnx+MhsKisIdkohI0C7YvysnD+nAQ+N/4f3pq8MdTqNRW5uzd51z5wHHOOfOC/i7zDn3XQPFKNKotUpP4t9nDyFnazGXvvQzJWXl4Q5JRCQoZsadx/VjSOdmXP3GdOaszg93SI1CbZc1bzKzZs6576uZfoiZHROa0ESajn7tM7nvhAH8uDSHuz6cG+5wRESClhgXy2Nn7k1mcjwXvTiZzVuLwx1S1KvtsuZM4AMz+8LM/m5m15jZLWb2gj9m3LHApNCHKdL4jRnUnt/v35Vnv1vKfyevCHc4IiJBa5WexBNnDmZdXhF/emUqZeogsFtqK85OdM6NAD4FZgOxQD7wIrCPc+5y59yGEMco0mRcd1Rvhndvzo3vzGLGys3hDkdEJGh7dWrGncf15euFG/n7p/PDHU5Uq604G2xm7YAzgPeAJ4HngZ+A5BDHJtLkxMXG8Mjpe9MyLZGLX5jCxkJ1EBCR6HHK0E6csW8nnvhqER/MUAeBXVVbcfYE8AXQG5gc8DfF/1dE6ll2agJPnjWYTVvUQUBEos+tx/ZlcOdmXP3fGcxbqw4Cu6K23poPOef2BJ52znUL+OvqnOvWQDGKNDn92mdyz9j+TFqSw90fqYOAiESPhLgYHj9jb9KT4rjohSnkbS0Jd0hRJ6gRApxz/xfqQERkR2P37sB5I7rwzLdLeXfaqnCHIyIStFYZSTx+5t6s3ryNv7w2VSMI1FGwwzeJSBjcMGpPhnZpxnVvzmTBuoJwhyMiErTBnbO55Zg+fDl/A49N0BjCdaHiTCSCxfsdBFIT47j4xSkUaoBhEYkiZ+7XmTGD2vHAZwv4ZuHGcIcTNVSciUS41hlJPHzaXizduIVr35yBc7o8ICLRwcy4Z2x/9miVxmWvTmVN3rZwhxQVVJyJRIFh3Ztz9RG9+XDGGp75dmm4wxERCVpKQhyPnzmYopIyLnnpZ4pL1QO9NirORKLExQd24/A+rbn7o7lMWZYT7nBERILWvWUafztxIFOXb1YP9CCoOBOJEmbG/ScNpH2zZC556WfdoFZEosrRA9py/ghviLr3pusGtTVRcSYSRTKT43nsjL3ZvLWEyzR+nYhEmetH9WZw52Zc9+YMflmvHujVUXEmEmX6tsvkzjH9+G7RJh76YmG4wxERCVp8bAyPnr43yfGxXPLSz2wrLgt3SBFJxZlIFDp5aEfG7t2eh8Yv5Ltf1D1dRKJHm8wk/nHKIBauL+T292eHO5yIpOJMJErdOaYf3Vqk8ufXprGhQO3PRCR6HNCzJZcc1J1Xf1qhEVCqoOJMJEqlJsbx6Bl7k7+thCten6bhUUQkqlx+WE/26ZLNDW/NZPGGwnCHE1FUnIlEsd5tMrh9dF++XrhRw6OISFSJi43hn6cNIiEuhktfnsr2ErU/q6DiTCTKnTK0I6MHtuPBzxYwafGmcIcjIhK0tpnJPHjyIOauyefOD+aEO5yIoeJMJMqZGXeP7U/n5qlc9upUNun+ZyISRQ7u3YqLDujGS5OW88EM3f8MVJyJNAppiXE8cvpe5G4t4YrXp6v9mYhElauO6MXenbK47s2ZLN24JdzhhJ2KM5FGom+7TG4+pg9fLdjA098uCXc4IiJBi4+N4eHT9yY2xvjzq1Ob/PibKs5EGpEz9+3E7/q05r5P5jFrVV64wxERCVr7rGTuO6E/01fm8eBnC8IdTlipOBNpRMyM+04YQPPURC57ZSpbi0vDHZKISNCO7NeW0/bpxJMTF/FtE77BtoozkUamWWoCD54ykCWbtnD7e+r9JCLR5ZZj+tC9ZRqXvzaNnC3F4Q4nLEJenJlZrJlNNbMPqpjW2cy+MLMZZjbBzDpUmp5hZivN7JFQxynSmAzv3oJLDurOa5NX8OGMNeEOJ2opf4k0vOSEWB46dS82by3hmjem41zT6+DUEGfO/gzMrWba/cDzzrkBwB3APZWm3wlMDGFsIo3WXw7ryaCOWVz31gxW5m4NdzjRSvlLJAz6tMvg+lG9+Xzuel74YVm4w2lwIS3O/CPJo4GnqpmlDzDef/wlMCZg2cFAa+B/oYxRpLGKj43hoVP3wjn4y6vTKC1r2r2f6kr5SyS8zh3ehYN7teSvH85l3tr8cIfToEJ95mwccA1Q3f8K04Gx/uPjgXQza25mMcADwFU1rdzMLjSzyWY2ecOGDfUUskjj0al5Cncd34/Jy3J5eLyGd6qjcYQwf4FymEhNzIy/nzSQjKR4LnulaQ3vFLLizMyOAdY756bUMNtVwIFmNhU4EFgFlAGXAB8551bW9B7OuX8554Y454a0bNmyvkIXaVTGDGrP2L3a8/D4hUxZlhvucKJCQ+QvUA4TqU2LtET+ccpAFqwr5N6P54U7nAYTF8J1jwBGm9koIAnIMLMXnXNnVszgnFuNf+RpZmnACc65zWY2DBhpZpcAaUCCmRU6564LYbwijdbtY/ry49Icrnh9Gh9eNpK0xFD+9BsF5S+RCDFyj5acP6IrT3+7hIN7t+LAno3/QCZkZ86cc9c75zo457oApwLjAxMbgJm18C8BAFwPPO0ve4ZzrpO/7FV4jW6V2ER2UXpSPA+ePIjlOVu5833dXqM2yl8ikeWaI3vRs3UaV/93OrlN4PYaDX6fMzO7w8xG+08PAuab2QK8xrN3NXQ8Ik3FPl2zufhA7/Yan85eG+5wopLyl0h4JMXHMu6UvcjdWswNb89s9LfXaJBrG865CcAE//EtAa+/AbxRy7LPAs+GLDiRJuTyw3oyccEGrn9rJnt1yqJVelK4Q4p4yl8ikaFPuwyu/F0v7v14Hm/+vIoTB3eofaEopRECRJqQhLgYxp0yiC1FpVz7xoxGf/QpIo3LH0Z2Y5+u2dz23mxW5DTe+zeqOBNpYvZonc71R/Xmy/kbeHHS8nCHIyIStNgY48GTB2LA5a9No6y8cR5gqjgTaYLOHtaFkXu04K4P57BoQ2G4wxERCVqHZincPqYvk5fl8sRXi8IdTkioOBNpgmJijPtPGkhSfCyXv6bRA0Qkuhy/V3uO7t+Wf3y2gFmr8sIdTr1TcSbSRLXOSOKu4/ozY2Uej09onEefItI4mRl3Hd+P7NQErnx9OkWljWv0ABVnIk3Y0QPacsyAtvzzi4XMXt34jj5FpPHKSkng3hP6M39dAeM+XxjucOqVijORJu7OMf3ISmmcR58i0rgd0rs1pwzpyJNfLeLn5Y1neDoVZyJNXLPUBO4d2595awt46IvGdfQpIo3fTcfsSdvMZK56fTrbihvHAaaKMxHhsD6tOWlwBx6fsIipjejoU0Qav/SkeP5+4gAWb9zC3z5tHIOjqzgTEQBuPrYPbTKSuLIRHX2KSNMwvEcLzhnWmWe+Xcp3izaGO5zdpuJMRADISIrnbycOZPHGLfz90/nhDkdEpE6uPao3XZqncPV/Z1BYVBrucHaLijMR+dX+e7TgrP068/S3S/h+0aZwhyMiErSUhDgeOHkga/K2cdeHc8Idzm5RcSYiO7h+VG86N0/hmjensyXKjz5FpGkZ3DmbPxzQjVd+XMFXCzaEO5xdpuJMRHaQkhDH304YwIqcbbq8KSJR5/LDetK9ZSrXvzmDgu0l4Q5nl6g4E5Gd7NutOecM68xz3y/lp6U54Q5HRCRoSfGx/O3EgazJ3869H0dn700VZyJSpWuO7E37rGSueWOGem+KSFQZ3LkZF4zoykuTlvPdL9HXe1PFmYhUKTUxjvtOGMCSjVt48DNd3hSR6HLl73rRpXkK1741I+raz6o4E5FqjejRgtP26cR/vlnSqIZGEZHGLznBu7y5Mjf62s+qOBORGt0wqjdtMpK45o0ZbC/R5U0RiR77dM3mnGFdePa7pfy4JHraz6o4E5EapSfFc88JA/hlfaHG3hSRqHP1Eb3omJ3MtW9GT/tZFWciUqsDe7bkpMEdeHLiYmauzAt3OCIiQUtNjOO+sdHVflbFmYgE5aZj+tAiLYGr35hOSVl5uMMREQna8B4tOH1fr/3stBWbwx1OrVSciUhQMpPjuXNMP+atLeBfExeHOxwRkTq57qjetEpP4ro3Z0T8AaaKMxEJ2u/6tmFU/zb884uFLN5QGO5wRESClpEUz53HRccBpoozEamT20b3JSkuhuvfmkl5uQt3OCIiQTu8T2uO7t+Wf36xkEURfICp4kxE6qRVehI3Hr0nk5bk8NrkFeEOR0SkTm4d3SfiDzBVnIlInZ08pCPDujXn7o/msi5/e7jDEREJWqv0JG46ug8/Lsnh1Z8i8wBTxZmI1JmZcffY/hSXlnPru7PDHY6ISJ2cNKQDw7s3554IPcBUcSYiu6Rri1T+clhPPpm9lk9mrQl3OCIiQTMz7j6+P8Vl5dzy7qxwh7MTFWcisst+P7IrfdpmcMu7s8nbVhLucEREgtalRSqXH96TT2evi7gDTBVnIrLL4mNjuO+EAWwsLOK+T+aFOxwRkTr5/f5d6dsug5vfnU3+9sg5wFRxJiK7pX+HTM4f0ZWXJy1nyrLccIcjIhK0uNgY7h07gE2FRdz/aeQM7aTiTER22+WH96RdZhI3vDUz4u+8LSISqH+HTM4Z3oUXflgWMUM7qTgTkd2WmhjH7WP6MX9dAf/5Zkm4wxERqZMrf9eL1ulJXP/WTEoj4ABTxZmI1IvD+7TmiL6tGff5AlbkbA13OCIiQUtLjOO20X2ZuyafZ75dGu5wQl+cmVmsmU01sw+qmNbZzL4wsxlmNsHMOvivDzKz781stj/tlFDHKSK777bRfYk14+Z3Z+FcZN55uy6Uv0SajiP6tuawPVvx4GcLWJkb3gPMhjhz9mdgbjXT7geed84NAO4A7vFf3wqc7ZzrCxwJjDOzrFAHKiK7p21mMlf+rhcT5m/go5lrwx1OfVD+EmkizIzbx/TDDG59d3ZYDzBDWpz5R5JHA09VM0sfYLz/+EtgDIBzboFzbqH/eDWwHmgZylhFpH6cM7wL/dpncNv7kdU1va6Uv0SanvZZyVx+WE++mLeeT2eH7wAz1GfOxgHXANW1rpsOjPUfHw+km1nzwBnMbB8gAVhUeWEzu9DMJpvZ5A0bNtRb0CKy62JjjHuOj7yu6btgHCHMX/505TCRCHPeiC7s2TaD296bQ0GYDjBDVpyZ2THAeufclBpmuwo40MymAgcCq4CygHW0BV4AznPO7ZQgnXP/cs4Ncc4NadlSB6YikaJ/h0zOHuZ1TZ+6PPrufdYQ+QuUw0QiUVxsDPeM7c+6gu088L8FYYkhlGfORgCjzWwp8CpwiJm9GDiDc261c26sc24v4Eb/tc0AZpYBfAjc6Jz7IYRxikgIXPm7nrRKT+Smd2ZRVh51nQOUv0SasEEdszhz3848//1SZq3Ka/D3D1lx5py73jnXwTnXBTgVGO+cOzNwHjNrYWYVMVwPPO2/ngC8jdfY9o1QxSgioZOeFM/Nx/Rh9up8XvxhWbjDqRPlLxG56oheZKcmcNM7syhv4APMBr/PmZndYWaj/acHAfPNbAHQGrjLf/1k4ADgXDOb5v8NauhYRWT3HN2/Lfv3aMH9/5vPhoKicIez25S/RJqOzOR4rj9qT6at2Mzrk1c06HtbY7gXEcCQIUPc5MmTwx1Go7Vt2zbmzJkT7jCiSmZmJj169Ah3GGG3aEMhR46byLED2vHgKYPqdd1mNsU5N6ReVxomymGhNXfuXLZu1c2R62LgwIHExcWFO4ywcs5xypM/sGB9AeOvPIjs1IR6W3dN+UsjBIhISHVvmcZFB3Tnramr+GHxpnCHIyISNDPjzuP6UbC9lL99Mq/B3lfFmYiE3KUH96B9VjK3vDtLA6OLSFTp1SadC/bvyqs/reDnBup9ruJMREIuOSGW20b3ZcG6Qp75VgOji0h0+fOhe9AmI4mb35nVIAOjqzgTkQZxeJ/WHNq7FeM+X8iavG3hDkdEJGipiXEN2vtcxZmINJjbRvelrNxx5wfqXCIi0WVU/zaM3KMFD/xvAesLtof0vVSciUiD6Zidwh8P7sFHM9cycYGGKxKR6GFm3DGmH0Wl5dzzUWg7B6g4E5EGdeGB3ejSPIXb3ptNcak6B4hI9OjaIpULD+jG21NX8dPSnJC9j4ozEWlQiXGx3HpsXxZv3KLOASISdS45uDvtMpO49d3ZIRuaTsWZiDS4g3u34tDerXjoi4Wsyw9t2w0RkfqUkhDHjUf3Yc6afF7+cXlI3kPFmYiExc3H9KGkzHHvxw13Y0cRkfowqn8bhnVrzgP/m0/uluJ6X7+KMxEJiy4tUvnDAV1D3nZDRKS+mRm3j+lLwfZS/v6/+fW+fhVnIhI2lx7cg7YhbrshIhIKPVunc/awzrzy43Jmrcqr13WrOBORsPHabuwZ0rYbIiKh8pfDepKdksAt786ivB4PMFWciUhYHd2/LcO6Nef+T+eTE4K2GyIioZKZHM+1R/bm5+WbeXvqqnpbr4ozEQkrM+O20X0pLCrl/hC03RARCaUTB3dgYMcs7vl4HgXbS+plnSrORCTserUJXdsNEZFQiokx7hjdl01binjoi4X1ss64elmLiMhu+sthPfllfSFO/QJEJMoM7JjF+SO60i4ruV7Wp+JMRCJCZnI8L1ywb7jDEBHZJTcf06fe1qXLmiIiIiIRRMWZiIiISARRcSYiIiISQVSciYiIiEQQFWciIiIiEUTFmYiIiEgEUXEmIiIiEkF0nzMJSkxMDCkpKeEOI6okJiaGOwQR8SUlJYU7hKhjZuEOocky10hux21mG4BldVikBbAxROHsjkiNCyI3tkiNCyI3tkiNC+oWW2fnXMtQBtNQ6pjDGsv+a0iRGhdEbmyKq+7qJX81muKsrsxssnNuSLjjqCxS44LIjS1S44LIjS1S44LIji1SRPI2itTYIjUuiNzYFFfd1VdsanMmIiIiEkFUnImIiIhEkKZcnP0r3AFUI1LjgsiNLVLjgsiNLVLjgsiOLVJE8jaK1NgiNS6I3NgUV93VS2xNts2ZiIiISCRqymfORERERCKOijMRERGRCNIoizMzO9LM5pvZL2Z2XRXTE83sNX/6JDPrEjDtev/1+WZ2RCTEZWZdzGybmU3z/55o4LgOMLOfzazUzE6sNO0cM1vo/51Tn3HVQ2xlAdvsvQaO6wozm2NmM8zsCzPrHDAt3NusptjCuc0uNrOZ/nt/Y2Z9AqaF7HcZaZS/QhJbWHJYpOavIGMLSw6L1PwVZGz1l8Occ43qD4gFFgHdgARgOtCn0jyXAE/4j08FXvMf9/HnTwS6+uuJjYC4ugCzwri9ugADgOeBEwNezwYW+/828x83i4TY/GmFYdxmBwMp/uP/C9iXkbDNqowtArZZRsDj0cAn/uOQ/S4j7U/5K2SxNXgOi9T8VYfYGjyHRWr+qkNs9ZbDGuOZs32AX5xzi51zxcCrwJhK84wBnvMfvwEcambmv/6qc67IObcE+MVfX7jjCqVa43LOLXXOzQDKKy17BPCZcy7HOZcLfAYcGSGxhVIwcX3pnNvqP/0B6OA/joRtVl1soRRMXPkBT1OBit5KofxdRhrlrxDEFqYcFqn5K9jYwpHDIjV/BRtbveWwxlictQdWBDxf6b9W5TzOuVIgD2ge5LLhiAugq5lNNbOvzGxkPcUUbFyhWLYh1p9kZpPN7AczOy6McV0AfLyLyzZkbBDmbWZml5rZIuBvwGV1WbaRUP4KTWyhWDbU6w7VbxEiN4dFav4KOrb6ymEa+Dw6rAE6Oec2mdlg4B0z61upSpeddXbOrTKzbsB4M5vpnFvUkAGY2ZnAEODAhnzfYFQTW1i3mXPuUeBRMzsduAmo9zZ50uCUv3ZN2PMXRG4Oi8T8BfWXwxrjmbNVQMeA5x3816qcx8zigExgU5DLNnhc/qnQTQDOuSl416t7NmBcoVg25Ot3zq3y/10MTAD2asi4zOww4EZgtHOuqC7Lhim2sG+zAK8Cx+3istFM+Ss0sYVi2ZCuO4S/xaBjC0MOi9T8FXRsAXYvhwXTEC6a/vDOBi7Ga3RX0Wivb6V5LmXHhquv+4/7smOjvcXUX4Pa3YmrZUUceI0RVwHZDRVXwLzPsnNj2iV4jUKb+Y/rJa56iK0ZkOg/bgEspFLjzRDvy73w/hPao9LrYd9mNcQW7m22R8DjY4HJ/uOQ/S4j7W8380STy1/BxhYwb+U8EbLf427GFbLfYh32Z4PnsN2MKxK2Wb3lsHoJOtL+gFHAAn8H3ui/dgdelQ2QBPwXr1Hej0C3gGVv9JebDxwVCXEBJwCzgWnAz8CxDRzXULxr5FvwjtBnByx7vh/vL8B5YdiXVcYGDAdm+j+ImcAFDRzX58A6f59NA96LoG1WZWwRsM3+GfA9/5KAxBfK32Wk/e1qngj1dtrVuAhx/goytrDksF2NK9S/xSBjC0sO29W4ImSb1VsO0/BNIiIiIhGkMbY5ExEREYlaKs5EREREIoiKMxEREZEIouJMREREJIKoOBMRERGJICrOJOKYWXMzm+b/rTWzVf7jQjN7LNzxiYhUR/lL6oNupSERzcxuAwqdc/eHOxYRkbpQ/pJdpTNnEjXM7CAz+8B/fJuZPWdmX5vZMjMba2Z/M7OZZvaJmcX78w32B1ueYmafmlnb8H4KEWmKlL+kLlScSTTrDhwCjAZeBL50zvUHtgFH+wnuYbxhUQYDTwN3hStYEZEAyl9SrbhwByCyGz52zpWY2UwgFvjEf30m0AXoBfQDPjMz/HnWhCFOEZHKlL+kWirOJJoVATjnys2sxP3WgLIc77tteGPVDQtXgCIi1VD+kmrpsqY0ZvOBlmY2DMDM4s2sb5hjEhEJhvJXE6biTBot51wxcCJwn5lNB6YBw8MalIhIEJS/mjbdSkNEREQkgujMmYiIiEgEUXEmIiIiEkFUnImIiIhEEBVnIiIiIhFExZmIiIhIBFFxJiIiIhJBVJyJiIiIRJD/BxlVGwyTe7JEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.patches as patches\n",
    "\n",
    "fn_f = lambda t: 5 + -1*t**2\n",
    "\n",
    "h = 0.1\n",
    "eps = 0.2*h\n",
    "t0,t1,t2,t3 = 0,h,h*2,h*3\n",
    "tru_t = np.linspace( t0, t0+3*h, 40 )\n",
    "tru_f = list(map( fn_f, tru_t ));\n",
    "\n",
    "plt.figure( figsize=(10,4))\n",
    "ax = plt.subplot( 1, 2, 1 )\n",
    "plt.plot( tru_t, tru_f )\n",
    "plt.plot( t0, fn_f(t0), 'C0.', markersize=16)\n",
    "plt.plot( t1, fn_f(t1), 'C1.', markersize=16)\n",
    "plt.plot( t2, fn_f(t2), 'C2.', markersize=16)\n",
    "plt.plot( [t1,t2], [fn_f(t1),fn_f(t2)], 'C3')\n",
    "yb = plt.ylim()\n",
    "poly = patches.Polygon([(t1,yb[0]),(t1,fn_f(t1)),(t2,fn_f(t2)),(t2,yb[0])],facecolor='0.8')\n",
    "ax.add_patch(poly)\n",
    "plt.title( 'Trapezoid Slope from Current Step\\n(Trapezoidal Rule)')\n",
    "plt.ylabel( 'f(t)')\n",
    "plt.xlabel( 'Time')\n",
    "\n",
    "\n",
    "ax = plt.subplot( 1, 2, 2 )\n",
    "plt.plot( tru_t, tru_f )\n",
    "plt.plot( t0, fn_f(t0), 'C0.', markersize=16)\n",
    "plt.plot( t1, fn_f(t1), 'C1.', markersize=16)\n",
    "plt.plot( t2, fn_f(t2), 'C2.', markersize=16)\n",
    "slope = (fn_f(t1)-fn_f(t0))/h\n",
    "plt.plot( [t0,t2], [fn_f(t0),fn_f(t0)+slope*2*h], 'C3');\n",
    "yb = plt.ylim()\n",
    "poly = patches.Polygon([(t1,yb[0]),(t1,fn_f(t0)+slope*h),(t2,fn_f(t0)+slope*2*h),(t2,yb[0])],facecolor='0.8')\n",
    "ax.add_patch(poly)\n",
    "plt.title( 'Trapezoid Slope from Previous Step\\n(Used by Adams-Bashforth)');\n",
    "plt.xlabel( 'Time');\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine steps taken by AB2\n",
    "\n",
    "Use the slope at previous steps to get the final slope estimate for the current step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAEGCAYAAABmRJHZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAByW0lEQVR4nO3deVxU1fsH8M+ZGRj2fd9lGYZBQBRxX9D8pobmkqa4kOYSpvTLrGyztNW0zUrDct8wzTShNM0Fl8JQZIcBFFH2fWdgZu7vD5YQAVFBQJ7368XLmXvPPfPcK8rDueeeh3EcB0IIIYSQ3orX1QEQQgghhHQlSoYIIYQQ0qtRMkQIIYSQXo2SIUIIIYT0apQMEUIIIaRXE3R1AA/KyMiIs7Oz6+owCCGkR7l69Wo+x3HGj9iHiUAg+AlAX9Av06TnUAKIlcvliwYMGJDbUoMelwzZ2dkhIiKiq8MghJAehTF261H7EAgEP5mZmbkYGxsX8Xg8WpeF9AhKpZLl5eVJsrOzfwIwuaU2lNkTQghpr77GxsallAiRnoTH43HGxsYlqBvRbLnNY4yHEEJIz8ajRIj0RPXft63mPJQMEUIIIaRXo2SIEEJIj5Cdnc0Xi8USsVgsMTIy8jAxMXFveF9dXc064jMOHDig6+LiInF2dpY4ODi4btiwwait9itXrrRYs2aNaVttLl++rH7w4EHdhvf79u3Tffvtt806Il7SMXrcBGpCCCG9k5mZmSIxMTEeqEtCtLS0FOvWrctp2F9bWwsVFZWH7l8mk7FXXnnF9u+//05wcHCoraqqYlKpVPVR446IiNCIiIjQfP7550sAYM6cOSUASh61X9JxaGSIEEJIjzV9+nQ7Pz8/G3d3d3FAQIBV85EaJycn16SkJFUA2Lx5s4Gbm5uLWCyW+Pn52crl8rv6Ki4u5snlcmZqaioHAHV1dc7Dw0MGAElJSaqDBw8WiUQiyZAhQ0TJycn3JEne3t7OYWFhGgCQlZUlsLS0dKuurmaffvqpxfHjx/XFYrHkxx9/1N+0aZPh/Pnzbdrqd/r06XYvvPCCtaenp9jKysptx44d+p10CQl60chQSV4Vrv6RhlGzncFXoRyQEEIexeuHo6yl2WUaHdmnyEy7csNzHrcf9LisrCzVa9euJQoEAqxcudKipTbXrl1TO3z4sEFERESiUCjk5s6da/PDDz8YLl++vKChjampqWLcuHHFNjY27sOGDSudOHFiyZIlSwr5fD4CAgJs5syZU7BixYqCr7/+2jAgIMD69OnTqfeLTU1NjXvrrbcyIyIiNHfv3p0OAJs2bTJs2N9Wvzk5OSoRERGJ169fV5s6darjggULih702pD26RVZQWZKMQ5/FoGEy1k4uy8RHEcPQxBCyJNi2rRpRQJB27/bnzhxQjs2NlbDw8PDRSwWSy5evKhz48YNYfN2Bw8evHXixAmpl5dXxaZNm8xmzpxpBwCRkZGaS5YsKQSAgICAwqtXr2p1ROxt9Tt58uRiPp+PAQMGVBcUFDz8/T9yX71iZCjnRimqK2oBAEn/ZEPfTAMDxtt1bVCEENKDPcwITmfR0tJSNrwWCAScUtn4FjKZjAEAx3FsxowZBd9//33G/frz9vau8vb2rlqyZEmho6OjG4C09sQhEAg4hUIBAKisrHzkCd1qamqNv7nTL/Gdq1eMDPUbZw2XoeaN7/85egOp11pckZsQQkgPZmdnJ7t+/bomAFy8eFEjIyNDCADjx48vDQkJ0c/IyBAAQE5ODr/55OiSkhJeSEiIdsP78PBwdQsLixoA8PT0rPjpp5/0ASAoKMjAy8urvPlnW1tby65cuaIJAPv27Wuc46Ojo6MoLy9v8edte/olna9XJEOMMYzyc4aFk17jttM74pF7q7TrgiKEENLh5s+fX1RUVMR3dHR0/eabb0xsbW2rAWDAgAHV7777bsbYsWNFIpFIMmbMGNHt27fvuvWkVCqxYcMGUzs7u75isViybt06y23btt0EgB9++CF9z549RiKRSHLgwAHDzZs33zMytnr16pxt27YZu7i4SPLz8xvvvEyYMKFMKpWqN0ygbnpMe/olnY/1tKE3Ly8v7mFrk1WX1+Lw+giU5FUBADR0VTFjtRe09NU6MkRCCOl2GGNXOY7zepQ+oqKi0jw8PPI7KiZCHqeoqCgjDw8Pu5b29YqRoQZqWip45mV3CDXqEvbKkhqEbo5GTbX8PkcSQggh5EnVq5IhANA308TTS/qCx6ub25Z/uxynd8SDU/asETJCCCGEdIxelwwBgLXYACNnixrf34zKx99H77tcBCGEEEKeQL0yGQIA1xGW8HjKuvF95J/piL+U2YUREUIIIaQr9NpkCACGTnOEnVvjQqA4vy8JGUm0wCchhBDSm/TqZIjHYxj3oisMLesW/FQqOfyxNQbFOZVdHBkhhBBCHpdenQwBgKqaAM+87A51nbq1t2QVcoRujm5csZoQQkj3wefzB4jFYknD19tvv232MP1Mnz7drqH46fPPP2979erVVtdY2bdvn27D5+zZs0evpbZvvvmmWUNMTWP86KOPTB4mvo4UFham8cILL1i3tM/S0tItKyurzWoUrbX5+uuvDUUikUQkEkmcnJxc9+7dqwfU1V5LS0vrUeVDekU5jvvRNlDDxAA3HP0yEopaJYpzKnFiaywmBXqAz+/1+SIhhHQbQqFQmZiYGN+RfR48ePBWW/vnzJlTAqAEAI4ePaonl8tLBgwYUN20zfr167PXr1+fDQAaGhqezWNUKpXgOA58Pr8jQ2+XkSNHVo4cObJDb3mkpqaqfPHFF+bXr19PMDQ0VJSUlPAaEqa9e/ca9evXr8rOzq7HjCrQT/p6Zn10MdbfpfF9RlIRwoKlVA+GEEK6uYKCAr6dnV3fqKgoIQBMmjSpzxdffGEE1CUmL774orWjo6PrkCFDRJmZmfcMAnh7ezuHhYVpAMDhw4d1JBKJi7Ozs2TIkCEioG6kY/78+TanTp3SPH36tN67775rJRaLJXFxcfcUem0qKSlJ1c7Oru/UqVPtRCKRa2pqquqcOXNs+vbt6+Lo6Oj66quvWjS0tbS0dHvppZesRCKRxM3NzSU2NlYIAJmZmYKnn37aoW/fvi59+/Z1+fPPPzUBYNSoUY4No0/a2tr9vv32W8PKykr23HPP2YlEIomLi4vk+PHj2gAQEhKi7ePj4wgA2dnZ/GHDhjk5Ojq6Pv/887ZNf8Y99dRTDq6uri6Ojo6uGzduNGrr3LKyslQ0NTWVurq6CgDQ1dVVisXimh07dujHxsZqzJ8/314sFkvKy8vZhQsXNAYOHOjs6urqMnz4cKdbt26pNFz3BQsWWIvFYomTk5Pr2bNnNe7zV91paGSoCScvUxTnVOLK8ZsAgPgLmdA31UC/p2y6ODJCCOmGtvo437PNZVIhRqzMg6ych12TnO7Z7/F8Pga9VICybAEOzHa4a9+Ss0n3+0iZTMYTi8WShvevvfZa1uLFi4u++uqrdH9//z7Lli3LKS4uFrz22mv5AFBVVcXz8vKq2LZt2+1Vq1aZr1692mL37t3pLfWdmZkpWL58ud25c+cSxWJxTU5Ozl3DOOPGjat46qmnin19fUsWLFjQrqdt0tPThdu2bbs5duzYNAD48ssvM0xNTRVyuRxDhw51Dg8PVx80aFAVAOjq6sqlUmn8d999Z7hixQrrs2fPpixdutR65cqVOU8//XR5cnKy6tNPP+1048aNuPPnz6cAwIULFzRefPFFOz8/v+L169ebMMYglUrjIyMj1SZOnOiUmpoa2zSe1atXWwwZMqR848aNWcHBwbo///xzY9Kzb9++NFNTU0V5eTnz9PSUzJ07t8jMzEzR0nkNHjy40sjIqNba2tpt2LBhZdOmTSvy8/MrWbBgQdGWLVtMNm7ceHvkyJGVMpmMBQYG2oSGhqZYWFjIf/zxR/1Vq1ZZHjp0KK3h7ycxMTH+jz/+0FqyZEmf5OTkuPZc145GyVAzXhPtUJxTCemVHADApV9SoGuigT7ubSbJhBBCHoPWbpNNnTq19Oeff9Z/4403bK9evdr4A5XH42HRokWFALBw4cKCadOmObbW97lz5zS9vb3LxGJxDQCYmpq2mAg8CHNz85qxY8dWNLzftWuXwc6dO43kcjnLy8tTiYqKUmtIhvz9/QsBYPHixYXvvvuuNQBcunRJJzk5Wb3h+PLycn5JSQlPV1dXmZWVJXjhhRf6BAcHpxoaGiouX76stWLFilwA8PT0rLawsKiJiYm5a37TP//8o33kyJEUAJg1a1bJ0qVLG89x/fr1pqGhoXoAkJ2drRIXF6dmZmZWgRYIBAKEhYUlnz9/XuPPP//UWb16tXVERITml19+edcaNdHR0cLk5GT1MWPGiIC624XGxsaNt8/8/PwKAWDChAnl5eXlvPz8fL6RkdEjX/cHRclQM4wx+MwTozS/Gtk3SgAO+HNbHKa/3h9GVtr374AQQnqLtkZyhFrKNvdrm8nbMxLUXgqFAlKpVE1NTU1ZUFAgcHBwaHG+CmOsoz6yXTQ0NJQNrxMTE1W/++4706tXryYYGxsrpk+fblddXd04XYXH+2/mCmOMAwCO43Dt2rUEDQ2Nu+ZsyOVyTJ8+3f7NN9/MHDhw4F3zlx5GSEiI9vnz57UjIiIStbW1ld7e3s5VVVVtTqXh8Xjw8fGp9PHxqZwwYULpokWL7JonQxzHMUdHx6rr168nttRH87+Px/3304DmDLVAoMLHhJfcoG1Yl1DLZQqEfh+NihJZF0dGCCGkJevWrTMViUTVO3fuvLFw4UI7mUzGgLqRiIanxnbu3Gno7e1d1lofo0ePrrhy5Yp2YmKiKgA0v00GAFpaWorS0tKH+tlZVFTEV1dXVxoYGChu374tOHfunG7T/bt37zYAgG3btul7enpWAMDw4cNLP/3008Yn0i5fvqwOAC+//LKVRCKpXLJkSePtumHDhpXv3bvXAKgbkcnKylJ1d3e/K1EaPHhw2c6dOw0B4Oeff9YpLS3lA0BxcTFfV1dXoa2trYyMjFSLiorSbOtc0tLSVC5evNg4xyciIkLD0tKypuEalZSU8AHA3d29urCwUHD69GlNAJDJZCwiIqJxtOrAgQP6AHDy5EktbW1thaGh4WMfFQJoZKhVGjqqeOZld/zy+VXUVitQXiTD71tiMHWlJwSqj/9pAEIIIffOGRozZkzJ0qVL8/fs2WN09erVBH19feXhw4fLVq9ebf7VV19lqqurK69cuaK5YcMGC0NDw9ojR47caK1vCwsL+aZNm9KmTp3qqFQqYWhoWHv58uXkpm3mzJlTGBAQYPfDDz+YHj58ONXV1bXdvyUPGTKkqm/fvpUODg59zc3NawYMGFDedH9RURFfJBJJVFVVueDg4BsAsHXr1tuLFi2yEYlEEoVCwQYNGlQ2dOjQ9K1bt5o6OjpWi8ViHQB47733Mt54443c+fPn24pEIgmfz0dQUFCaurr6XSNKn332Web06dPtHR0dXb28vMrNzc1rAGD69OklW7duNba3t3e1t7ev9vDwaPH2WIOamhq2atUqq5ycHBWhUMgZGBjU/vjjj+kAMH/+/PwVK1bYvv7668qIiIiE4ODg1MDAQJuysjK+QqFgAQEBOV5eXtUAoKamxrm4uEjkcjnbunXrzfZey47GetrTUl5eXlxERMRj+7xbcQUI/S4KDZfJcYAJ/veiKxiva4byCCHkYTDGrnIc5/UofURFRaV5eHjkd1RMj4OGhoZnZWVlZFfHcT+WlpZuERERCebm5vKujuVx8fb2dm6YaP04Pi8qKsrIw8PDrqV9NDJ0H7auhhg+U4QLB6WoAYedsRlYteYOyuUK6GuoYO4QWywd6QBNIV1KQgghpCfqtJ/gjLHtAHwB5HIc17eF/QzANwAmAqgE8ALHcdc6K55H4e5jheyMMrwVdRPFPA6K+ry9sLIWQedv4ERsNn5dNowSIkII6WZ6wqgQAGRkZMR0dQyP25UrVzpsAv2j6swJ1DsBjG9j/wQATvVfSwBs6cRYHlmMMUMJn4Oi2d0xmVyJWwWVCApL7ZrACCGEEPJIOi0Z4jguDEBhG02eBbCbq/MPAD3GmHlnxfOo9oWno7UbuTK5Env/aXENL0IIIYR0c135aL0lgNtN3t+p33YPxtgSxlgEYywiLy/vsQTXXFFl2yVWiiprHlMkhBBCCOlIPWKdIY7jtnIc58VxnJexsXGXxKCv0XYBXn0N1ccUCSGEEEI6UlcmQxkArJu8t6rf1i3NHWILoaDly8XngBHK0sccESGE9D58Pn9AQ4FSsVgsSUpKUvX09BQ/SB/r1q0zKSsra/XnX1ZWlkAgEPT//PPP7/rt29LS0k0kEknEYrFEJBJJ9u7dqwcAKSkpKoMGDRI5ODi4Ojo6un744YcmLfUbFRUl9Pb2dhaLxRJ7e3vX2bNn2wJ1CykePHhQt6VjyOPRlcnQbwDmszqDAZRwHJfVhfG0aelIB9gaatyTEPE5QE/JYJelh+tnbrdyNCGEkI7QUJus4cvZ2bkmMjLynlIPtbWtT20ICgoyLS8vb/Xn3+7du/U9PDwqDh06ZNB83/nz56WJiYnxhw4dSn3jjTesAUBFRQVffPHFndTU1Lh///03Ydu2bSZXr15Va37syy+/bBMYGJiTmJgYf+PGjbhXX301F6hbvTk0NJSSoS7UackQY+wAgL8BODPG7jDGXmSMvcQYe6m+ye8AbgBIAfAjgGWdFUtH0BQK8OuyYVg6yh4GmqpgDDDQVMV4fR3MLRNCFQwXf07GrfPhQE2bC3cSQgjpQBoaGp5AXX2tAQMGOI8ZM8bRycmpb2lpKW/06NGOzs7OEicnJ9cff/xR/6OPPjLJzc1VGTVqlGjQoEGilvo7dOiQwcaNG2/n5OSopKamtjhHori4mK+jo6MAAFtb29rhw4dXAoC+vr7SwcGhKj09/Z65E7m5uSq2traNE0y9vb2rqqur2aeffmpx/PhxfbFYLPnxxx/1S0tLeTNmzLBzc3NzcXFxaRyB2rRpk+HYsWMdvL29nW1tbfu+9tpr3faho56m0xbG4Thu9n32cwBe7qzP7wyaQgFWjnPGynHOjdtqaxQ4+sU15N4qAwNwMjgP06+/AMMXNwNaXTO/iRBCHofZIbOdm28bazu2cJHboryK2greopOLnJrv93XwzZ/jMqcgrzJPEHgm0KHpvgO+B+677kzTchzW1tayU6dO3bWuSXx8vEZkZGScWCyu2blzp56ZmVntuXPnUgCgoKCAb2hoqNiyZYvp+fPnpS2t9pySkqKSl5en4uPjUzl58uSi3bt3G6xduzanYf+oUaNEHMexO3fuqG7fvv2e0h5JSUmq8fHxGqNGjSpvvu/ll1/OmThxosjT07Ni7NixJS+//HKBkZGR4q233sqMiIjQ3L17dzoALF++3NLHx6f00KFDafn5+XwvLy+XyZMnlwJAdHS0ZkxMTJyWlpbS09NT8uyzz5Y8rhWcn2Q9YgJ1d6aiysfEZe7Q1BcCAGo5DYRKp6EqaApQQGsPEUJIR2p6m6x5IgQA7u7uFWKxuAYA+vfvX3XhwgWdgIAAyxMnTmi1pwjo7t27DSZPnlwEAPPmzSv85Zdf7rpVdv78eWlycnJcRERE/KpVq2xKSkoaf46WlJTwpk2b5vDZZ5/dNjAwUDbv+5VXXimIiYmJmzZtWmFYWJj2wIEDxVVVVffUdjp37pzOV199ZS4WiyXDhw93lslkLCUlRRWoK9xqZmam0NLS4p555pmic+fOabXnupG20ZLJHUBTV4hnlrnjyMZrkMsUKFOYIvTWPEz9aSL4c/YDVgO6OkRCCOlwbY3kaKpoKtvab6xhLG/PSNCD0tDQaExC3N3dZdeuXYv/5ZdfdN977z3L06dPl27cuLHNuam//PKLQV5ensqRI0cMgLpbWzExMUI3N7e7CrK6urrKDA0Na69du6bm4+NTKZPJ2DPPPOMwY8aMQn9//+LW+rezs6v9v//7v4L/+7//K3BycnKNiIhQb96G4zgcPnw4xcPD467PvHjxomZd8Yb/NH9PHg6NDHUQY2ttjFsgAeq/L3NqRDiV9wK4a7u7NjBCCOml0tLSVLS1tZXLli0rXLlyZfb169c1AEBTU1PRdESnQXR0tLCiooKfm5sbnZGREZORkRGzfPny7F27dt0zkTojI0Nw584doaOjY41SqcSsWbNsRSJR9QcffJDTvG2Dw4cP68hkMgYA6enpguLiYr6trW2Njo6OoumEbh8fn9IvvvjCVKmsy+suXbrUmDBdvHhRJycnh19eXs5+//13vZZux5EHR8lQB7LvZ4yhUx0b36eWD8KvJS/WvZHR9yshhDxOV69eVe/Xr5+LWCyWfPzxxxZr1qzJAgB/f//88ePH3zOBeteuXQYTJ04sarpt1qxZRQ2jREDdnCGxWCwZNWqU85o1a+5YW1vLT506pXX06FHDixcvajc88t/So/InTpzQcXZ2dnV2dpaMGzdOtHbt2js2NjbyCRMmlEmlUvWGCdSfffZZplwuZ2KxWOLo6Oj67rvvNi5I7O7uXjF58mQHV1dX10mTJhXRfKGOwermMfccXl5eXERERFeH0SqO43B2byISLv03EiuZZAKf1PmA2wxg9GqAhjUJIY8ZY+wqx3Fej9JHVFRUmoeHR35HxUQezKZNmwybTrQmDyYqKsrIw8PDrqV9NDLUwRhjGDXbGZYivcZt0SE5uKk2Djj/GfDbCjSWvSeEEEJIl6NkqBPwBTyMX+oGXZO627wCjuGXqP8hx301ELkHCJ5NaxERQgh5IIGBgQU0KtQ5KBnqJGqaKvB92QNCjboH9tQVDD+dH4QSny+AlNPAybe7OEJCCCGEAJQMdSo9Uw2MX9IXPF7dHCE9GfBVqCMqnwsGxrzXxdERQgghBKBkqNNZiQ0wyu+/RVqNS5X45LAGatUMAHkN8Msi4M7VLoyQEEII6d0oGXoMJMMt0O8p68b35lm1+HDTFSjLc4HbV4BdvoD0ZBdGSAghhPRelAw9JkOmOcLO3ajxvVFSBT47mgfuxVOAkQg4MBu4uqsLIySEkO7vzTffNHN0dHQViUQSsVgsOXPmjObjjiEkJETbx8fHsbX9CxcutDYxMXFXKP6r/rFp0yZDfX19j4a1g8aPH29fVlbGA4APPvjA1MHBwVUkEkmGDBkikkql9xR5BVo/93Xr1pk09EUeDl28x4THYxi3UAJDq7oyMnwwqIYXYvPZQuCFUMDBBzgeCPz9fRdHSggh3dPp06c1T548qRcTExMvlUrjz549K7W3t6+5/5GPj0KhwIkTJ/TMzc1rfv/9d+2m+yZNmlSUmJgYn5KSEqeiosJt375dHwAGDBhQef369QSpVBo/ZcqUoldffdWqeb9tnXtQUJBp0xWsyYOji/cYqaoJ8Mwyd2jo1CX96hxD4YkM7LuSB8wOBgYuBuxGdHGUhBDSPWVkZKgYGBjI1dXVOQAwNzeX29nZ1QLAhQsXNAYOHOjs6urqMnz4cKdbt26pAEBsbKxw6NChImdnZ4lEInGJi4sTKpVKLF261MrJyclVJBJJfvzxR32gbsTH29vbefz48fZ9+vRxnTx5cp+GkhiHDx/W6dOnj6tEInE5fPiwXmsxhoaGajs5OVUtWrQob//+/feU8QCA2tpaVFZW8gwMDBQAMGnSpDJtbW0lAAwfPrw8KyvrnpGh1s79o48+MsnNzVUZNWpU44raR44c0enXr59YIpG4TJgwwb6h9IilpaXbSy+9ZCUSiSRubm4usbGxwof7m3jyUDL0mGkbqGFigDv4KnWX3kDJQ/TBFPwRmwc8sxEwd69reG0PrUVECOnWbs6Y6dz8K3/rVmMAUFZU8FraX7hnjyEAyPPyBM333e/zpkyZUpqZmalqZ2fXd+7cuTahoaFaACCTyVhgYKDNsWPHUuPi4hL8/f3zV61aZQkAfn5+fV566aXcpKSk+IiIiEQbG5va3bt368XExKgnJCTE/fXXX9I1a9ZYNSRPCQkJ6t9///3tlJSUuPT0dOGpU6e0Kisr2fLly+1+++23lNjY2ITc3FyV1mLcv3+/wcyZMwvnzJlT9Ndff+k21CIDgOPHj+uLxWKJmZmZR3FxsWD27NnFzY8PCgoyfuqpp0rae+7vvvturomJSe358+el4eHh0qysLMEnn3xiHhYWJo2Pj0/o379/5Ycffmja0I+urq5cKpXGL126NHfFihXWzT+nt6JkqAuY9tHBWH+Xxvc2cj5+2x6HSyl5dRty4upume30BcrzuihKQgjpXnR1dZWxsbHx33333S1jY2O5v7+/w6ZNmwyjo6OFycnJ6mPGjBGJxWLJhg0bzDMzM1WKiop4OTk5qvPnzy8GAA0NDU5bW1t54cIF7ZkzZxYKBAJYW1vLBw0aVH7x4kUNAHBzc6twcHCo5fP5cHV1rUxNTVW9fv26mpWVlczNzU3G4/EwZ86cgpbiq66uZmfOnNH18/MrNjAwUPbr16/iyJEjOg37G26T5eXlRbm4uFStWbPGrOnxmzdvNoiKitJYu3ZtdnvPvXm7c+fOaaampqp5e3uLxWKxJDg42DA9Pb1xpMnf378QABYvXlwYGRmp9ZB/FU8cQVcH0Fs5eZmiJLcS4b/dBAD0lfHx7vZIfPvyIPS1dAWe3wccXghsGwfM/QUwdOjiiAkh5G59Dv2c1No+nqamsq39AmNjeVv7Wz1OIICvr2+Zr69vmbu7e9WePXsMBw8eXOno6Fh1/fr1xKZti4qKHvgXfqFQ2Fiwk8/nQy6Xt7uY5JEjR3TKysr4ffv2dQWAqqoqnpqamnL27Nl3jfTweDxMnjy5+Pvvvzdp2Hb06FHtjRs3ml+4cCGp4VZYcy2de2Bg4F2JGcdxGD58eOnx48dvttQHj/ffJWGM9azipJ2IRoa60IAJdhB5m0KgwsPgec6o0RHAf/sV3MyvAMQTAf/jQHUJsO1/tBYRIaTXi4qKEsbExDTOc4mMjFS3srKqcXd3ry4sLBScPn1aE6i7bRYREaGmr6+vNDMzq9mzZ48eAFRVVbGysjLeyJEjyw4fPmwgl8uRmZkpuHLlitaIESNanZfQr1+/6oyMDNW4uDghAAQHB7c4F+jAgQMGX3/99a2MjIyYjIyMmLS0tJiLFy/qtPSk14ULF7Tt7OxkAHDp0iX1FStW2B47dizF0tKyxeKVrZ07AGhqaioa5gWNHj26IiIiQqthPlBpaSkvOjq68bjdu3cbAMC2bdv0PT09aS5GPRoZ6kKMMfjME6P/07YwtNTCHpEunvvhb8zbFo5fAobC1Hog8OIpYN9zQMltwGpAV4dMCCFdprS0lB8YGGhTWlrK5/P5nJ2dnWzXrl231NTUuODg4NTAwECbsrIyvkKhYAEBATleXl7Ve/fuvbl48WLbDz/80EJFRYU7dOhQ6rx584ovX76s5eLi4soY49auXXvHxsZGHh0d3eLnamhocN9+++0tX19fR3V1deWgQYPKy8vL+U3blJWV8cLCwnR37dp1q2Gbjo6O0svLqzw4OFgXaJwzpKVUKmFubl6zf//+NAB4/fXXrSsrK/kzZsxwAAALC4uaM2fOpLTn3AHA398/f/z48SJTU9Oa8PBwaVBQUNqsWbPsa2pqGAC8//77Ge7u7jIAKCoq4otEIomqqioXHBx8o6P+bno6xnE9a5TMy8uLi4iI6OowOk30nWLM3voPrPQ18PPSIdDVUAHkMkBQn9gX3gAM7Ls2SEJIj8MYu8pxnNej9BEVFZXm4eGR31ExkcfL0tLSLSIiIsHc3LzF0acnXVRUlJGHh4ddS/voNlk3426lh63zvXAzvwIv7voXVTWK/xKhO1eB7wYCZz8FelgSSwghhHRXlAx1Q8McjfD1rH64ml6El/dfQ62ibp0LmLsD7s8D5z8DflsBKHplck8IIeQhZGRkxPTWUaH7oWSom5roZo6PpvTFmcRcvHk4GkolB/BVgGe/B0a+AUTuAYJn01pEhBBCyCOiZKgbmzPIFq+NE+FIZAY+/j0BHMcBjAFj3gF8vwJSTgNRB7o6TEIIIaRHo6fJurnlYxxRUFGDbRdvwlBLFctG19cG9FoIWPQHzD3q3iuVAI9yW0IIIeRB0U/Pbo4xhjW+Ekz2sMDnJ5IQfCX9v50W/epGigpSgS1DaC0iQggh5CFQMtQD8HgMG2d4YKTIGG//GoMTsfes1A7Iq4FdvkDSiccfICGEPCZvvvmmmaOjo6tIJJKIxWLJmTNnNAHA29vbOSwsTKOzPjcyMlKtX79+YlVV1f5r1qwxba3d9u3b9e3t7V0biqY+jI8++sjE3t7edfLkyX327dun+/bbb5vd/6g6SUlJqj/88EOLi0I2WLdunYlQKOxfUFDQuFZSSEiItra2dj+xWCwRiUSSoUOHijIyMgQAsGXLFgORSCQRiUQST09P8d9//63eUr9ff/21YUM7Jycn17179+oBwKZNmwzT0tJarefWHVAy1EOoCnj4YW5/uFvpITA4EpdTmyz1YehQtzijkahuUvXVnV0WJyGEdJbTp09rnjx5Ui8mJiZeKpXGnz17Vmpvb1/zOD7bxMRE/s0336QvXbo0p612O3bsMNqyZcut8PBwaXv6ra2tvWfbtm3bjE+dOiX97bffbs6ZM6fkk08+uec34JaOA4Dk5GThwYMH20yGDh8+bNC3b9+KhmSlgZeXV3liYmK8VCqN9/T0rNi4caMJADg6OsouXbqUJJVK4996663MpUuX2jbvMzU1VeWLL74w//vvv5OkUml8REREgpeXVyUA7N271yg9PZ2SIdIxNFQF2PHCQNgaaGDhzn9xLin3v51aJsALoYDDWOD4K0DM4a4LlBBCOkFGRoaKgYGBvKF2l7m5udzOzu6erCAoKMigYXQiICDAsmG7hoaG54svvmjt6OjoOmTIEFFmZqYAAOLi4oQjRoxwcnV1dRkwYIBzZGSkWvM+LS0t5aNGjapUUVFpdZG3VatWmV+9elVr6dKldkuXLrWqrKxkzz33nJ1IJJK4uLhIjh8/rg3UjZSMGTPGcfDgwaKhQ4c6N+3Dz8/P5s6dO8IJEyY4rV271mTTpk2G8+fPtwGA6dOn2/n5+dm4u7uLAwICrEJDQ7XEYrFELBZLXFxcJEVFRbx33nnHMiIiQkssFkvWrl1r0jzGuLg4YWVlJX/dunUZP//8c4tJk1KpRFlZGV9fX18OAOPGjaswNjZWAICPj09Fdna2avNjsrKyVDQ1NZW6uroKoK6wrFgsrtmxY4d+bGysxvz58+3FYrGkvLycXbhwQWPgwIHOrq6uLsOHD3e6deuWClA3urdgwQJrsVgscXJycj179mynjfQ1RxOoexh9TVXsXzwY/tuvYNGuCHz5fD9M9rCo2ynUAmYfAC5vApwndm2ghJAn2vcvnem0+kAv/zCmxQmQU6ZMKf30008t7Ozs+g4fPrx09uzZhc8880x50zZpaWkqH3zwgeXVq1cTjI2N5SNGjBDt2bNHb968ecVVVVU8Ly+vim3btt1etWqV+erVqy12796dvmjRItutW7fecnNzk505c0YzICDA5p9//mnXyE5TGzduzAoLC9PZuHHj7ZEjR1a+//77powxSKXS+MjISLWJEyc6paamxgJAXFycRnR0dJypqamiaR/79+9PP3/+vO758+el5ubm8uaV6bOyslSvXbuWKBAIMGbMGMdNmzbd+t///ldRUlLC09DQUH788ccZX3zxhenZs2fvKufRYPfu3fpTp04tHD9+fPnixYvVbt++LbC2tpYDQEMSVVxcLFBXV1d8/fXXd5of/+233xr5+PiUNN8+ePDgSiMjo1pra2u3YcOGlU2bNq3Iz8+vZMGCBUVbtmwxabgmMpmMBQYG2oSGhqZYWFjIf/zxR/1Vq1ZZHjp0KA2oK26bmJgY/8cff2gtWbKkT3JyctyD/j08DBoZ6oGMtYUIXjoY/W318UpwJPb8nfbfTr4KMOI1QFUDqC4FTrwFyMpb7YsQQnoKXV1dZWxsbPx33313y9jYWO7v7+/QPFm4ePGi5uDBg8ssLCzkKioqeP755wvPnz+vBdRVbF+0aFEhACxcuLDgypUrWiUlJbzIyEitGTNmOIjFYsmyZctsc3NzO+SWzuXLl7XmzZtXAACenp7VFhYWNTExMWoAMGLEiNLmiVB7TJs2rUggqBvHGDx4cPmqVausP/roI5P8/Hy+isr9wz5y5Ijh/PnzC/l8PiZOnFi0Z88e/YZ9DbfJsrOzo/38/AqWL19u1fTY48ePa+/du9fom2++uSdJEggECAsLS96/f3+qk5NT9erVq61Xrlxp0bxddHS0MDk5WX3MmDEisVgs2bBhg3lmZmZj4H5+foUAMGHChPLy8nJefn4+v3kfnYFGhnooHTUV7F7ojeX7r+G9Y3EorKhF4FhHMMb+a3TrMhD+A5D+N+B3CNAy7rqACSGkAwgEAvj6+pb5+vqWubu7V+3Zs8cwMDCw4GH6YoxBoVBAW1tbnpiYGN/RsbZFQ0ND+TDHaWlpNR73ySefZE+ZMqXk2LFjuiNGjBCHhoYmt3XslStX1G/duiUcP368CABqa2uZlZVVzdtvv53XvO306dOLGwrHAkB4eLj6smXLbENDQ5PNzMxaTOJ4PB58fHwqfXx8KidMmFC6aNEiuy+//DKzaRuO45ijo2PV9evXE1vq466fYS287yyUDPVgaip8bJk7AG8ejsZXp6UoqqzBGl8JeLz6bx7n8cCs/cChBcC2ccDcX+omWxNCyCNq7VZWZ4qKihLyeDy4ubnJACAyMlLdysrqrgnUI0aMqHjjjTess7KyBMbGxvJDhw4ZLFu2LBeomwuzY8cO/SVLlhTt3LnT0Nvbu8zAwEBpZWVVs337dv2FCxcWKZVKhIeHqw8ZMqTqUeMdNmxY+d69ew0mT55cFh0dLczKylJ1d3evDg8P75C5MHFxcUJvb+8qb2/vqqtXr2rExsaq2dnZ1ZSXl7c4mrJ7926D1157LfPTTz9tnJBtaWnpJpVK75kDdPbsWS1bW1sZACQnJ6vOmDHDYfv27Tfd3d1lLfWdlpamcufOHZXhw4dXAkBERISGpaVlDQBoaWkpSkpK+ADg7u5eXVhYKDh9+rTmU089VSGTyVhMTIzQy8urGgAOHDigP2nSpLKTJ09qaWtrKwwNDR949OxhUDLUw6nwedg4wwN6GqrYfukmSqpq8flz7lDh198BdZ4AvBAC7J/5X0Jk4dm1QRNCyEMoLS3lBwYG2pSWlvL5fD5nZ2cn27Vr162mbWxtbWvff//9jFGjRok4jmNPPfVU8dy5c4sBQF1dXXnlyhXNDRs2WBgaGtYeOXLkBgAcOHDgxuLFi23Xr19vLpfL2dSpUwubJ0Pp6emCgQMHSioqKviMMS4oKMg0ISEh1sDAoNURnjfeeCN3/vz5tiKRSMLn8xEUFJTWMPm7I3z++ecmly9f1mGMcc7OzlXPPfdcCY/HA5/P55ydnSV+fn7577//fuOTNkePHjU4fvz4XaNHEyZMKNq1a5fBkCFDKhrmDHEcB21tbcX27dvTAODdd981Ly4uFqxYscIWAAQCARcbG5vQtJ+amhq2atUqq5ycHBWhUMgZGBjU/vjjj+kAMH/+/PwVK1bYvv7668qIiIiE4ODg1MDAQJuysjK+QqFgAQEBOQ3JkJqaGufi4iKRy+Vs69atNzvqWt0P43pY9XMvLy8uIiKiq8PodjiOw3dnUvDFKSnGik3w/Zz+UFNp8stBQSpwbDkwbSugZ911gRJCugRj7CrHcV6P0kdUVFSah4dH/v1bdk8aGhqelZWVkV0dB2mZt7e3c8NE687oPyoqysjDw8OupX00gfoJwRjDirFO+HBKX5xJysX8bVdQUtXkiVNDB2DB73WJkFIJpJ7pumAJIYSQbqRTkyHG2HjGWBJjLIUxtrqF/TaMsbOMsUjGWDRjjJ4Hf0TzBtti0yxPRN4uwqyt/yCvrMnt3YaJaNd2AXumAmc/AXrYyCAhhDwsGhXq3q5cuZLUWaNC99NpyRBjjA/gewATAEgAzGaMSZo1exfAzxzHeQKYBWBzZ8XTm0zysMBP/gORll+BGT9cxu3CZt9bnnOBfnOB8+uB35YDipZXMiWEEEJ6g84cGfIGkMJx3A2O42oABAN4tlkbDoBO/WtdAJkgHWKUyBh7Fw1CUWUtnvvhMpKyy/7byVcBnv0OGPkGELkXODCb1iIihBDSa3VmMmQJ4HaT93fqtzX1AYC5jLE7AH4HsKKljhhjSxhjEYyxiLy8e5ZDIK0YYKuPn5cOAccBM4P+xrX0ov92MgaMeQfw/bpuPaL8pC6LkxBCCOlKXT2BejaAnRzHWQGYCGAPY+yemDiO28pxnBfHcV7GxrRw4INwNtPGLwFDoaehgjk/hiNM2iyZ9FoA/F80YFm/sn5V8WOPkRBCCOlKnZkMZQBo+gy3Vf22pl4E8DMAcBz3NwA1AEadGFOvZG2ggUMvDYGdkSZe3PUvQqKb3Y3UrL/kUQeBbwcAdx77WmqEENIub775ppmjo6OrSCSSiMViyZkzZzSBuseyw8LCOq2w55YtWwxEIpFEJBJJPD09xX///bd6S+22b9+ub29v7zpo0CDRw37WRx99ZGJvb+86efLkPvv27dN9++23zdp7bFJSkuoPP/zQZtX6devWmQiFwv4FBQWN66+EhIRoa2tr9xOLxRKRSCQZOnSoKCMjQwC0/9xbsmnTJsO0tLRuXbEe6Nxk6F8AToyxPowxVdRNkP6tWZt0AGMBgDHmgrpkiO6DdQITbTUELxmMftZ6WHEgEvvCb93byMqrrtjrLl8g6cTjD5IQQtpw+vRpzZMnT+rFxMTES6XS+LNnz0rt7e1r7n/ko3N0dJRdunQpSSqVxr/11luZS5cutW2p3Y4dO4y2bNlyKzw8vF2FXmtr732AZdu2bcanTp2S/vbbbzfnzJlT8sknn2S35zgASE5OFh48eLDNZOjw4cMGffv2rdi7d69e0+0NtcmkUmm8p6dnxcaNG02A9p97S/bu3WuUnp7ee5MhjuPkAJYDOAkgAXVPjcUxxtYxxibXN3sNwGLGWBSAAwBe4HraKpA9iK66CnYvHAQfZxO882ssvj+bgrsut6ED8OIpwEgEBM8Gru7sslgJIaS5jIwMFQMDA3nDKs7m5uZyOzu7e7KCoKAgA5FIJHFycnINCAhonKuqoaHh+eKLL1o7Ojq6DhkyRJSZmSkA6spajBgxwsnV1dVlwIABzpGRkWrN+xw3blyFsbGxAgB8fHwqsrOz7ylhsWrVKvOrV69qLV261G7p0qVWlZWV7LnnnrMTiUQSFxcXyfHjx7WButGSMWPGOA4ePFg0dOhQ56Z9+Pn52dy5c0c4YcIEp7Vr15ps2rTJcP78+TYAMH36dDs/Pz8bd3d3cUBAgFVoaKiWWCyWiMViiYuLi6SoqIj3zjvvWDasJL127VqT5jHGxcUJKysr+evWrcv4+eefW0yalEolysrK+Pr6+vL2nrtcLsf06dPtnJycXEUikWTt2rUmO3bs0I+NjdWYP3++vVgslpSXl7MLFy5oDBw40NnV1dVl+PDhTrdu3VIB6kb2FixYYC0WiyVOTk6uZ8+e7bRRvpZ06pwhjuN+5zhOxHGcA8dxH9dvW8Nx3G/1r+M5jhvGcZwHx3H9OI77szPjIYC6Kh9B8wZgSj8LbDiZhI9DE6BUNkmItEyAF0IBh7HA8VeAzOtdFishhDQ1ZcqU0szMTFU7O7u+c+fOtQkNDdVq3iYtLU3lgw8+sDx37pw0Pj4+LjIyUnPPnj16AFBVVcXz8vKqSElJiRs2bFjZ6tWrLQBg0aJFtps3b06Pi4tL2LBhw52AgACbtuL49ttvjXx8fEqab9+4cWNW3759K3fv3n0jKCjozvr1600YY5BKpfH79++/sWTJErvKykoGAHFxcRrHjh1L/ffff+96emX//v3pJiYmtefPn5c2LaXRICsrS/XatWuJP/30050vvvjCbNOmTbcSExPj//nnn0QtLS3lxx9/nNEwwtPS8bt379afOnVq4fjx48tv3rypdvv27cayXA1JlIWFhfuFCxe0ly9ffs9q462d+99//62RlZWlkpycHCeVSuNffvnlggULFhQ1XI/ExMR4FRUVBAYG2hw7diw1Li4uwd/fP3/VqlWNyWpVVRUvMTExftOmTbeWLFnSp62/g47W1ROoSRdQ4fPw5cx+eGGoHX66eBOvH46GXNGkvI5QC5h9oK7Iq0W/LouTEEKa0tXVVcbGxsZ/9913t4yNjeX+/v4OmzZtMmza5uLFi5qDBw8us7CwkKuoqOD5558vPH/+vBZQV1V90aJFhQCwcOHCgitXrmiVlJTwIiMjtWbMmOEgFosly5Yts83NzW31ts7x48e19+7da/TNN9/cuV+8ly9f1po3b14BAHh6elZbWFjUxMTEqAHAiBEjSk1NTR+4COm0adOKBIK6/GXw4MHlq1atsv7oo49M8vPz+Soq978bdeTIEcP58+cX8vl8TJw4sWjPnj36Dfsakqjs7OxoPz+/guXLl1s1PbatcxeLxbLbt28L/f39rQ8fPqyjr69/z7lFR0cLk5OT1ceMGSMSi8WSDRs2mGdmZjYG7efnVwgAEyZMKC8vL+fl5+e3WHC2M1Ch1l6Kx2N4f5IE+hqq+Oq0FKXVtfh2tud/9cz4KoD4mbrXt/8FLm8CpmypS5QIIaSLCAQC+Pr6lvn6+pa5u7tX7dmzxzAwMLDgYfpijEGhUEBbW1uemJgYf7/24eHh6suWLbMNDQ1NNjMze6Rq6hoaGq0WeG2LlpZW43GffPJJ9pQpU0qOHTumO2LECHFoaGhyW8deuXJF/datW8Lx48eLAKC2tpZZWVnVvP322/fM1Z0+fXrxjBkzHBre3+/cjY2NFbGxsfG//vqrzg8//GB88OBBg0OHDqU1bcNxHHN0dKy6fv16YkvxsYYqCa2870w0MtSLMcbwylNOWDvZFafic+C//QpKq1uYlFeQDCSG1E2sLqf57YSQrhEVFSWMiYkRNryPjIxUt7KyumsC9YgRIyrCw8O1s7KyBHK5HIcOHTIYPXp0OVA3F2bHjh36ALBz505Db2/vMgMDA6WVlVXN9u3b9RvatPS0VHJysuqMGTMctm/fftPd3V3WfH9Lhg0bVr53714DoG5UJCsrS9Xd3b364a/A3eLi4oTe3t5VH3/8cba7u3tFbGysmq6urqK8vLzFEZXdu3cbvPbaa5kZGRkxGRkZMbm5udE5OTkqUqn0njlAZ8+e1bK1tZUB7Tv3rKwsgUKhwAsvvFD86aefZsTExGgAgJaWlqKkpIQPAO7u7tWFhYWC06dPawKATCZjERERjfOzDhw4oA8AJ0+e1NLW1lYYGho+UsL5IGhkiMB/qB30NFTw2s9RmL31H+xa6A0jLeF/Dfr5Aer6wKEFwLZxwNxf6iZbE0LIY1RaWsoPDAy0KS0t5fP5fM7Ozk62a9euux6NtbW1rX3//fczRo0aJeI4jj311FPFc+fOLQYAdXV15ZUrVzQ3bNhgYWhoWHvkyJEbAHDgwIEbixcvtl2/fr25XC5nU6dOLRwyZEhV037fffdd8+LiYsGKFStsAUAgEHCxsbEJbcX7xhtv5M6fP99WJBJJ+Hw+goKC0homf3eEzz//3OTy5cs6jDHO2dm56rnnnivh8Xjg8/mcs7OzxM/PL7/pvKGjR48aHD9+/K7RowkTJhTt2rXLYMiQIRUNc4Y4joO2trZi+/btae0997S0NJUXX3zRTqlUMgBYt27dHQCYP39+/ooVK2xff/11ZUREREJwcHBqYGCgTVlZGV+hULCAgIAcLy+vagBQU1PjXFxcJHK5nG3duvVmR12n9mA97eEtLy8vLiIioqvDeCKdTcxFwL6rMNdVx54XvWGl32wy/50IYP/MuteLTgMG9o8/SELIQ2GMXeU4zutR+oiKikrz8PC4Z1JtT6GhoeFJxVq7J29vb+eNGzfe7sxCrVFRUUYeHh52Le2j22SkkY/YBHtfHISCchme2/I3knPK7m5g5VX36H3f5wC9di8zQQghhHRrdJuM3MXLzgAHlw7B/O1XMCPob+xc4I1+1nr/NTB0ACZ+Xve6NBO4GQZ4zOqSWAkhXeeL530HdFbfrx0M6ZRl8GlUqPu6cuVKlxbIpJEhcg8Xcx0cfmkIdNRU4PfjP7iY3Mqo+OVvgV+XAmc+BnrY7VZCCCGkASVDpEW2hpo4/NIQ2BhoYOHOf/FHTNa9jcatAzznAmGfA8eWA4qWl4cnhJCO0lW1ySIjI9X69esnVlVV7b9mzRrT1tp1RG2yrhYSEqLt4+Pj2Nr+hQsXWpuYmLgrFP897LVp0yZDfX19D7FYLHF0dHQdP368fVlZGQ8APvjgA1MHBwdXkUgkGTJkiKilp9das27dOpOGfjoT3SYjrTLRUcPBJUOwcNe/eHn/NXw81Q2zvZsszMpXASZ/B+hYAufXA+U5wIydtBYRIb1AZ93KakvT2mTq6upcVlaWQCaTPZbFaExMTOTffPNN+uHDh/XbatdQm+zpp58ub0+/tbW1aM9iid2FQqHAiRMn9MzNzWt+//137UmTJjVOLp00aVLR7t270+tf99m+fbv+K6+8UjBgwIDK1157LUFbW1u5fv1641dffdUqNDT0Rns+LygoyHTx4sWF2traD7UuU3vRyBBpk66GCva86I0RTsZ460gMtpxLvbsBY4DP28Ckb4CybEBJo0OEkM7RlbXJLC0t5aNGjapUUVFpdU5AR9QmKy0t5Y0ePdrR2dlZ4uTk5Prjjz/qA0BrNb1iY2OFQ4cOFTk7O0skEolLXFycUKlUYunSpVYNdcIa+ggJCdH29vZ2Hj9+vH2fPn1cJ0+e3EeprMsxDh8+rNOnTx9XiUTicvjwYb3WzjE0NFTbycmpatGiRXn79+9vsbZZbW0tKisreQYGBgoAmDRpUllDMjN8+PDyrKyse0aGWjrvjz76yCQ3N1dl1KhRooaRtiNHjuj069dPLJFIXCZMmGBfUlLCq//7cXvppZesRCKRxM3NzSU2NlbY/DPaQskQuS8NVQF+nO+FyR4WWH8iEZ/+noB7lmQY8AKw5GzdekS11UDhY10ighDSC3SX2mSt6YjaZEeOHNExMzOrTUpKik9OTo6bNm1aqUwmY63V9PLz8+vz0ksv5SYlJcVHREQk2tjY1O7evVsvJiZGPSEhIe6vv/6SrlmzxqoheUpISFD//vvvb6ekpMSlp6cLT506pVVZWcmWL19u99tvv6XExsYmtFWOZP/+/QYzZ84snDNnTtFff/2l23Rk7vjx4/pisVhiZmbmUVxcLJg9e3Zx8+ODgoKMn3rqqXtqm7V03u+++25uQ5228PBwaVZWluCTTz4xDwsLk8bHxyf079+/8sMPP2y8ZamrqyuXSqXxS5cuzV2xYoX1g/zdUTJE2kVVwMPXz/fDvMG2CAq7geUHIlEhk9/diF//7+ePN4CfxtatS0QIIR2kO9QmexAPU5usf//+VRcuXNAJCAiwPHHihJahoaGitZpeRUVFvJycHNX58+cXA4CGhganra2tvHDhgvbMmTMLBQIBrK2t5YMGDSq/ePGiBgC4ublVODg41PL5fLi6ulampqaqXr9+Xc3Kykrm5uYm4/F4mDNnTovlTaqrq9mZM2d0/fz8ig0MDJT9+vWrOHLkiE7D/kmTJhUlJibG5+XlRbm4uFStWbPGrOnxmzdvNoiKitJYu3ZtdnvOu3mbc+fOaaampqp5e3uLxWKxJDg42DA9Pb1xlMnf378QABYvXlwYGRn5QPM1aM4QaTcej2Hds66w0FPHhpOJkGaX4Yd5A+Bg3Ox7btgrwM3zwE7fujlEzuO7JF5CyJOnK2uTdaTWapO5u7vLrl27Fv/LL7/ovvfee5anT58unTlzZnFLNb2KiooeeEBDKBQ2Duvz+XzI5fJ2z7k6cuSITllZGb9v376uQN1Im5qamnL27Nl3jfTweDxMnjy5+Pvvvzdp2Hb06FHtjRs3ml+4cCGppVW4WzrvjRs33vXkDsdxGD58eOnx48dbvPXA4/13ORhjD/SIM40MkQfCGEPAaAfseXEQCipq8Ox3l3AitlmSb+hQtzijiRgIng1E7OiaYAkhT5SurE32MB6mNllaWpqKtra2ctmyZYUrV67Mvn79ukZrNb309fWVZmZmNU1uA7KysjLeyJEjyw4fPmwgl8uRmZkpuHLlitaIESMqWvvMfv36VWdkZKjGxcUJASA4OLjFuUAHDhww+Prrr2811DZLS0uLuXjxok5LT3tduHBB287OTgYAly5dUl+xYoXtsWPHUiwtLeX39tzyeQOApqamomFe0OjRoysiIiK0GuYDlZaW8qKjoxu/H3bv3m0AANu2bdP39PRs9XxbQiND5KEMczTC8RXDsWzvVby09yoCRjvgtXEiCPj1/ya0TAD/EODQC8BfawHJs4BGi/++CCGkXbqyNll6erpg4MCBkoqKCj5jjAsKCjJNSEiINTAwaPUpp4epTXb16lX1t956y4rH40EgEHCbN2++paamxrVW02vv3r03Fy9ebPvhhx9aqKiocIcOHUqdN29e8eXLl7VcXFxcGWPc2rVr79jY2Mijo6Nb/EwNDQ3u22+/veXr6+uorq6uHDRoUHnzYq9lZWW8sLAw3abXW0dHR+nl5VUeHBysCzTOGdJSKpUwNzev2b9/fxoAvP7669aVlZX8GTNmOACAhYVFzZkzZ1Lud94A4O/vnz9+/HiRqalpTXh4uDQoKCht1qxZ9jU1NQwA3n///YyG4rFFRUV8kUgkUVVV5YKDg9v1tFoDqk1GHolMrsDa4/HYH56OYY6G2DTLE4ZNi7wqausmUxuL6hZmVCoAPuXghDxuVJuMapM9ySwtLd0iIiISzM3NWxx5Aqg2GelEQgEfn0x1w+fPuePftCL4fnsR128X/9eAr4LsUoZzu3+E/Ox64MAsQNau5TcIIYSQx4KSIdIhZnpZ40jAUPB5DDN/+Bv7w9PBcRxklRUI+WY9roYew4GjsShKuATsfAYoz+3qkAkhvQyNCj25MjIyYtoaFbofSoZIh+lrqYvjy4djsIMh3v41Bm8cjkbkqZMoyambYJ2bU4y9twcjMTkX2DYOKEi9T4+EEEJI56NkiHQofU1V7HhhIALHOuHQ1Tv4MM0I/Wa8AL6gbp5QjawWobcd8KdUC7U/TaRbZoQQQrocJUOkw/F5DCvHibDN3wu3i6qwMk4HTovehp6peWObmHx97L89EAV5hV0YKSGEEELJEOlEY11McXzFcJjrqmH5n3konxAI5yEjG/fn5xVj79uvInb/RnD/bu/CSAkhhPRmrSZDjLHfGWN2jzEW8gSyNdTEr8uGYUo/S3x1Ph2/6o7G8AXLIFCpW0FdLpPh5LFz+GPrD6g5sa7u8XtCCGnFm2++aebo6OgqEokkYrFYcubMGU0A8Pb2dg4LC9PorM/dsmWLgUgkkohEIomnp6e4tYUZt2/frm9vb+/aUFi0JwoJCdH28fFxbG3/woULrU1MTNwViv8qZmzatMlQX1/fQywWSxwdHV3Hjx9v37AY4wcffGDq4ODgKhKJJEOGDBFJpdJ7CrW2Zt26dSYtLerY0dr6gB0A/mSMvcMY65A6LaR3Ulfl48uZHvjwWVeEpeTjjRh1DH71QxhY/ldHL6HUBHv3nUPuziV1axMRQkgzp0+f1jx58qReTExMvFQqjT979qzU3t6+5v5HPjpHR0fZpUuXkqRSafxbb72VuXTpUtuW2u3YscNoy5Ytt8LDw6Xt6be2tmf9f6dQKHDixAk9c3Pzmt9//1276b6G2mQpKSlxKioqXMOq3gMGDKi8fv16glQqjZ8yZUrRq6++atXezwsKCjItLy/v9GSo1dXvOI47xBj7A8B7ACIYY3sAKJvs/7KzgyNPDsYY5g2xg8RCF8v2XYX/0XR8/PwqWESGIPbsKQBAUY0G9p/MwKg7M9Bv9UEwFeF9eiWEdJU7qy8M6Ky+rT4bcbWl7RkZGSoGBgbyhlWcW3uUOigoyOCLL74wa1iBesuWLRlA3aKLs2fPzj9//ryOsbFx7S+//HLDwsJCHhcXJ3zppZdsCgsLBWpqasqffvrplqen511lM8aNG9dY3sHHx6di+fLl94xurFq1yvzq1ataS5cutXv66aeLv/rqq4z58+fbRkdHa/D5fHz++ee3J02aVLZp0ybDo0eP6ldWVvIUCgVrWrm+tLSUN3nyZPusrCxVpVLJ3njjjczFixcXXbhwQWPlypXWlZWVPH19ffm+ffvSbG1ta2NjY4VLliyxLSgoEPD5fO7QoUM3XFxcZAEBAVZnzpzRZYxxr7/+etbixYuLQkJCtNetW2dhYGBQm5SUpO7m5lZ59OjRmzweD4cPH9Z5/fXXrdXV1ZXe3t6tPtkSGhqq7eTkVPXcc88V7d+/32DSpEllzdvU1taisrKSZ2BgoACApm2GDx9efvDgQcPmx7R03jk5OSq5ubkqo0aNEunr68vDw8OlR44c0Vm3bp1FTU0Ns7W1lQUHB6fp6uoqLS0t3SZNmlR05swZHaFQyB04cOBG3759Za2dR3P3y7ZqAFQAEALQbvZFyAMbYKuPkBUj4GGlh9d+TcQl87H437KVUFGrG3FWcDyciZXj+KaNqK6gJ80IIf+ZMmVKaWZmpqqdnV3fuXPn2oSGht5TmTwtLU3lgw8+sDx37pw0Pj4+LjIyUrNJ7S6el5dXRUpKStywYcPKVq9ebQEAixYtst28eXN6XFxcwoYNG+4EBATYtBXHt99+a+Tj41PSfPvGjRuz+vbtW7l79+4bQUFBd9avX2/CGINUKo3fv3//jSVLlthVVlYyAIiLi9M4duxYatNECKgrhmpmZlablJQUn5ycHDdt2rRSmUzGAgMDbY4dO5YaFxeX4O/vn79q1SpLAPDz8+vz0ksv5SYlJcVHREQk2tjY1O7evVsvJiZGPSEhIe6vv/6SrlmzxurWrVsqAJCQkKD+/fff305JSYlLT08Xnjp1SquyspItX77c7rfffkuJjY1NyM3NbfVu0P79+w1mzpxZOGfOnKK//vpLVyaTNRZ6rS/HITEzM/MoLi4WzJ49u7j58UFBQcZPPfXUPdeupfN+9913c01MTGrPnz8vDQ8Pl2ZlZQk++eQT87CwMGl8fHxC//79Kz/88EPThj50dXXlUqk0funSpbkrVqywbv4ZbWlrztB4ANcBaADoz3Hc+xzHrW34epAPIaQpY20h9i0ahCUj7bH771t4P06I8e+sh7GdfWOb5Ct/Y8+qAGRF/NWFkRJCuhNdXV1lbGxs/HfffXfL2NhY7u/v77Bp06a7RhkuXryoOXjw4DILCwu5iooKnn/++cLz589rAXVVzRctWlQIAAsXLiy4cuWKVklJCS8yMlJrxowZDmKxWLJs2TLbtpKB48ePa+/du9fom2++uXO/eC9fvqw1b968AgDw9PSstrCwqImJiVEDgBEjRpSampoqmh/Tv3//qgsXLugEBARYnjhxQsvQ0FARHR0tTE5OVh8zZoxILBZLNmzYYJ6ZmalSVFTEy8nJUZ0/f34xUFdjTFtbW3nhwgXtmTNnFgoEAlhbW8sHDRpUfvHiRQ0AcHNzq3BwcKjl8/lwdXWtTE1NVb1+/bqalZWVzM3NTcbj8TBnzpyCls6nurqanTlzRtfPz6/YwMBA2a9fv4ojR47oNOxvuE2Wl5cX5eLiUrVmzRqzpsdv3rzZICoqSmPt2rXZzftu6bybtzl37pxmamqqmre3t1gsFkuCg4MN09PTG0fo/P39CwFg8eLFhZGRkfckym1pq0jUOwBmcBwX9yAdEtIeAj4Pb090QT9rPbx+KAp+P1fimwVvwfLyUVw/GQoAKC0sQvCGLzHcNxZec1aA8ejhR0K6i9ZuZXU2gUAAX1/fMl9f3zJ3d/eqPXv2GAYGBrb4w/t+GGNQKBTQ1taWJyYmxt+vfXh4uPqyZctsQ0NDk83MzO75Yf0gNDQ0Wizw6u7uLrt27Vr8L7/8ovvee+9Znj59unTmzJnFjo6OVdevX09s2raoqOiB/1MUCoWNT6nw+XzI5XLWVvumjhw5olNWVsbv27evK1A30qampqacPXv2XSM9PB4PkydPLv7+++9NGrYdPXpUe+PGjeYXLlxIaqlYbUvnvXHjxqymbTiOw/Dhw0uPHz9+s6X4eE1+RjDGHuhpnFYvJMdxIygRIp1tops5ji0fBl11FczfHYkbjk9j0sq3INTQBAAowRAWcgq/rXkJNSfWAJ/bAx/o1f159hNatJGQXiQqKkoYExPTOJkwMjJS3crK6q4J1CNGjKgIDw/XzsrKEsjlchw6dMhg9OjR5QCgVCqxY8cOfQDYuXOnobe3d5mBgYHSysqqpmGyr1KpREtPiiUnJ6vOmDHDYfv27TcbqqTfz7Bhw8r37t1rAADR0dHCrKwsVXd39+q2jklLS1PR1tZWLlu2rHDlypXZ169f13B3d68uLCwUnD59WhMAZDIZi4iIUNPX11eamZnVNLkNyMrKyngjR44sO3z4sIFcLkdmZqbgypUrWiNGjKho7TP79etXnZGRoRoXFycEgODgYIOW2h04cMDg66+/vpWRkRGTkZERk5aWFnPx4kWdlp72unDhgradnZ0MAC5duqS+YsUK22PHjqVYWlq2OM+rpfMGAE1NTUVJSQkPAEaPHl0RERGhFRsbKwTq5hlFR0c3fj/s3r3bAAC2bdum7+np2er5toTKh5Mu52iijWPLh+P1Q1H4+PcEXHczx7sffomzW75EVkoSVJgCw+THwbtcDfDqf5mqLAAufQPE/wYsOg0IH2hElBDSA5WWlvIDAwNtSktL+Xw+n7Ozs5Pt2rXrVtM2tra2te+//37GqFGjRA0TqOfOnVsMAOrq6sorV65obtiwwcLQ0LD2yJEjNwDgwIEDNxYvXmy7fv16c7lczqZOnVo4ZMiQqqb9vvvuu+bFxcWCFStW2AKAQCDgYmNjE9qK94033sidP3++rUgkkvD5fAQFBaW1NCrS1NWrV9XfeustKx6PB4FAwG3evPmWmpoaFxwcnBoYGGhTVlbGVygULCAgIMfLy6t67969NxcvXmz74YcfWqioqHCHDh1KnTdvXvHly5e1XFxcXBlj3Nq1a+/Y2NjIo6OjW/xMDQ0N7ttvv73l6+vrqK6urhw0aFB5eXk5v2mbsrIyXlhYmG7T662jo6P08vIqDw4O1gUa5wxpKZVKmJub1+zfvz8NAF5//XXryspK/owZMxwAwMLCoubMmTMp9ztvAPD3988fP368yNTUtCY8PFwaFBSUNmvWLPuamhoGAO+//35GQ3JaVFTEF4lEElVVVS44OPhGW9e5Ocb1sHVdvLy8uIiIiK4Og3QCjuPw44Ub+OyPRNgba2HzbA9k/XUEKn9/hYGGd6DCa2FUWaAGDHsF8Hn78QdMSA/CGLvKcZzXo/QRFRWV5uHhkd9RMT1uGhoanlSs9clkaWnpFhERkdBWsdaoqCgjDw8Pu5b20SQM0m0wxrBkpAP2LhqEoooaTN3yDyo9JmCQZUnLiRAAyKuBf7c93kAJIYQ8USgZIt3OUAcjhAQOh8hMGwH7roFXU9r2AVUPNXeSENLL0KjQkysjIyOmrVGh+6FkiHRL5rrqCF4yGPMG26KQa3s+UBWEqJZVtdmGENIhlEqlst1PHxHSXdR/37Zyi4GSIdKNCQV8fDilL3JEc1HNtbzsR62Sh8hcE2x5azFK83Mfc4SE9DqxeXl5upQQkZ5EqVSyvLw8XQCxrbXp1KfJ6hdu/AYAH8BPHMd91kKbmQA+AMABiOI4zq8zYyI9j+S5d1H9w1+QFd2EEP/V8ZGDj5IaISIKrKDkirHnjUA4+U2G0sEA/7P7H3iMcn1COpJcLl+UnZ39U3Z2dl/QL9Ok51ACiJXL5Ytaa9BpT5MxxvgApADGAbgD4F8AszmOi2/SxgnAzwDGcBxXxBgz4TiuzV/v6WmyXkpWjurzX6Hmnx+hpShFOV8HzOtFJFU7IOzQz+CU/41+xtuVomSwEV4esAKjrUeDMfollpCOeJqMkCdVZ2b23gBSOI67wXFcDYBgAM82a7MYwPccxxUBwP0SIdKLCbWg9r/3oPXuLez6XyS8a4Mw9Mpg3LAZheff/wzaRsaNTSVpOnA9KcM7ISsx9/e5+Df73y4MnBBCSHfXmcmQJYDbTd7fqd/WlAiAiDF2iTH2T/1ttXswxpYwxiIYYxF5eXmdFC7pCXg8hgXD+uDk/42ExFwHbxyOxluXyzHu7fVw8Brc2E6niOG5f+zAS8qHtEgKoG4dI0IIIaS5rr7nKwDgBGA0gNkAfmSM6TVvxHHcVo7jvDiO8zI2Nm6+m/RCtoaaOLB4MNY964qItEJM2hqJypHz4PPCEvAFdVPhOFkt+oWrwPB8Hmpl1TiachQBpwMQX3DfEkSEEEJ6kc5MhjIAWDd5b1W/rak7AH7jOK6W47ibqJtj5NSJMZEnCI/HMH+IHU7+30i4W+ni7V9j8VWmKca8/hH0TM0b28WdPY19b69EdU4hovOi8XzI81h5biVSi1O7MHpCCCHdRWcmQ/8CcGKM9WGMqQKYBeC3Zm2Oom5UCIwxI9TdNnugeiKEWBtoYN+iQfh4al9cu1WE2UfuQDD9NTgPGdnYpuBOOnKDfsfXRm/iJfeXcDnzMqYem4rN1zd3YeSEEEK6g05LhjiOkwNYDuAkgAQAP3McF8cYW8cYm1zf7CSAAsZYPICzAF7nOI6WEyYPjDGGOYNscfLVkehvq4/3/0jFPvUR8JqzFAIVVQCAvEaGsJ+2os/FKvw24Qhe6PsCJIYSAEB5TTmyK7K78hQIIYR0ESrUSp44HMfh4L+38XFoAhQch9e8dMA7tweFGf/N59c3t4Dv/62GiZ09AGDz9c34KeYnzHSeiUVui2CkbtRV4RPSKejRekJaR8kQeWJlFlfhrSMxOC/Nw2AbLcyoiUDa32cb9/MFAoyavwj9/vcMsiqyEBQdhGMpx6DKV4Wf2A8L+i6ArlC3C8+AkI5DyRAhraNkiDzROI7D4at3sC4kHrUKJV6xK0Nt2CHUVv9Xy8zJeyj+91Ig1DS1kFaShs1Rm3Hi5gk8ZfsUvhz9ZRdGT0jHoWSIkNZRMkR6heySarz9awzOJOZiuLESozNPouROWt1OxjDj3Y9h09e9sb20SAoBE8Bezx4Z5Rk4lXYKs8SzoCZQ65oTIOQRUTJESOu6ep0hQh4LM101bPP3wpczPRBTLsSXwqeh4lb3tNmgKTPvSoQAQKQvgr1e3Xyi07dO44urX2DikYk4kHgAtYrae/onhBDSc1EyRHoNxhim9bfCqVdHYpjYHF+WuyLGfTZMRk1u8zh/V3/seHoHrLWt8Un4J/D91RfHU48/pqgJIYR0NkqGSK9joqOGrfMG4JtZ/XBdYQzf7y9jy7lUyBXKVo/xMvPCzvE78cNTP0BfTR+RuZGN+3rarWZCCCF3ozlDpFfLK5PhvaOxOBGXDQ8rXWyY4QGRqXabx3Ach2pFNdQF6rieex0f/vMhlvdbjtHWo8EYe0yRE/JgaM4QIa2jkSHSqxlrC7Flbn985+eJ20VV8N10Ed+fTWlzlIgxBnWBOgCgWlGNank1As8GYs7vc3A58zKNFBFCSA9DI0OE1Csol2HNb3EIjc5CX0sdbHjOAy7mOvc9rlZZi+Opx7ElaguyK7IxxnoMvhnzzWOImJD2o5EhQlpHyRAhzfwRk4X3jsWipKoWy32csMzHASr8+w+i1ihqcEh6CHzGxyzxLCg5JZKLkuFs4PwYoiakbZQMEdI6SoYIaUFhRQ3WHo/DseuZcDHXwcYZ7nC1eLDVqE/cPIHXw17HONtxWOaxDI76jp0ULSH3R8kQIa2jOUOEtMBAUxXfzPJE0LwByCuT4dnvLuHLU1LUyFufS9TcMMthCPAIwOXMy5j22zS8deEt3C69ff8DCSGEPFY0MkTIfRRX1mDd8XgcicyA2Ewbnz/nDncrvfYfX12M7XHbcSDhAMy1zHHs2WP01Bl57GhkiJDWUTJESDv9lZCDt3+NQW6ZDM/1t8Kqp51hqtP+8hx5lXnIrsiGm7EbKmsrsTV6K+ZK5sJI3agToyakDiVDhLSOkiFCHkBJVS2+P5uCnZfSwOcxLB1ljyUj7aGhKnigfi7cuYAVZ1ZAla+K2eLZWNh3IXSFDzYniZAHQckQIa2jZIiQh5BeUIn1JxIRGpMFUx0hXvufM6b3twKf1/7bX7dKb2Hz9c344+Yf0FTRxHzX+VjsthgC3oMlVoS0ByVDhLSOkiFCHsHVW4X4MCQB128XQ2Kug3efccFQxwe77SUtkmLz9c0okZVg+9PbwRiDklOCx+j5BtJxKBkipHWUDBHyiDiOw/HoLKz/IxEZxVV4ysUEqye4wNFE64H6kSlkEPKFyK3Mhf8f/pjvOh/TnaZDla/aSZGT3oSSIUJaR796EvKIGGOY7GGBv14bhTfHi/HPjUI8/XUY3j8Wi8KKmnb3I+QLAQDlteUw0TDBJ+GfYNKvk/Br8q+QK+WdFT4hhPR6NDJESAfLL5fh69NSHLhyGxqqfKwY4wj/oXYQCvjt7oPjOPyd+Tc2RW5CXEEc7HXtcdD3INQE7X96jZCmaGSIkNZRMkRIJ0nOKcMnvyfgbFIerA3UsXq8Cya6mT3QGkMcx+HM7TNIKkzCsn7LAACx+bFwNXSltYrIA6FkiJDWUTJESCe7kJyHj0MTkJhdhgG2+njnGRf0t9F/qL6SCpPw3PHn4GbkhuWeyzHEfAglRaRdKBkipHWUDBHyGCiUHA5fvY2Nf0qRVybDJA8LvPG0M6wNNB6oH7lSjt9Sf8MPUT8gqyILA0wHINAzEP1N+3dS5ORJQckQIa2jZIiQx6hCJkfQ+VRsvXADSg5YOKwPlvk4QEdN5YH6qVHU4LD0MLZGb0WlvBKnZ5yGjqpOJ0VNngSUDBHSOkqGCOkCWSVV2HAyCUeuZcBQUxX/N06E2QOtIeA/2AOeVfIqxBfEY4DpAHAch02Rm/BMn2fgqO/YSZGTnoqSIUJaR8kQIV0o5k4JPgqNR/jNQjiaaOGdiS4Y7Wz8UPOA0kvTMTNkJiprKzHRfiKWeSyDjY5NJ0RNeiJKhghpHSVDhHQxjuNwKj4Hn/6RiJv5FRjuaIR3nnGBi/mD3/Yqri7G9rjtOJBwALXKWvib+sFfNA8Gfcw7IXLSk1AyREjrKBkipJuokSuxL/wWvvkrGSVVtZg5wBqv/U8EE50HX1sorzIPP8X8BMk5IwwoEkMo0ofWKEuoOzzcU2yk56NkiJDWUTJESDdTUlmLb88kY9ffaVDh8/DSKAcsHmEPddX2L9oIALU5Fcj56tpd2/KMymAzwR0GEkt6JL+XoWSIkNZRMkRIN3WroAKf/ZGIP2KzYaajhtefdsZUT0vweO1LYuQFVSj58xaqovOAZv/MiwwqYfW0K/TcLMHa2R/p2SgZIqR1lAwR0s39m1aIj0LiEXWnBH0tdfD602KMdDJq98hObX4Vys7dRuW1XEB59793zlgFhmPtoe5uTEnRE46SIUJaR8kQIT2AUsnheHQmPj+RhIziKrhZ6mLZaAc87WrW/pGiomqUhd1Bxb/ZgPzuf/cKPR70fGyhM8ACTED1m59ElAwR0jpKhgjpQWRyBY5GZuCH8zdwM78C9saaCBjlgCmellBp5xpFirIalF3IQMU/meBqlHf3r6mA0RgH6HhbgKk82Bwl0r1RMkRI6ygZIqQHUig5/BGbhc1nUxGfVQoLXTUsGWmP5wfatHuitaKiFuWXM1F+KQNcteKufbXqShiO7gOtwRbgCQWdcQrkMaNkiJDWUTJESA/GcRzOSfOw5WwqrqQVwlBTFQuG2WHeEDvoqrevxIeyWo7yf7JQfiEDyorau/YxdQG0h1mA562HXal7cTDpIIplxdAT6uF55+exoO8CaKg8WH010jUoGSKkdZQMEfKE+DetEJvPpuBsUh60hALMHWyLF4f3gbG2sF3HK2sUqPg3G2Vhd6AsqWncXsWq8X99NiBbWIAa/LddyBfCStsK+yfup4SoB6BkiJDWdepMScbYeMZYEmMshTG2uo120xljHGOM/qES8pAG2hlgxwJvhAYOx2hnY2wNS8Xw9Wfw3tFY3C6svO/xPFU+tIdZwvz1gdCf5gS+Qd1ij4cNTyNbJf+uRAgAZAoZ7pTdwY7YHZ1yPoQQ8rh02sgQY4wPQApgHIA7AP4FMJvjuPhm7bQBhAJQBbCc47g2h31oZIiQ9rmZX4Gg86n45dodKDngWQ8LBIx2gJOpdruO5xQcqqLz8L/rz6KEV95qO32hPsJmhXVU2KST0MgQIa3rzJEhbwApHMfd4DiuBkAwgGdbaPchgPUAqjsxFkJ6nT5GmvhsujvC3vDBC0Pt8EdsNsZ9FYYluyNw/XbxfY9nfAYNTxOU8irabFdcXYya7NaTJUII6e46MxmyBHC7yfs79dsaMcb6A7DmOC60rY4YY0sYYxGMsYi8vLyOj5SQJ5i5rjre85Xg0uoxCBzrhPCbhZjy/SXM+ekfXErJx/1Gh/WEem3u11FoIvfrSOTviUfJib+hKG87eSKEkO6my1ZXY4zxAHwJ4LX7teU4bivHcV4cx3kZGxt3fnCEPIEMNFWxcpwIl1aPwdsTxUjOKcecn8IxZfNlnIzLhlLZclL0vPPzEPJbnoStqlTBM0UjAQDVcQUoOyfH7Vf24c5rH6PszBlwNTUtHkcIId1JZyZDGQCsm7y3qt/WQBtAXwDnGGNpAAYD+I0mURPSubSEAiwZ6YCwN3zw8dS+KKqowdI9V/H012E4cu0OahV3L8S4oO8CWGlb3ZMQCflCWGtbYa7x83dtFxi7ACojUbD/JlInL0LpX2c7/ZwIIeRRdOYEagHqJlCPRV0S9C8AP47j4lppfw7AKppATcjjJVcoERpTt4BjUk4ZrPTVsXSkPWZ4WUOtfhXqytpK7Ijd0eo6Q7XZFSg9e7vForACYxXojneCvDABlZcvQ8f3GahJJO2urUY6Bk2gJqR1nbrOEGNsIoCvAfABbOc47mPG2DoAERzH/das7TlQMkRIl1EqOZxJzMXmcym4ll4MIy0hXhzeB3MH20BbrX0LOLZVFJapylAZvg/yW+FQ7WMHHd9noPvMM1C1s+uEsyHNUTJESOto0UVCyF04jkP4zUJ8fzYFF5Lzoa0mgP8QOywYZgdDrfYt4NhWUVgmlEORcxkVYcFQsTSHw58nwRiDsqoKPHX1zjglAkqGCGkLJUOEkFbF3CnB5nMpOBGXDaGAh1kDbeA/1A59jDTbdXxbRWF52gKoOTLoTxsIcAokj/aBmtgZOs/4QnvcU+Dr6HTGKfValAwR0jpKhggh95WSW44fzqfiaGQG5EoOQx0MMdvbBk+7mkFVcP/nMP4rCpsJrlp+1z79552h5qSJwu3bURISgtr0dDAVFWiNHgXDpS9Bva9rZ51Wr0LJECGto2SIENJuuaXVOHT1Dg5cScedoioYaqriOS8rzB5oA7t2jBY1FoW9mAFleS34+kKYrfIC49clVBzHoTomBiUhISj9/Q9Yff0VNLy8ILtxE7VZmdAcNAhMIOjs03wiUTJESOsoGSKEPDClkkNYch72h6fjr8RcKJQchjsaYba3DcZJTO87WtRQFJavqQKNfiYttuHkcoDHA+PxkPPpZyjctQt8IyPoTJgAXd9noObuTk+kPQBKhghpHSVDhJBHklNajZ//vY3gf28jo7gKRlqqmOFljdkDbWBj2DHV7JUyGcrPn0dpSCjKz50DV1MDNVdX2B0+RAlRO1EyREjrKBkihHQIhZJDmDQP+8LTcSYxB0oOGOFkBD9vGzwlMYUKv2PWeFWUlaHs1Gkoy0ph4O8PjuOQEfgK1Pt5QGfiRKiYm3fI5zxpKBkipHWUDBFCOlxWSRV+/vcODv6bjsySahhrCzHTywqzBtrA2qBjRosaKIqLkb5kKaqjowEAGl5e0PH1hc74p8HX0+vQz+rJKBkipHWUDBFCOo1CyeFcUi4OXEnHmcRccABGOhljtrcNxrqYdNhoEQDU3LqFktBQlB4PQc3Nm7D8+ivojB8PRXk5GI8HnkbHJmE9DSVDhLSOkiFCyGORWVyFg//exsF/byO7tBom2kI8P9Aazw+0hpV+xyUqHMdBlpAAVXt78NTUkP9DEPK3boX22LHQ9X0GmkOHgqm0b0XtJwklQ4S0jpIhQshjJVcocTYpDweupONsUi4AYLSobrRojNgEgg4cLQKAqpgYFB86jNKTJ6EsKQFfTw86vr4wfeftXjX5mpIhQlpHyRAhpMtkFFfh4JV0HIy4jZxSGcx01DBzoDVmDbSGhV7HlubgampQfvESSkNCwCkUsPrmawBAUXAw1D09oebs3KGf191QMkRI6ygZIoR0OblCib8S6+YWnZfmgQHwcTbBbG8b+IhNwOd17AgOx3FgjEFeVITkkaOA2loInZzqJl4/MxGqVlYd+nndASVDhLSOkiFCSLdyu7Cybm5RxG3klclgrqvWOLfIXLfjC7nKCwtReuIESkNCUXXtGgDAYuNG6Po+0+Gf1ZUoGSKkdZQMEUK6pVqFEn8l5GD/ldu4kPzfaNHkfhYY62IKLWHHl+WozchAye+/Q2/KFAiMjVFy7BhKQkKh6/sMtMY+Bb5W+wrUdkeUDBHSOkqGCCHd3u3CShy4ko4j1zKQXVoNoYCHMWIT+LpbYIzYBOqq/E753OIjvyL/u+9Qm5kJJhRCa4wPdH19oTVmTI+bfE3JECGto2SIENJjKJUcrqYXISQqE6Ex2cgvl0FDlY+xLqbwdTfHKJEx1FQ6NjHiOA5VkddRGhKC0j/+gMDcDPZHjgAAZMnJULW3B+N3TjLWkSgZIqR1lAwRQnokhZJD+M0ChERn4URsNgoraqAtFGCcxBS+HuYY7mh834KxD4qrrYU8NxcqlpZQVlRAOmw4+Do60Jk4ETq+vlBzlXTbESNKhghpHSVDhJAeT65Q4nJqAUKiM3EiNhul1XLoqqvgaVdT+LpbYKiDYYevX6SsqUH56dMoCQlF+YULQG0tVO3sYPb+GmgOGdKhn9URKBkipHWUDBFCnig1ciUupuQhJCoLf8bnoFwmh4GmKsb3NYOvuzkG9THs8Ef1FcXFKP3zT5SGhML0nXeg5ixC5bVIVMdEQ3vCBKiYmHTo5z0MSoYIaR0lQ4SQJ1Z1rQLnpXkIic7CXwk5qKxRwFhbiIl9zeDrYYEBNvrgdXBi1CD3669R8EMQwBg0Bg+Crq8vtMeNA19Hp1M+734oGSKkdZQMEUJ6haoaBc4k5iIkOhNnEnMhkythrquGiW7m8HU3Rz9rvQ6f7yO7cQOlIaEoCQlBbXo6VKyt4fDnSTDGwCmVYLyOvXXXFkqGCGkdJUOEkF6nXCbHXwk5OB6VhTBpHmoUSljpq+MZd3NMcreAq4VOhyZGHMehOjYW8txcaI8dC06hwA3fSVD38ICOry80Bw8CE3T8uklNUTJESOsoGSKE9GolVbU4FZ+DkOhMXEzOh1zJwc5QA77uFvD1MIezqXaHjxgpSkuR89l6lP35J5Tl5eAbGUFnwgTo+82GsE+fu9oqKyqQv307ivcfgKK4GHw9Pej5zYbRwoXgabZ/EUhKhghpHSVDhBBSr6iiBifjshESnYXLqflQcoCjiRZ83c3h624BRxOtDv08pUyG8vPnURoSivJz52C1ZTO0hg1DbVYWlJWVUDEzw83nZ6H29m1wMlnjcUwohIq1NfocDG53QkTJECGto2SIEEJakF8uwx+x2QiJysSVtEJwHOBkogUfsQlGOxvDy9agQ9cxUpSVgaeuDiYQIHfjRhT8tA18Q0MoiosBheKe9kwohMGiF2GyYkW7+qdkiJDWUTJECCH3kVNajd9jsvBXQi6u3CxEjUIJTVU+hjsZwcfZBKOdTWCmq9ZhnyfPy0PpHyeQs359i4lQA76BPkSXL7erT0qGCGkdJUOEEPIAKmRyXE4twNmkXJxLzEVmSTUAQGymXTdqJDJGf1t9qHTAIo8JLhKgrf+jGYNLQny7+qJkiJDWde7jC4QQ8oTRrC/5MU5iCo7jkJxbjrOJuTiXlIcfw25gy7lUaKsJMNLJGKOcjTFaZAwTnYcbNeLr6UFRVNT6fn29hzwLQkhTlAwRQshDYoxBZKoNkak2lo5yQFl1LS6l5ONcUh7OJuUiNCYLANDXUgejRSbwERujn7V+u1fA1vObjcKftt01ebrxs4VC6M2e3aHnQ0hvRbfJCCGkE3Ach4SsMpyT5uJcYh6uphdBoeSgq66CkSJj+DgbY6TIGEZawlb7UFZU0NNkhDwGlAwRQshjUFJVi4vJ+XVzjZLykF8uA2OAu6UuRjvXPaHmbqV3z6hR4zpDBw5AUVQMvr4e9GbTOkOEdCRKhggh5DFTKjnEZ5XibGIuzibl4vrtYig5wEBTFaNExhjtbIyRTsbQ11TtsM+kZIiQ1lEyRAghXayoogZhyXk4l5SH89I8FFbUgMeAftZ6GO1sAh9nE7ha6DxSUVlKhghpHSVDhBDSjSiUHGIySuqeUJPmIfpOMTgOMNIS4j1fFzzbz/Kh+qVkiJDW0dNkhBDSjfB5DP2s9dDPWg+vjhMhv1yGMGndqJHpQz6iTwhpGyVDhBDSjRlpCTGtvxWm9bfq6lAIeWJ1XGGdFjDGxjPGkhhjKYyx1S3sX8kYi2eMRTPG/mKM2XZmPIQQQgghzXVaMsQY4wP4HsAEABIAsxljkmbNIgF4cRznDuAwgM87Kx5CCCGEkJZ05siQN4AUjuNucBxXAyAYwLNNG3Acd5bjuMr6t/8AoHFgQgghhDxWnZkMWQK43eT9nfptrXkRwB8t7WCMLWGMRTDGIvLy8jowREIIIYT0dp06Z6i9GGNzAXgB2NDSfo7jtnIc58VxnJexsfHjDY4QQgghT7TOfJosA4B1k/dW9dvuwhh7CsA7AEZxHHdvNUJCCCGEkE7UmSND/wJwYoz1YYypApgF4LemDRhjngCCAEzmOC63E2MhhBBCCGlRpyVDHMfJASwHcBJAAoCfOY6LY4ytY4xNrm+2AYAWgEOMseuMsd9a6Y4QQgghpFP0uHIcjLE8ALe6Og4ARgDyuzqIR9CT4+/JsQMUf1fqybEDjxa/LcdxNOmSkBb0uGSou2CMRfTkOj89Of6eHDtA8Xelnhw70PPjJ6S76hZPkxFCCCGEdBVKhgghhBDSq1Ey9PC2dnUAj6gnx9+TYwco/q7Uk2MHen78hHRLNGeIEEIIIb0ajQwRQgghpFejZIgQQgghvRolQy1gjI1njCUxxlIYY6tb2C9kjB2s3x/OGLNrsu+t+u1JjLGnH2vgePjYGWN2jLGq+sUvrzPGfnjcsdfHcb/4RzLGrjHG5Iyx55rt82eMJdd/+T++qBs//1FiVzS59l2y+Gg74l/JGItnjEUzxv5ijNk22del174+hkeJv0uvfztif4kxFlMf30XGmKTJvi79P4eQJwLHcfTV5AsAH0AqAHsAqgCiAEiatVkG4If617MAHKx/LalvLwTQp74ffg+J3Q5AbA+49nYA3AHsBvBck+0GAG7U/6lf/1q/J8Rev6+8B1x7HwAa9a8DmnzvdOm1f9T4u/r6tzN2nSavJwM4Uf+6S//PoS/6elK+aGToXt4AUjiOu8FxXA2AYADPNmvzLIBd9a8PAxjLGGP124M5jpNxHHcTQEp9f4/Lo8TeHdw3fo7j0jiOiwagbHbs0wBOcRxXyHFcEYBTAMY/jqDrPUrs3UF74j/LcVxl/dt/UFd8Gej6aw88WvxdrT2xlzZ5qwmg4cmXrv4/h5AnAiVD97IEcLvJ+zv121psw9XVYCsBYNjOYzvTo8QOAH0YY5GMsfOMsRGdHWwLHuX69YRr3xY1xlgEY+wfxtiUDo2sfR40/hcB/PGQx3aGR4kf6Nrr367YGWMvM8ZSAXwOIPBBjiWEtE3Q1QGQbiMLgA3HcQWMsQEAjjLGXJv9Rko6jy3HcRmMMXsAZxhjMRzHpXZ1UC1hjM0F4AVgVFfH8jBaib/bX3+O474H8D1jzA/AuwC6ZG4WIU8iGhm6VwYA6ybvreq3tdiGMSYAoAugoJ3HdqaHjr1+mL0AADiOu4q6uQeiTo+4ldjqPcj16wnXvlUcx2XU/3kDwDkAnh0ZXDu0K37G2FMA3gEwmeM42YMc28keJf6uvv4Pev2CAUx5yGMJIS3p6klL3e0LdaNlN1A3GbFhMqNrszYv4+5JyD/Xv3bF3ZMZb+DxTqB+lNiNG2JF3UTODAAG3e3aN2m7E/dOoL6Jugm8+vWvH1v8jxi7PgBh/WsjAMloNoG2O8SPugQhFYBTs+1deu07IP4uvf7tjN2pyetJACLqX3fp/zn0RV9PyleXB9AdvwBMBCCt/4/znfpt61D32yQAqAE4hLrJilcA2Dc59p3645IATOgpsQOYDiAOwHUA1wBM6qbXfiDq5kVUoG40Lq7JsQvrzysFwIKeEjuAoQBi6n+oxQB4sZte+9MAcuq/R64D+K27XPtHib87XP92xP5Nk3+fZ9EkWerq/3Poi76ehC8qx0EIIYSQXo3mDBFCCCGkV6NkiBBCCCG9GiVDhBBCCOnVKBkihBBCSK9GyRAhhBBCejVKhghpBWPMsEkl82zGWEb963LG2Oaujo8QQkjHoEfrCWkHxtgHqKtsvrGrYyGEENKxaGSIkAfEGBvNGAupf/0BY2wXY+wCY+wWY2waY+xzxlgMY+wEY0ylvt2A+gK4VxljJxlj5l17FoQQQhpQMkTIo3MAMAbAZAB7AZzlOM4NQBWAZ+oTom9RV4JjAIDtAD7uqmAJIYTcjarWE/Lo/uA4rpYxFgOAD+BE/fYYAHYAnAH0BXCKMYb6NlldECchhJAWUDJEyKOTAQDHcUrGWC3330Q8Jer+jTHU1SEb0lUBEkIIaR3dJiOk8yUBMGaMDQEAxpgKY8y1i2MihBBSj5IhQjoZx3E1AJ4DsJ4xFoW6yuNDuzQoQgghjejRekIIIYT0ajQyRAghhJBejZIhQgghhPRqlAwRQgghpFejZIgQQgghvRolQ4QQQgjp1SgZIoQQQkivRskQIYQQQnq1/wckXZISNAaphwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fn = lambda t, y, params: -params[0]*y\n",
    "params = (5,)\n",
    "y0 = 1\n",
    "h = 0.1\n",
    "eps = 0.2*h\n",
    "t0,t1,t2,t3 = 0,h,h*2,h*3\n",
    "tru_t = np.linspace( t0, t0+3*h, 10 )\n",
    "tru_y = y0*np.exp(-params[0]*tru_t)\n",
    "\n",
    "# use explicit trapezoidal to get first two points\n",
    "yprime = fn(t0,y0,params);\n",
    "yhatFE = y0 + h*yprime;\n",
    "yprimehatFE = fn(t1,yhatFE,params);\n",
    "y1 = y0 + h/2*(yprime+yprimehatFE);\n",
    "\n",
    "# Compute f at i-1\n",
    "fim1 = fn(t1,y1,params);\n",
    "# Compute f at i-2\n",
    "fim2 = fn(t0,y0,params);\n",
    "y2 = y1 + h/2*(3*fim1-fim2);\n",
    "\n",
    "# Compute f at i+1\n",
    "fim0 = fn(t2,y2,params);\n",
    "y3 = y2 + h/2*(3*fim0-fim1);\n",
    "\n",
    "plt.plot( tru_t, tru_y.T, 'C0' );\n",
    "plt.plot( [t0,t1], [y0,y1], 'C1--' );\n",
    "plt.plot( [t1,t2], [y1,y2], 'C2--' );\n",
    "plt.plot( [t2,t3], [y2,y3], 'C3--' );\n",
    "plt.plot( [t0-eps,t0+eps], [y0-eps*fim2,y0+eps*fim2], 'C4', linewidth=3)\n",
    "plt.plot( [t1-eps,t1+eps], [y1-eps*fim1,y1+eps*fim1], 'C5', linewidth=3)\n",
    "plt.plot( [t2-eps,t2+eps], [y2-eps*fim0,y2+eps*fim0], 'C6', linewidth=3)\n",
    "\n",
    "plt.plot( t0, y0, 'C0.', markersize=16)\n",
    "plt.plot( t1, y1, 'C1.', markersize=16)\n",
    "plt.plot( t2, y2, 'C2.', markersize=16)\n",
    "plt.plot( t3, y3, 'C3.', markersize=16)\n",
    "\n",
    "plt.legend( bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0., \n",
    "           labels=['True Solution','Explicit Trapezoidal Step','First AB2 Step',\n",
    "                   'Second AB2 Step','Slope 1 for first AB2 Step',\n",
    "                   'Slope 2 for first AB2 step\\nSlope 1 for second AB2 step',\n",
    "                   'Slope 2 for second AB2 step']);\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( 'Y');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}