[EXTRA_02]_Exercício_EDO_Solução_Numérica_Euler_e_ODEINT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/andrebelem/ensino/blob/main/%5BEXTRA_02%5D_Exerc%C3%ADcio_EDO_Solu%C3%A7%C3%A3o_Num%C3%A9rica_Euler_e_ODEINT.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e97Yi7bo1yE9"
   },
   "source": [
    "![](http://observatoriooceanografico.org/wp-content/uploads/elementor/thumbs/logo-o2-2-p31vla9ib17i85olifd4l1xk0e5lktrkal6u0s4c6o.png)<br>\n",
    "# Notebook no Python para prática e resolução de EDO pelo método de Euler e comparação com um *solver* (`odeint`) <br>\n",
    "\n",
    "Criado por [Y.C.L.F. Paoli](mailto:ypaoli@id.uff.br),  [M.G.C. Silva](mailto:marinagaldez@id.uff.br) @ [Observatório Oceanográfico](http://www.observatoiriooceanografico.uff.br) 2021 (versão 1) para a monitoria da disciplina de **Modelagem e Otimização de Sistemas Agrícolas e Ambientais**<br>\n",
    "[![DOI](https://zenodo.org/badge/342919527.svg)](https://zenodo.org/badge/latestdoi/342919527)\n",
    "<br>\n",
    "Você pode usar o Python para resolver EDOs usando agora um solver. Talvez você precise relembrar os conceitos apresentados em aula e que estão no Guia Rápido (veja no canal do [Github](https://github.com/andrebelem/ensino)). <br>\n",
    "   \n",
    "Histórico:<br>\n",
    "Jun 2021: revisado e publicado no Github por [A.L.Belem](mailto:andrebelem@id.uff.br)<br>\n",
    "Ago 2022: adicionado novos comentários no Markdown e mais algumas funções   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0zJKgUDo2Wos"
   },
   "source": [
    "(relembrando) Uma solução de uma EDO é uma função $y = f(x)$  cujas derivadas satisfazem a equação. Mas cautela - não está garantido que uma solução exista, e caso exista, normalmente ela não é única. Em uma equação diferencial, os coeficientes podem ser constantes ou funções da variável independente. Lembra disso em cálculo ? O material na internet sobre isso é [bem farto](https://www.google.com/search?q=EDO+matem%C3%A1tica) (pode conferir). Mas vamos ao que interessa:<br>\n",
    "\n",
    "Imagine este exercício: **Calcule** a solução aproximada do seguinte Problema de Valor Inicial (P.V.I.) usando Python simples (note que você pode fazer isso diretamente no colab).\n",
    "\\begin{cases}\n",
    "y'= f(y,t) = y + 2 t       \\\\\n",
    "y(0) = 1  \n",
    "\\end{cases}\n",
    "\n",
    "Note que em EDO, a solução numérica pelo Método de Euler é \n",
    "$y_{k+1} = y_k + h . f(x_k,y_k)$. Para essa equação, a solução exata é $y = 3. e^t - 2(t + 1)$<br>\n",
    "<br>\n",
    "\n",
    "Além da solução, neste exercício vamos plotar os resultados e usar algums pacotes acessórios, o `numpy` e o `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "P6wcFHnzeqEk"
   },
   "outputs": [],
   "source": [
    "# primeiro, importe os pacotes\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "dIWvos7qeqak"
   },
   "outputs": [],
   "source": [
    "# Vamos definir a equação em uma função EDO\n",
    "\n",
    "def f(y,t):\n",
    "    dy = y + 2*t\n",
    "    return dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "VZ706H4xffAQ"
   },
   "outputs": [],
   "source": [
    "# definição dos parâmetros do método de Euler\n",
    "# definimos aqui o passo de tempo. => teste diferentes passos (0.1, 0.01, ...)\n",
    "h = 0.1 \n",
    "tmax = 2 # t máximo igual a 2 significa que o cálculo será feito entre t=0 e t=2, isto é, tmax = último valor de t\n",
    "# essa função gera vários valores de t entre 0 e tmax, separados por um passo=h, que serão colocados no array t\n",
    "t = np.arange(0, tmax+h, h) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-8ZucXvbgCBc"
   },
   "source": [
    "Obs: O que chamamos de *passo* é a distância entre dois valores consecutivos da variável $x$ (neste caso, $t$), isto é, \n",
    "$h = t_{n+1} - t_{n}$<br>\n",
    "Note que não estamos definindo uma unidade de tempo (segundos ou horas ou dias) porque isso é intrínsico na sua equação (ou seja, faz parte dela. É a equação que define o tempo). Veja o gráfico abaixo<br>\n",
    "<center><img src=\"https://drive.google.com/uc?id=1ZtuexbeNkz_scTGkqgZeFtghcnrfuvQl\" width=\"400\"></img></center>\n",
    "<br>\n",
    "\n",
    "Imagine agora que eu vamos dizer qual o ponto de partida (ou seja, a condição inicial)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "NYa_5CbDiolI"
   },
   "outputs": [],
   "source": [
    "# condição inicial\n",
    "y0 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "ojfMhsP-i7_0"
   },
   "outputs": [],
   "source": [
    "# Esta é a Solução usando o método de Euler\n",
    "\n",
    "y_euler = np.zeros_like(t) # essa função cria um vetor composto por zeros do mesmo tamanho do vetor t\n",
    "y_euler[0] = y0 # define que o primeiro valor do vetor será igual à condição inicial, isto é, o y0.\n",
    "for n in range (t.size-1): # t.size dá o tamanho do array t\n",
    "    y_euler[n+1] = y_euler[n] + h * f(y_euler[n], t[n])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "TlqUhgHMlOfO"
   },
   "outputs": [],
   "source": [
    "# Aqui temos a solução exata para este exercício\n",
    "yexato = 3 * np.exp(t) - 2 * (t + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "id": "1-kauQpEkb91",
    "outputId": "b911094c-3b55-454b-e7c3-fd211793232b"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsqUlEQVR4nO3dfZzVc/7/8cdLpdSEyJJC+ZEoo8kIuagscn25u2nXbonNRSS2tVqk7LbschNy2S5y0cYquYhZiWzlYnynGlIkS1bEZFCmdDHT6/fH+0ymaaY5M3PO+ZyZ87zfbnPrnM/5XLx8zvF5fT7vS3N3REQkc20XdQAiIhItJQIRkQynRCAikuGUCEREMpwSgYhIhmsadQB10bZtW+/YsWPUYYiINCjz5s372t13q7y8QSaCjh07UlBQEHUYIiINipl9WtVyFQ2JiGQ4JQIRkQynRCAikuEaZB1BVTZu3Mjy5ctZt25d1KFIArRo0YIOHTrQrFmzqEMRafQaTSJYvnw5rVu3pmPHjphZ1OFIPbg7xcXFLF++nE6dOkUdjkij12iKhtatW8euu+6qJNAImBm77rqrnu5EUqTRJAJASaAR0XcpkjopSwRm9pCZFZnZe5WWX2FmH5jZIjP7W6riERFpSH74AYYNgy++SPy+U/lEMBE4qeICM+sLnAkc4u5dgdtSGE/CjR07lq5du5KdnU337t3Jz8/f5vp9+vSpsWPc6NGjad++Pd27d9/8991339U6tsLCQl588cVabyci6eHmm2H8ePjgg8TvO2WVxe4+28w6Vlp8KXCLu6+PrVOUqngS7c0332T69OnMnz+f5s2b8/XXX7Nhw4aE7Puqq65ixIgR9dpHYWEhBQUFnHLKKQmJSURS58MP4a9/hV/9Co47LvH7j7qOoDNwjJnlm9l/zOyw6lY0syFmVmBmBStXrkxhiPFZsWIFbdu2pXnz5gC0bduWPffcE4BXXnmFnJwcDj74YAYPHsz69evrfbxx48YxePBgABYuXEi3bt1Yu3Ytb7/9NkceeSQ5OTn06tWLJUuWsGHDBkaNGsWTTz5J9+7defLJJ6tcT0TSjzsMHQo77AC3JanMJOrmo02BXYAjgMOAf5nZvl7F/JnuPgGYAJCbm7vN+TWHD4fCwsQG2r073HFH9Z+feOKJ3HTTTXTu3Jnjjz+e/v3707t3b9atW8egQYN45ZVX6Ny5M7/5zW+47777GD58eNzHHjduHI8//jgAbdq0YdasWVx55ZX06dOHadOmMXbsWB544AFatmxJly5dmDNnDk2bNmXmzJn88Y9/ZOrUqdx0000UFBRw9913A7B69eoq1xOR9PKvf8HMmXD33bDHHsk5RtSJYDnwdOzC/7aZbQLaAul3y1+DrKws5s2bx5w5c5g1axb9+/fnlltuIScnh06dOtG5c2cABg4cyD333FOrRFBV0dB2223HxIkTyc7O5uKLL+aoo44CYNWqVQwcOJClS5diZmzcuLHKfca7nohEZ/VquOoqOPRQuOSS5B0n6kTwDNAXmGVmnYHtga/ru9Nt3bknU5MmTejTpw99+vTh4IMP5pFHHiEnJydpx1u6dClZWVl8UaEZwQ033EDfvn2ZNm0ay5Yto0+fPlVuG+96IhKdUaPgyy/h2WehSZPkHSeVzUcnA28CB5jZcjO7EHgI2DfWpPQJYGBVxUINwZIlS1i6dOnm94WFheyzzz4ccMABLFu2jI8++giAxx57jN69e9f7eKtWrWLYsGHMnj2b4uJipkyZsnl5+/btAZg4ceLm9Vu3bs3333+/xfZVrSci6WHBgtBK6JJL4LBqa08TI2WJwN0HuHs7d2/m7h3c/UF33+Du57t7N3fv4e6vpiqeRCspKWHgwIEcdNBBZGdns3jxYkaPHk2LFi14+OGH+fnPf87BBx/MdtttxyVVPONddNFF1TYlHTdu3BbNR5ctW8ZVV13F0KFD6dy5Mw8++CDXXnstRUVFXHPNNYwcOZKcnBxKS0s376Nv374sXrx4c2VxdeuJSPQ2bYJLL4W2bWHs2OQfzxriDXhubq5Xvmi+//77HHjggRFFJMmg71Qy1YQJcPHF8Oij8OtfJ26/ZjbP3XMrL4+6+aiIiFRQVATXXgu9e8P556fmmEoEIiJp5A9/gO+/h3vvhVQNuaVEICKSJubMgYkTYcQIOOig1B1XiUBEJA1s3BgqiPfZB66/PrXHjrofgYiIEPo/LVoU+gy0apXaY+uJICIPPPAA3377bdRhiEga+OwzGD0aTj8dzjgj9cdXIkggM+P8CtX8paWl7Lbbbpx22mlbrHfTTTfRpk0b2rRpA9R9iOiOHTvy9dd164g9ceJELr/88rjXf+2119hpp5226M8wc+bMbW4zaNCgzR3dRKR6V14ZBpe7665ojp+5RUNlZZCXF7rv5eTAySfXuw93q1ateO+99/jhhx/YYYcdePnllzf33q1o1KhRW7xvKENEH3PMMUyfPj1p+y8tLaVp08z9SUpmeuEFmDYN/vIX6Ngxmhgy84mgrAz69YMBA+DGG8O//fqF5fV0yimn8MILLwAwefJkBgwYsPmzNWvWMHjwYHr27ElOTg7PPvtslUNEf/PNN5x11llkZ2dzxBFH8O677wJQXFzMiSeeSNeuXbnooouo2Bnw9ttvp1u3bnTr1o07qhls6eGHH6Zz58707NmT119/ffPylStXcu6553LYYYdx2GGHbfFZTZYtW0a3bt02v7/tttsYPXr0VuvNmzeP3r17c+ihh9KvXz9WrFgBhMl5hg8fTm5uLnfeeWfcxxVpDNauhSuugC5d4He/izAQd29wf4ceeqhXtnjx4q2WVev5592zstzD01j4y8oKy+uhVatW/s477/i5557rP/zwgx9yyCE+a9YsP/XUU93dfeTIkf7YY4+5u/u3337r+++/v5eUlPjDDz/sQ4cO3byfyy+/3EePHu3u7q+88oofcsgh7u5+xRVX+JgxY9zdffr06Q74ypUrvaCgwLt16+YlJSX+/fff+0EHHeTz58/fIrYvvvjC99prLy8qKvL169d7r169Nh9zwIABPmfOHHd3//TTT71Lly5b/bfNmjXLd9xxRz/kkEM2/3300Uf+ySefeNeuXTevd+utt/qNN97o7u4DBw70p556yjds2OBHHnmkFxUVubv7E0884RdccIG7u/fu3dsvvfTSKs9nrb5TkQbouuvC5efVV1NzPKDAq7imZuZz+IIFsGbNlsvWrAmTGFQqz6+t7Oxsli1bxuTJk7cq6pkxYwbPPfcct8Vml1i3bh3/+9//ttrH3LlzN88NcNxxx1FcXMzq1auZPXs2Tz/9NACnnnrq5jqGuXPncvbZZ9Mq1tTgnHPOYc6cOVuMfJqfn0+fPn3YbbfdAOjfvz8ffvghADNnzmTx4sWb1129ejUlJSVkZWVtEVdVRUPLli2r8ZwsWbKE9957jxNOOAGAsrIy2rVrt/nz/v3717gPkcZmyRL4299C7+G+faONJTMTQU5OaJ9VUvLjslatwuwzCXDGGWcwYsQIXnvtNYqLizcvd3emTp3KAQccsMX6Nc1tnGybNm3irbfeokWLFrXetmnTpmzatGnz+3Xr1m21jrvTtWtX3nzzzSr30SrVbeVEIuYOl10GLVsmb9ax2sjMOoKTT4bDD4esrNCHOysrvD/55ITsfvDgwdx4440cfPDBWyzv168f48eP31y2v2DBAmDrIaKPOeYYJk2aBITWOm3btmXHHXfk2GOP5Z///CcAeXl5m5ufHnPMMTzzzDOsXbuWNWvWMG3aNI455pgtjn344Yfzn//8h+LiYjZu3MhTTz21+bMTTzyR8ePHb35fWIvp3XbffXeKioooLi5m/fr1VVYmH3DAAaxcuXJzIti4cSOLFi2K+xgijc0TT8Crr4YK4t13jzqaTH0iaNIEXnoptBoqLAxPAgloNVSuQ4cODBs2bKvlN9xwA8OHDyc7O5tNmzbRqVMnpk+fTt++fbnlllvo3r07I0eOZPTo0QwePJjs7GxatmzJI488AsCNN97IgAED6Nq1K7169WLvvfcGoEePHgwaNIiePXsCYUjryhPitGvXjtGjR3PkkUey8847073C089dd93F0KFDyc7OprS0lGOPPZb7779/q/jnzJmzxXbXX389P/vZzxg1ahQ9e/akffv2dOnSZavttt9+e6ZMmcKwYcNYtWoVpaWlDB8+nK5du9b63Io0dKtWwdVXQ25uGGE0HWgYaklb+k6lMRo2LMw//PbbIRnELQFN3qsbhjoznwhERCIwbx7cc0+oH6h1EujXD/LzQ8OWVq1CcfZLLyWkJCOVU1U+ZGZFsWkpK3/2OzNzM2ubqnhERFJp/XoYNAh+8hP4859ruXFeXkgCJSWhprmkJLzPy0tIbKmsLJ4InFR5oZntBZwIbN2OspYaYjGXVE3fpTQ2o0bBe+/Bgw/CzjvXcuNtNXlPgFTOWTwb+KaKj8YB1wD1+j+/RYsWFBcX6wLSCLg7xcXFdWrOKpKOXn8dbr0VLroI6jSSTHmT94oS2OQ90joCMzsT+Nzd37EapuIxsyHAEGBza5mKOnTowPLly1m5cmUyQpUUa9GiBR06dIg6DJF6KymBgQPDPAO3317HnZQ3ea9cR5CgJu+RJQIzawn8kVAsVCN3nwBMgNBqqPLnzZo1o1OnTgmNUUSkvq65Bj7+GGbNgtat67iTJDd5j/KJ4P8BnYDyp4EOwHwz6+nuX0YYl4hIQsyYAffdB1ddFSajr5cmTcIQOPUcBqcqkSUCd18I/KT8vZktA3LdvW4D7IuIpJHvvoPBg+HAA2Hs2Kij2bZUNh+dDLwJHGBmy83swlQdW0Qk1YYNgy+/hEcfhR12iDqabUvZE4G7D6jh844pCkVEJKmmTYPHHgtNRmvVcSwimTnonIhIkhQVhTGEevSA66+POpr4KBGIiCSIe0gCq1eHIqFmzaKOKD4aa0hEJEEeewyeeSZ0HmtIg+vqiUBEJAE++yzMP3zMMaG5aEOiRCAiUk+bNoWmomVlMHFiwvp5pYyKhkRE6um++2DmTLj/fth336ijqT09EYiI1MPSpfD738NJJ8GQIVFHUzdKBCIidVRWFgaUa94c/vGPMAV6Q6SiIRGROrr1VnjzTZg0Cdq3jzqautMTgYhIHbz7bug5/LOfwYBtjpuQ/vREICJSSxs2wG9+A7vsEiqK4yoSSsDk88miRCAiUks33QTvvAPPPQdt45lpPcmTz9eXioZERGohPx9uvhkuuABOPz3OjZI8+Xx9KRGIiMTp229DfUCHDjBuXC02TPLk8/WloiERkTi4w6BBYSiJOXNgp51qsXH55PMlJT8uS+Dk8/WlJwIRkTjcdluoE7jtNjjiiFpuXD75fFZWqFnOykro5PP1pScCEZEazJkDI0eGpqLDhtVhB0mefL6+zN1TcyCzh4DTgCJ37xZbditwOrAB+C9wgbt/V9O+cnNzvaCgIInRiogEX30VSnaysqCgAHbcMeqI6s7M5rn7VnOmpbJoaCJwUqVlLwPd3D0b+BAYmcJ4RES2qawsVA5/9x1Mndqwk8C2pCwRuPts4JtKy2a4e2ns7VtAh1TFIyJSkxtvhFmz4N574eCDo44medKpsngwkB6NakUk4+XlwdixcOGFobVQY5YWicDMrgNKgUnbWGeImRWYWcHKlStTF5yIZJxPP4Xzz4dDDoHx46OOJvkiTwRmNohQifwr30bNtbtPcPdcd8/dbbfdUhafiGSW9evhF7+A0lKYMgV22CHqiJIv0uajZnYScA3Q293XRhmLiAjAiBHw9tuhcni//aKOJjVS9kRgZpOBN4EDzGy5mV0I3A20Bl42s0Izuz9V8YiIVPbkk3D33XD11XDOOVFHkzopeyJw96pG7H4wVccXEdmWDz6Aiy6CXr3glluijia1Iq8jEBGJ2po1odfwDjvAv/4FzZpFHVFqaYgJEclo7nDJJbB4McyY0bCnnKwrJQIRyWh//zs8/jiMGQPHHx91NNFQ0ZCIZKz588Mgcv36wfXXRx1NdPREICIZ6dtvQ73AbruFJ4Lt4rktTuN5h+tDiUBEMk7FSWZmz24c8w7Xh4qGRCTjVJxk5sgj49wozecdrg8lAhHJKNOnw7XX1mGSmTSfd7g+lAhEJGMsWADnnReK9ydODLNGxq183uGK0mje4fpQIhCRjPD553DaabDLLvD881tf02uU5vMO14cqi0Wk0SspCUng++9h7lxo164OO0nzeYfrQ4lARBq18ukmFy4M9QPZ2fXYWZMmIaOcdlrC4ksHSgQi0qhdfXVIAPfeCydVnjVdANURiEgjNn483HUXXHUVXHpp1NGkLyUCEWmUXngBhg+HM8+EW2+NOpr0pkQgIo1OYSH07x9afE6a1Cjqc5NKiUBEGpXyZqJt2oTew7VuJpqBUjlV5UNmVmRm71VYtouZvWxmS2P/tklVPCLS+JSUwOmnw6pVoWhozz2jjqhhSOUTwUSgcp39tcAr7r4/8ErsvYhIrZWVwS9/Ce+8E2YZq1cz0QyTskTg7rOBbyotPhN4JPb6EeCsVMUjIo3L734XegyPH98oOvumVNR1BLu7+4rY6y+B3aMMRkQaprvvhjvvDM1EL7ss6mganrTpUObubmZe3edmNgQYArD33nunLC4RSW8vvABXXglnnBFnM9FGOrlMfUSdCL4ys3buvsLM2gFF1a3o7hOACQC5ubnVJgwRyRyFhWE00e7d4Z//jON63ognl6mPqIuGngMGxl4PBJ6NMBYRaUDKm4nuvHMtRhNtxJPL1Ecqm49OBt4EDjCz5WZ2IXALcIKZLQWOj70XEdmm1avr2Ey0EU8uUx8pKxpy9wHVfPTTVMUgIg3fmjVw6qlhNNHnnqtlM9HyyWVKSn5c1kgml6mPqIuGRETi9sMPoVL4jTdCnUCtm4k24sll6iPqymIRkbhs2BDmGZ41Cx55BH7+8zrspBFPLlMfSgQikvZKS8PkMi++CA88AL/+dT121kgnl6kPFQ2JSForK4OBA+Hpp+GOO2DIkKgjanyUCEQkbW3aBBdfHOoDbr45dByTxFMiEJG05B4u/A8+CDfcANdqSMqkqXUiMLNWZpbZNSsiklTu8Ic/hDGERoyAMWOijqhxqzERmNl2ZvZLM3vBzIqAD4AVZrbYzG41s/2SH6aIZJIxY8K4QZddBn/7W2jpKckTzxPBLOD/EeYK2MPd93L3nwBHA28BfzWz85MYo4hkkL/+NSSCCy4IQ0orCSRfPM1Hj3f3jWY23917lC9092/M7HN3P9fMmiUxRhHJEHfdFeoCBgyAv/8dtlMtZkrEkwjONrMeQGszOxBY4u6bYp9NALLdfWPSIhSRjPD3v4fK4bPPDh3GttnHS0NJJ1Q8ieB1oAVwEXA7YdC474AvgB+SF5qIZIrHHw/NRE8+GSZPhmbbKmPQUNIJV2MicPfPgUfN7L/u/jqAme0KdCRUHIuI1NmUKaHDWJ8+MHUqNG9ewwYVh5KGLYeSVm/hOomn1ZABlCeB2Otid5/n7msqriMiUhvTp4f6gCOOCCOJ7rBDHBtpKOmEi6vVkJldYWZbzA9pZtub2XFm9gg/Ti4jIhKXF16Ac88N4769+GIYCDQu5UNJV6ShpOslnkRwElAGTDazL2L9Bz4GlgIDgDvcfWISYxSRRmbSJDjzTOjWLRTt77RTLTbWUNIJF08dwTrgXuDeWDPRtsAP7v5dkmMTkUZo/HgYNgz69oVnnoEdd6zlDjSUdMLFPQy1mb0MjHD3d5IYj4g0Uu6ho9iYMeFp4IknoEWLOu5MQ0knVG26a/wBuMPMHjazdokMwsyuMrNFZvaemU02s7r+PEQkDW3aFJ4CxoyBQYNCS6E6JwFJuLgTgbvPd/e+wHTg32Z2o5nFU8e/TWbWHhgG5Lp7N6AJcF599ysi6WHjxjCRzN13w9VXh9FEm2pKrLRSqw7csWaiS4D7gCuApWZWn7mCyjUFdjCzpkBLQmc1EWng1q6Fs84K8wn85S9w220aNiIdxf2VmNnrwOfAOKA9MAjoA/Q0swl1DSDWYe024H/ACmCVu8+o4vhDzKzAzApWrlxZ18OJSIp8913oAJyXF6aXHDlSA8ilK3P3+FY06wos9io2MLP33f3AOgVg1gaYCvQHvgOeAqa4++PVbZObm+sFBQV1OZyIpMCXX4Yk8P77oalonSaal4Qzs3nunlt5eW3qCBZVlQRiTq1zZHA88Im7r4wNXvc00Kse+xORCH3yCRx9NHz0Ueg5rCSQ/hJSZePuH9dj8/8BR5hZS8Igdj8FdLsv0gAtXBieBNatg1deCUNHVEmjh6aVyOvu3T3fzKYA84FSYAFheGsRaUDefBNOOQVatoQ5c6Br12pW1OihaSct6u/d/UZ37+Lu3dz91+6+PuqYRCR+L70Exx8PbdvC669vIwnAlqOHum85eqhEIi0SgYg0XE8+CaefDp07w9y50LFjDRto9NC0o0QgInXiDnfe+eMw0q+9BrvvHseGGj007SgRiEitrV8Pv/0tDB8OZ5xRyxFENXpo2om8slhEGpaiIjjnnFAXcN11cNNNtewtrNFD044SgYjErbAwPAF8/XUYPbR//zruSKOHphUVDYlIXJ56Co46KtQNzJ1bjyQgaUeJQES2adMmGDUKfvGLUIpTUAA9ekQdlSSSioZEpFolJWEI6WeegQsugPvug+bNo45KEk2JQESq9MknYSaxRYvgjjvCxDIaPbRxUiIQka289hr87Gc/Dgl04olVrKTxghoNJQIR2cL998MVV8B++8Fzz8H++1exksYLalRUWSwiQJhS8rLL4NJLwxPAW29VkwRA4wU1MkoEIsLXX4eL/333wTXXhCeBbfYU1nhBjYqKhkQy3MKFoZPYihXw2GNw/vlxbFQ+XlBJyY/LNF5Qg6UnApEM5R7qA3r2DGMHzZ4dZxIAjRfUyOiJQCQDffNNGDTu6afhhBPg0Udhjz1qsQONF9SoKBGIZJjZs+FXvwoTzN96K1x9dS0HjSun8YIajbQoGjKznc1sipl9YGbvm9mRUcck0tiUlsLo0dC3b+gd/MYbMGJEHZOANCrp8kRwJ/Bvd/+ZmW0PtIw6IJHG5H//C08Bc+eGISPuuQdat446KkkXkScCM9sJOBYYBODuG4ANUcYk0phMnQoXXRSeCOJuFSQZJR0eCjsBK4GHzWyBmf3DzFpVXsnMhphZgZkVrFy5MvVRijQwa9fCxReHoSL23z80/d8qCZSVwfTp8Kc/hX/LyiKJVaJl7h5tAGa5wFvAUe6eb2Z3Aqvd/YbqtsnNzfWCgoKUxSjS0Lz7bphLePHi0EHsT3+C7bevtJKGicg4ZjbP3XMrL0+HJ4LlwHJ3z4+9nwJotHOROnAP5f89e0JxMcyYAX/9axVJADRMhGwWeSJw9y+Bz8zsgNiinwKLIwxJpEEqLoazzoLLL4fjjgtPBSecsI0NNEyExEReWRxzBTAp1mLoY+CCiOMRaVBeey2U/xcVwbhxYe6AGpuFapgIiUmLRODuhcBW5VYism2rVsHIkWGwuM6dQ8lOTk6cG5cPE1G5jkDDRGSctEgEIlJ7zz4LQ4eGweKuuipUCLfaqr3dNmiYCIlRIhBpYL78MkwcM2UKZGfDtGlw2GF13JmGiRDSoLJYROLjDg8+CAceCM8/D2PHQkFBPZKASIyeCEQagI8+giFDYNYsOPZYmDABDjig5u1E4qFEIJLGNm6E228Pg8U1bx4SwIUXVmoRpEnkpZ6UCETS1Lx5YYygwkI45xwYPx723LPSSuodLAmgOgKRNLN2Lfz+96F38FdfhUHjpk6tIgmAegdLQigRiKSRmTPh4IPhttvC08DixeFpoFrqHSwJoEQgkga++goGDQpDQjRpEnoKP/AA7LxzDRuW9w6uSL2DpZaUCEQitG4d3Hwz7LcfTJoEf/xjGCOod+84d6BJ5CUBVFksEgF3ePJJuPZa+PRTOPNM+NvfwjARtaLewZIASgQiKfbmm2HC+LfeCtfthx8O8wjXmXoHSz2paEgkRZYtg/POg169wlPAQw+FnsF9+6KZwiRSeiIQSbLVq0M9wLhxoSPYqFGheWhWVmwF9QWQiOmJQCRJSktDy5/99oNbboFf/AI+/BDGjKmQBEB9ASRySgQiSTBjRmjZeckl0KUL/N//waOPQocOVaysvgASMSUCkQRavBhOOSWU9KxdG3oE/+c/kLutaZfUF0AiljaJwMyamNkCM5sedSwitfXRRzB4cJgf4I03Qs/g8l7BZjVsrL4AErF0qiy+Engf2DHqQETitWRJmBdg0iTYfvswcfz110PbtrXYifoCSMTSIhGYWQfgVGAscHXE4YjUaPFi+POf4YknoEWLMFXkiBGwxx513KH6AkiE0iIRAHcA1wCtq1vBzIYAQwD23nvv1EQlUsnChSEBPPUUtGwZmoH+7nfwk101J4A0XJEnAjM7DShy93lm1qe69dx9AjABIDc311MTnUhQWBj6ej39NLRuDSNHhqeAtm1RPwBp8NKhsvgo4AwzWwY8ARxnZo9HG5JIMG8enHVWuMmfORNuuCH0EB47tkI9gPoBSAMXeSJw95Hu3sHdOwLnAa+6+/kRhyUZ7u23Q3F9bm5o/jlmTBgW4qabYJddKq2sfgDSwEVeNCSSLtxhzpwwHMS//x0u+H/+c2gJtNNO29iwvB9AScmPy9QPQBqQtEoE7v4a8FrEYUiG+eEHmDwZ7roL3nknFPnccgtcdlmoD6hReT+AynUE6gcgDURaJQKRVFq+HO69FyZMgOJiyO5axr8vz6NvmwVs3zUHWp4MxFHZq34A0sApEUhGcQ89f++6Kwz/4A5nnAHDhpbR5+Z+2MQ6tvxRPwBpwCKvLBZJhXXr4JFHQuXv0UeHQeGuuioMDTFtGvRdl4e9rZY/kpmUCKRR++KL0ORz773D5PDr1sH994dioVtvhU6dYiuq5Y9kMBUNSaPjHm7m77or9AAuK4PTT4dhw+C446oZBE4tfySDKRFIo/H116H1z8SJMH9+aPI5bBgMHQr77lvDxmr5IxlMiUAatI0b4cUXQ/n/9OlQtrGMy/fN497TFpD9mxx2OCfO1jtq+SMZzNwb3rA9ubm5XlBQEHUYEqHCwnDn/89/wsqVsPvu8OtfljHqjX60XqQxf0SqYmbz3H2raZJUWSwNxldfhQngDzkkFOnfdx/07h2eBJYvh1uPywtJQC1/RGpFRUOS1tavDxf6iRPD9bysDHr2hHvugfPOqzTuz7Za/qh9v0i1lAgk7bj/ONn75MnwzTew555h4peBA+HAA6vZUC1/ROpEiUDSQllZ6PE7dSo8+3QZXT/L47CmC/jDkTkccu3JHN+vSc3F/Gr5I1InSgQSmQ0bYNasMNnLM89AURHssH0Zb+zYj67N82m6YQ22oBXcfjj0e4kax/1Ryx+ROlEikJRauzYM7zB1Kjz/PKxaFW7cTz0Vzj0XTiOPlhfmw/pY8U7FCt94yvk15o9IrSkRSNKtWgUvvBDu/PPyQjJo0wbOPhvOOQdOOCFMAA/An1ThK5JqSgSSeGVlFD+exyfTFjB9eQ5/ffdk1m1swh57hMrec84JzT6bNatiW1X4iqScEoEkxNq1YXavGXll/Pwf/ei6Jp8c1nCgteKXHQ7n60kvccRRTdiupp4rqvAVSbnIE4GZ7QU8CuwOODDB3e+MNiqpiTssXBjqZmfMCElg/Xo4q2ke2Z5PS8IdfSsvofO3+XRelQfbxVnGrwpfkZSKPBEApcDv3H2+mbUG5pnZy+6+OOrAZEtffQUvvxwu/C+/DF9+GZZ37RqmdTzxRDjujQVs/+d6lvGrwlckpSJPBO6+AlgRe/29mb0PtAeUCKJUVsb6Z/JYPn0Bs1fncPd/T2b+O+GuvG3bUMF74onh3/btK2xXqjJ+kYYm8kRQkZl1BHKA/Co+GwIMAdh7771TG1iGKCqC11+HN+aU8YsH+9FldT4dWcNPaEX2zofz8p9f4oSTmpCTQ/Vl/SrjF2lw0mb0UTPLAv4DjHX3p7e1rkYfrT93WLIE5s4NF//XX4elS8NnZzWdzqRNA2i5qcJdfVZWGO8hnuKasjKV8YukoepGH02LJwIzawZMBSbVlASkFsovyAsWsLFbDm/vejKvv9WEuXPDcA7FxWG1XXcN8/j+9rdw1FHQ86UFNP1TPcr5VcYv0qBEngjMzIAHgffd/fao42kMyspgyeIy2gzoxy4f5tNs4xrW0YofOJyRvMR+nZtwxhnh4n/UUdC5c6XpG79ROb9IJok8EQBHAb8GFppZYWzZH939xehCajjWr4dFi8LUjAsWhH/ffRf6rs1jMvk0jzXjbE0JfVrkU/z3PHY+v4Y7dZXzi2SUyBOBu88FqppOXGCL4p0fuuQwf/eTWfBuk80X/kWLwnSNAK1bh465v/0tDPxsAVnT1oSeGTFN169h52WFQA2JQG35RTJK5IlAtrZuHXz4IXywqIzc6/rR7rN8ti9dQ2mseOdKXmKXtk3o0QNOOilc/Hv0CBO0b27NMz0HZtSjeEfl/CIZQ4kgFSrc1ZOTs/nu+ttv4f33w98HH/z4+pNPQqueU8njZPLZoVLxzlf35bHrwNO2LNevTMU7IhInJYIkW7+2jNLj+9F8QT7brVvDhqatWJR1OKc1e4kvV/5Y1NK8eai0zc2F888Ps3D1nrOArPu2Lt5pu7wQTMU7IpIYSgTxqOaOHqC0FD77DJYtC3fy5f+Wv875PI9/kk/T2F19i9ISDlydz7XH5bGx32l06RIu+h07VnGNzsqBR1W8IyLJlTmJYBsX8+q4w6pvymh2Wj+av5NPk3Vr2NCsFf/d9XCu6PwS/13WhOXLw67LbbcddOgAnTqF4Rd+8+kCsl7b8q6+pa/hyt6FMEKtd0QkepmRCMrKoF+/zRdUb9WKjTmH8+H4l1hR1IQvvoAVK8Jf5dc/XReaYZbf0TffUMLeX+bTY6c82h19Gh07hot+p07hrn6vvSqNsz89B/6vjnf1Kt4RkRTIjESQlxeSQOxibCUlrJ+Tz7Xd83ihQlPK1q2hXTvYc89w473nnnD2ewvImrHlHX0Wa7j1V4VwfRxFLvW9q1fxjogkWWYkggVbT3/YijXc/ItCrhl6Gu3ahQSQlVXFttNz4PV6ltPrrl5E0lhmJIIqpj/cLqsVB/+6Oxxbw7aJKKfXXb2IpLHMSAT1uZjrjl5EGrm0GYa6Nuo0DLWGRhaRDJfWw1CnhIpnRESqVN08UyIikiGUCEREMpwSgYhIhlMiEBHJcEoEIiIZrkE2HzWzlcCnddy8LfB1AsNJFMVVO4qrdhRX7aRrXFC/2PZx990qL2yQiaA+zKygqna0UVNctaO4akdx1U66xgXJiU1FQyIiGU6JQEQkw2ViIpgQdQDVUFy1o7hqR3HVTrrGBUmILePqCEREZEuZ+EQgIiIVKBGIiGS4RpUIzOwkM1tiZh+Z2bVVfN7czJ6MfZ5vZh0rfDYytnyJmfVLcVxXm9liM3vXzF4xs30qfFZmZoWxv+dSHNcgM1tZ4fgXVfhsoJktjf0NTHFc4yrE9KGZfVfhs6ScLzN7yMyKzOy9aj43M7srFvO7ZtajwmfJPFc1xfWrWDwLzewNMzukwmfLYssLzayW47rXO64+Zraqwnc1qsJn2/z+kxzX7yvE9F7s97RL7LNknq+9zGxW7DqwyMyurGKd5P3G3L1R/AFNgP8C+wLbA+8AB1Va5zLg/tjr84AnY68Piq3fHOgU20+TFMbVF2gZe31peVyx9yURnq9BwN1VbLsL8HHs3zax121SFVel9a8AHkrB+ToW6AG8V83npwB5gAFHAPnJPldxxtWr/HjAyeVxxd4vA9pGdL76ANPr+/0nOq5K654OvJqi89UO6BF73Rr4sIr/H5P2G2tMTwQ9gY/c/WN33wA8AZxZaZ0zgUdir6cAPzUziy1/wt3Xu/snwEex/aUkLnef5e5rY2/fAjok6Nj1imsb+gEvu/s37v4t8DJwUkRxDQAmJ+jY1XL32cA321jlTOBRD94CdjazdiT3XNUYl7u/ETsupO63Fc/5qk59fpeJjislvy0Ad1/h7vNjr78H3gfaV1otab+xxpQI2gOfVXi/nK1P5OZ13L0UWAXsGue2yYyrogsJWb9cCzMrMLO3zOysBMVUm7jOjT2GTjGzvWq5bTLjIlaE1gl4tcLiZJ2vmlQXdzLPVW1V/m05MMPM5pnZkAjiOdLM3jGzPDPrGluWFufLzFoSLqZTKyxOyfmyUGSdA+RX+ihpv7HMmaGsATCz84FcoHeFxfu4++dmti/wqpktdPf/piik54HJ7r7ezC4mPE0dl6Jjx+M8YIq7l1VYFuX5Sltm1peQCI6usPjo2Ln6CfCymX0Qu2NOhfmE76rEzE4BngH2T9Gx43E68Lq7V3x6SPr5MrMsQvIZ7u6rE7nvbWlMTwSfA3tVeN8htqzKdcysKbATUBzntsmMCzM7HrgOOMPd15cvd/fPY/9+DLxGuFNISVzuXlwhln8Ah8a7bTLjquA8Kj26J/F81aS6uJN5ruJiZtmE7+9Mdy8uX17hXBUB00hccWiN3H21u5fEXr8INDOztqTB+YrZ1m8rKefLzJoRksAkd3+6ilWS9xtLRsVHFH+Ep5uPCUUF5ZVMXSutM5QtK4v/FXvdlS0riz8mcZXF8cSVQ6gg27/S8jZA89jrtsBSElRxFmdc7Sq8Pht4y3+snPokFl+b2OtdUhVXbL0uhMo7S8X5iu2zI9VXfp7KlhV5byf7XMUZ196EOq9elZa3AlpXeP0GcFIK49qj/LsjXFD/Fzt3cX3/yYor9vlOhHqEVqk6X7H/9keBO7axTtJ+Ywk7uenwR6hV/5BwUb0utuwmwl02QAvgqdj/GG8D+1bY9rrYdkuAk1Mc10zgK6Aw9vdcbHkvYGHsf4aFwIUpjutmYFHs+LOALhW2HRw7jx8BF6Qyrtj70cAtlbZL2vki3B2uADYSymAvBC4BLol9bsA9sZgXArkpOlc1xfUP4NsKv62C2PJ9Y+fpndh3fF2K47q8wm/rLSokqqq+/1TFFVtnEKHxSMXtkn2+jibUQbxb4bs6JVW/MQ0xISKS4RpTHYGIiNSBEoGISIZTIhARyXBKBCIiGU6JQEQkwykRiIhkOCUCEZEMp0QgkiBm1sHM+kcdh0htKRGIJM5PCWPdizQo6lkskgBmdjTwLPAd8D1wjoeB70TSnhKBSIKY2b+BEe5e5TSIIulKRUMiiXMA8EHUQYjUlhKBSALExtJf5WHmO5EGRYlAJDE6Al9EHYRIXSgRiCTGB0BbM3vPzHpFHYxIbaiyWEQkw+mJQEQkwykRiIhkOCUCEZEMp0QgIpLhlAhERDKcEoGISIZTIhARyXD/H1KIgHu5Z8a7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Vamos colocar a solução exata em um gráfico (lembre que lá em cima você deu o import)\n",
    "plt.figure()\n",
    "plt.title('')\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel(f'$y(t)$')\n",
    "plt.plot(t, yexato, 'b', label='Sol. Exata')\n",
    "plt.plot(t, y_euler, 'ro', markersize=5, label='Método de Euler')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "p0aG-zYHXIgk"
   },
   "source": [
    "Por que a diferença entre a solução exata e a solução pelo método de Euler aumenta com o tempo ? (resposta está nas suas aulas de EDO !) e como não queremos que esse *erro* que aumenta com o tempo, temos que usar outra solução para nosso problema."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q81hrkOT1AJD"
   },
   "source": [
    "### Usando um `solver` (ODEINT)\n",
    "\n",
    "Para calcular a mesma solução aproximada do seguinte Problema de Valor Inicial (P.V.I.) podemos usar agora um solver. Note que as condições iniciais são as mesmas:\n",
    "\\begin{cases}\n",
    "y'= f(y,t) = y + 2 t       \\\\\n",
    "y(0) = 1  \n",
    "\\end{cases}\n",
    "\n",
    "A diferença aqui é que usamos um solver do módulo `scipy.integrate`, conhecido como `spi` (veja em no [SciPy.org](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html) onde há vários métodos disponíveis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "Z2irsbhf1nvX"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate as spi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "id": "9QijTtzlmgMD"
   },
   "outputs": [],
   "source": [
    "#definição da EDO\n",
    "\n",
    "def f(y,t):\n",
    "    dy = y + 2*t\n",
    "    return dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "gbobBrt0hTLL"
   },
   "outputs": [],
   "source": [
    "# definição dos parâmetros do iniciais\n",
    "# *** note que estamos escrevendo aqui de novo alguns parâmetros que já definimos lá em cima\n",
    "# mas não seria necessário caso você já tivesse rodado a mesma linha com essas definições.\n",
    "h = 0.1 \n",
    "tmax = 2 # t máximo igual a 2 significa que o cálculo será feito entre t=0 e t=2, isto é, tmax = último valor de t\n",
    "# essa função gera vários valores de t entre 0 e tmax, separados por um passo=h, que serão colocados no array t\n",
    "t = np.arange(0, tmax+h, h) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "2Sf8jFVdhTmY"
   },
   "outputs": [],
   "source": [
    "#condição inicial\n",
    "y0 = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "DWg7u5izhlrv"
   },
   "outputs": [],
   "source": [
    "# aqui é a como usamos a Solução via odeint\n",
    "\n",
    "y_solver = spi.odeint(f, y0, t) #integra a equação diferencial e salva no array y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "OcMkuIodiEl6"
   },
   "outputs": [],
   "source": [
    "# mais uma vez a solução exata\n",
    "yexato = 3 * np.exp(t) - 2 * (t + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "id": "H78nw8gMiJxi",
    "outputId": "b1889095-a629-4f6e-ba3b-c68436ffa08f",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArOklEQVR4nO3deXhU5fn/8ffNGAgEiiAUWVSgiIgsCQbjT6hgtUBqFBGtYFFQEZRvtbhgtaVWa91ta0VlsSoi1mJB0aJToBQrIkQDCaAgioIaBIkoahIg2/P74wyYhCUJmZkzyXxe18WVmbPenhnnPudZzTmHiIjErwZ+ByAiIv5SIhARiXNKBCIicU6JQEQkzikRiIjEuaP8DuBItGrVynXs2NHvMERE6pRVq1Z96ZxrXXl5nUwEHTt2JCsry+8wRETqFDP75GDLVTQkIhLnlAhEROKcEoGISJyrk3UEB1NcXExubi579uzxOxQBEhMT6dChAwkJCX6HIiJVqDeJIDc3l2bNmtGxY0fMzO9w4ppzjp07d5Kbm0unTp38DkdEqlBvEsGePXuUBGKEmXHMMceQl5fndygi9UZpWSnBTUGyt2WT0jaF9C7pBBoEwnLsepMIACWBGKLPQiR8SstKGTx7MJlbMykoKiCpYRJp7dNYOGphWJJB1CqLzewpM9thZu9WWn6dmb1vZu+Z2QPRikdEpK4IbgqyMjeT/KJ8HI78onwyt2YS3BQMy/Gj2WpoJjCk/AIzOwsYCvR2zp0CPBTFeMLu7rvv5pRTTqFXr14kJyeTmZl52O0HDhwYtY5x8+fPZ/369VE5l4iEV/a2bAqKCiosKygqIGd7TliOH7WiIefcG2bWsdLia4H7nHN7Q9vsiFY84bZixQoWLFjA6tWradSoEV9++SVFRUV+h7Xf/PnzycjIoHv37n6HIiI19MOyFChOgob5+5clNUwi+djksBzf734EXYEfm1mmmf3PzPoeakMzG2dmWWaWFYuVkNu2baNVq1Y0atQIgFatWtGuXTsAlixZQkpKCj179uTKK69k79691T5uXl4ew4cPp2/fvvTt25fly5cDMHToUGbNmgXA9OnT+cUvfgHAE088Qd++fenduzfDhw+nsLCQt956i1deeYVJkyaRnJzMRx99dNDtRCT2OAf/vCedwPY0mhzVFMNo2rApae3TSO+SHpZzWDSnqgw9ESxwzvUIvX8XWApcD/QF5gCdXRVBpaamuspFKhs2bODkk08GYOJEyMkJb+zJyfDww4den5+fT//+/SksLOScc87hkksuYcCAAezZs4cTTzyRJUuW0LVrVy6//HL69OnDxIkTGThwIA899BCpqamHPO6ll17KhAkT6N+/P59++imDBw9mw4YNfPHFF/Tr14+nn36aq666ipUrV9KyZUt27tzJMcccA8DkyZNp06YN1113HWPGjCEjI4OLLroI4JDbhVP5z0REjsycOTBiBDwypZROg4LkbM8h+djkI2o1ZGarnHMH/OD43WooF3gx9MP/tpmVAa2A2Lvlr0LTpk1ZtWoVy5YtY+nSpVxyySXcd999pKSk0KlTJ7p27QrA6NGjeeyxx5g4cWK1jvuf//ynQtn+t99+S35+Pm3atOEPf/gDZ511Fi+99BItW7YE4N1332Xy5Mns2rWL/Px8Bg8efNDjVnc7EfHPt9/CDTfAqafChGsDBAIZZHTNCPt5/E4E84GzgKVm1hVoCHxZ24Me7s49kgKBAAMHDmTgwIH07NmTZ555hpSUlFods6ysjJUrV5KYmHjAunXr1nHMMcfw+eef7182ZswY5s+fT+/evZk5cyavv/76QY9b3e1ExD+33w7bt8PLL0MgPF0GDiqazUefB1YAJ5lZrpldBTwFdA4VEf0DGF1VsVCs2rhxIx9++OH+9zk5OZxwwgmcdNJJbNmyhU2bNgHw7LPPMmDAgGofd9CgQUyZMqXCcQHefvttgsEg2dnZPPTQQ2zevBmA7777jrZt21JcXMxzzz23f79mzZrx3Xff7X9/qO1EJDZkZ8OUKXDNNdD3kLWn4RG1ROCcG+mca+ucS3DOdXDOPemcK3LOjXLO9XDO9XHO/Tda8YRbfn4+o0ePpnv37vTq1Yv169dzxx13kJiYyNNPP83FF19Mz549adCgAddcc80B+48dO/agTUkfeeQRsrKy6NWrF927d2fatGns3buXq6++mqeeeop27drxpz/9iSuvvBLnHHfddRdpaWn069ePbt267T/OiBEjePDBB0lJSeGjjz465HYi4r+yMrj2WmjVCu6+O/Lni2plcbhUVVkssUGficiRmTEDxo+HWbPgssvCd9xDVRb73XxURETK2bEDbr0VBgyAUaOic04lAhGRGPLrX8N338Hjj0O0huxSIhARiRHLlsHMmXDzzRDNQQCUCEREYkBxsVdBfMIJMHlydM/tdz8CERHB6//03nten4GkpOieW08EIiI+++wzuOMOOO88OP/86J9ficAn0RyCurItW7bQo0ePw27z+eef7x+X6HDuueeecIUlErd+9StvcLlHHvHn/HGbCErLSlnwwQLu+t9dLPhgAaVlpX6HdERKSyMTd7t27Zg7d26V2ykRiNTOq6/CSy/B734HHTv6E0NcJoJ9076NnDeS37/+e0bOG8ng2YNrlQwKCgo499xz6d27Nz169GDOnDlA1UNQT5s2jUmTJu1/P3PmTH75y18CMHv2bE477TSSk5MZP378/h/9pk2bctNNN9G7d29WrFhR4Xg5OTmcfvrp9OrVi2HDhvH1118DsGrVKnr37k3v3r157LHHvr8WpaVMmjSJvn370qtXL6ZPnw5UfGqYOXMmF154IUOGDOHEE0/klltuAeDWW29l9+7dJCcn7x8GW0Sqr7AQrrsOunWDm27yL464TATBTUEyt4Z32rd///vftGvXjjVr1vDuu+8yZMgQ9uzZw5gxY5gzZw7r1q2jpKSEqVOnVthv+PDhvPTSS/vfz5kzhxEjRrBhwwbmzJnD8uXLycnJIRAI7B8TqKCggLS0NNasWUP//v0rHO/yyy/n/vvvZ+3atfTs2ZM777wTgCuuuIIpU6awZs2aCts/+eSTNG/enHfeeYd33nmHJ554Yv+4ReXl5OTs/++YM2cOn332Gffddx+NGzcmJydH4xWJHIF77oHNm70+Aw0b+hdHXCaCSEz71rNnTxYvXsyvf/1rli1bRvPmzdm4ceMBQ1C/8cYbFfZr3bo1nTt3ZuXKlezcuZP333+ffv36sWTJElatWkXfvn1JTk5myZIlfPzxx4A3yunw4cMPiOGbb75h165d+we123e+Xbt2sWvXLs4880wALivXZ33RokXMmjWL5ORk0tLS2LlzZ4XB8/Y5++yzad68OYmJiXTv3p1PPvnkiK+ViMDGjfDAA17v4bPO8jeWuGw+mtI2haSGSeQXhW/at65du7J69Wpee+01Jk+ezNlnn83QoUOrte+IESN44YUX6NatG8OGDcPMcM4xevRo7r333gO2T0xMJBCmMWmdc0yZMuWA+Qi2bNlS4f2+mdfAS0QlJSVhOb9IPHIOJkyAJk3goRiYqT0unwjSu6ST1j6Npg3DN+3b559/TpMmTRg1ahSTJk1i9erV1R6CetiwYbz88ss8//zzjBgxAvDuwOfOncuOHd40zl999VWVd+HNmzenRYsWLFu2rML5jj76aI4++mjefPNNgArFOIMHD2bq1KkUFxcD8MEHH1BQUHDgwQ8hISFh/74iUj3/+Af8979e0VCbNn5HE6dPBIEGARaOWkhwU+2mfStv3bp1TJo0iQYNGpCQkMDUqVMrDEFdUlJC3759DzoEdYsWLTj55JNZv349p512GgDdu3fnj3/8I4MGDaKsrIyEhAQee+wxTjjhhMPG8cwzz3DNNddQWFhI586defrppwF4+umnufLKKzEzBg0atH/7sWPHsmXLFvr06YNzjtatWzN//vxq/3ePGzeOXr160adPH9UTiBxGaVkpwU1BVmzJ5vGHUzg1NZ3x4yM420wNaBhqiRh9JiKefS0VM7dmkr+3AIqS6NsujRUTFtbqBrSmNAy1iIhPyrdUxBw0ymfDd7VrqRhO0Zyq8ikz2xGalrLyupvMzJlZq2jFIyISLZFoqRhO0XwimAkMqbzQzI4DBgGf1vYEdbGYq77SZyHyvZS2KRzlKo4kV9uWiuEUzTmL3wC+OsiqvwC3ALX65UhMTGTnzp36AYoBzjl27txJYmKi36GIxIQf7EineHMaR5WFr6ViOPnaasjMhgJbnXNrrIqpeMxsHDAO4Pjjjz9gfYcOHcjNzSUvLy8SoUoNJSYm0qFDB7/DEPFdfj5cOSbACWULeeDGIB98E56WiuHkWyIwsybAb/CKharknJsBzACv1VDl9QkJCXTq1CmsMYqI1NYtt8DHH8PSpQEG9M4AMvwO6QB+thr6EdAJWGNmW4AOwGozO9bHmEREwmbRIpg6FSZO9Cajj1W+PRE459YBP9z3PpQMUp1zX/oVk4hIuOzaBVdeCSefDHff7Xc0hxfN5qPPAyuAk8ws18yuita5RUSi7frrYft2mDULGjf2O5rDi9oTgXNuZBXrO0YpFBGRiHrpJXj2Wbj9dkg9oB9v7FHPYhGRMNqxA8aPhz59YPJkv6OpHiUCEZEwcc5LAt9+6xUJJST4HVH1xOXooyIikfDsszB/Pjz4IJxyit/RVJ+eCEREwuCzz7z5h3/8Y7jhBr+jqRklAhGRWior85qKlpbCzJkQpgkEo0ZFQyIitTR1KvznPzBtGnTu7Hc0NacnAhGRWvjwQ5g0CYYMgXHj/I7myCgRiIgcodJSGD0aGjWCv/0Nqhg7M2apaEhE5Ag9+CCsWAHPPQft2/sdzZHTE4GIyBFYu9brOXzRRTDysOMmxD49EYiIVFNpWSnBTUGycrOZ9UAKLY5JZ+rUQJ0tEtpHiUBEpBpKy0oZPHuwNwn93gLok0SvQWm0aLkQqGPtRStR0ZCISDUENwW9JFCUD+agUT4f780kuCnod2i1pkQgIlIN2duyKSgqqLCsoKiAnO05/gQURkoEIiLVkHxsCg1KkyosS2qYRPKxyf4EFEZKBCIi1bD+5XRKP0mjIU0xjKYNm5LWPo30Lul+h1ZrqiwWEanCsmXw298EGD5sIaNHBFnzRQ7JxyaT3iWdQIO6XVEMYM656JzI7CkgA9jhnOsRWvYgcB5QBHwEXOGc21XVsVJTU11WVlYEoxUR8XzxBaSkQNOmkJUFP/iB3xEdOTNb5Zw7YM60aBYNzQSGVFq2GOjhnOsFfADcFsV4REQOq7TU6yy2axfMm1e3k8DhRC0ROOfeAL6qtGyRc64k9HYl0CFa8YiIVOX3v4elS+Hxx6FnT7+jiZxYqiy+Eqj7DXJFpF4IBuHuu+Gqq2DMGL+jiayYSARm9lugBHjuMNuMM7MsM8vKy8uLXnAiEnc++QRGjYLevWHKFL+jiTzfE4GZjcGrRP6FO0zNtXNuhnMu1TmX2rp166jFJyLxZe9e+PnPoaQE5s6Fxo39jijyfG0+amZDgFuAAc65Qj9jEREBuPlmePttr3K4Sxe/o4mOqD0RmNnzwArgJDPLNbOrgEeBZsBiM8sxs2nRikdEpLI5c+DRR+HGG+HCC/2OJnqi9kTgnDvYiN1PRuv8IiKH8/77MHYsnHEG3Hef39FEl+91BCIifiso8CaYadwYXngBEhL8jii6NMSEiMQ15+Caa2D9eli0qG5POXmklAhEJK498QTMng133gnnnON3NP5Q0ZCIxK3Vq+H662HwYJg82e9o/KMnAhGJK/vmHX5rczZP3Z1Cqx+mM3t2gAZxfFusRCAicWP/vMO5meQXFUD/JE69uH7MO1wbcZwDRSTe7J93uPj7eYc35tePeYdrQ4lAROJGfZ53uDaUCEQkbvygMAVXVD/nHa4NJQIRiQtbt8L9V6fTKC+NJkfVv3mHa0OVxSJS7+XnQ0YG5H8bYMX/LWRr4yA52+vXvMO1oUQgIvXavukm162DBQsgpXeAFDLI6Jrhd2gxQ4lAROq1G2/0EsDjj8OQyrOmC6A6AhGpx6ZMgUcegRtugGuv9Tua2KVEICL10quvwsSJMHQoPPig39HENiUCEal3cnLgkksgJQWeew4C8V0XXCUlAhGpV7Zu9VoItWgBr7wCSUlV7xPvojlV5VNmtsPM3i23rKWZLTazD0N/W0QrHhGpf/Lz4bzz4JtvvKKhdu38jqhuiOYTwUygcp39rcAS59yJwJLQexGRGisthUsvhTVrvFnGevXyO6K6I2qJwDn3BvBVpcVDgWdCr58BLohWPCJSv9x0E/zrX15LofT47ihcY37XEbRxzm0Lvd4OtPEzGBGpmx59FP76V6+Z6IQJfkdT9/idCPZzzjnAHWq9mY0zsywzy8rLy4tiZCISy159FX71Kzj/fDUTPVJ+J4IvzKwtQOjvjkNt6Jyb4ZxLdc6ltm7dOmoBikjsKS0rZcEHC5gw5y4u+s0CeieX8ve/q5nokfJ7iIlXgNHAfaG/L/sbjojEun2zjK34LJPC4gLsvCSSTkgjsXF8zzJWG9FsPvo8sAI4ycxyzewqvATwUzP7EDgn9F5E5JCCm4KszM2ksMSbZcwl5JPzpWYZq42oPRE450YeYtXZ0YpBROq+zE9Ds4zZ98v2zTKmEUWPjN91BCIi1bZ7N7wyPQU0y1hYKRGISJ1QVAQXXQRrX0yne/M0mjbULGPh4ndlsYhIlUpKvMllXnsNpk8PcNXYhQQ3aZaxcFEiEJGYVloKo0fDiy/Cww/DuHEAATK6apaxcFHRkIjErLIyGD8e/v53uPder+OYhJ8SgYjEJOe8H/4nn4Tf/Q5u1ZCUEVPjRGBmSWamwjgRiRjn4Ne/9sYQuvlmuPNOvyOq36pMBGbWwMwuNbNXzWwH8D6wzczWm9mDZtYl8mGKSDy5805v3KAJE+CBB8Cs6n3kyFXniWAp8CO8uQKOdc4d55z7IdAfWAncb2ajIhijiMSR++/3EsEVV3hDSisJRF51Wg2d45wrNrPVzrk++xY6574ys63OueFmlhDBGEUkTjzyiFcXMHIkPPEENFAtZlRU5zIPM7P7gGZmdrKZld9nBoBzrjgi0YlI3HjiCa9yeNgweOYZjSQaTdV5IlgOJAJjgT/jDRq3C/gc2B250EQkXsye7TUTTU+H55+HBJUxRFWVicA5txWYZWYfOeeWA5jZMUBHvIpjEZEaKS0rJbgpSPa2bAo/TuH+q9MZODDAvHnQqJHf0cWfKhOBmZnzLN+3zDm3E9hZeZsIxSgi9ci++QQyt2ZSUFSA25tE8wlpvPTHhTRurPIgP1Sr1ZCZXWdmx5dfaGYNzewnZvYM3qQyIiJVCm4Kkrk1k/yifBwOGuVTcmwmy7ZrPgG/VCcRDAFKgefN7PNQ/4GPgQ+BkcDDzrmZEYxRROqR7G2h+QTKKSz25hMQf1SnjmAP8DjweKiZaCtgt3NuV4RjE5F6KDcrBbc3CRrl71+m+QT8Ve1Wuma2GOjunNumJCAiNeUc3HEHzJiUTqu9aSQlaD6BWFGTYah/DTxsZluA3zjntoUrCDO7Aa95qgPWAVeEnkREpB4oK/P6CDz6KIwZE2DavQtZvEXzCcQKq2ljHzMbDtwOvAg84JyrVV8CM2sPvIn3tLHbzF4AXjtcvUNqaqrLysqqzWlFJEqKi2HMGG8o6Rtv9MYQUo9hf5jZKudcauXlNfo4zMyAjcBU4DrgQzO7LAzxHQU0NrOjgCZ4ndVEpI4rLIQLLvCSwD33wEMPKQnEoprUESwHtgJ/AdoDY4CBwGlmNuNIAwh1WHsI+BTYBnzjnFt0kPOPM7MsM8vKy8s70tOJSJTs2gWDB0MwCNOnw223aQC5WFXtoiEzOwVYf7COY2a2wTl38hEFYNYCmAdcAuwC/gnMdc7NPtQ+KhoSiW3bt3tJYMMGeO45uPhivyMSCEPRkHPuvcP0Hj73iCODc4DNzrm80OB1LwJn1OJ4IuKjzZuhf3/YtAkWLFASqAvCMnm9c+7jWuz+KXC6mTXBG8TubEC3+yJ10Lp13pPAnj2wZAmcfrrfEUl1+F5t45zLBOYCq/GajjYgNLy1iNQdK1bAmWd69QDLlikJ1CVheSKoLefc74Hf+x2HiFRf+RFE3bYU7hubTvu2ARYvho4d/Y5OaiImEoGI1C0HG0G06eg0/jd5Ie3aqmNYXeN70ZCI1D3BTUEycyuOIOraZ7L6O40gWhcpEYhIjb2Tm02+RhCtN5QIRKRGduyAfz6SAkVJFZZrBNG6S4lARKotJwdSU2Hz4nR6tEijaUONIFofqLJYRKrln//0Bo9r2RKWLwvQO3khwU0aQbQ+UCIQkcMqK/PmEbjrLjjjDHjxRWjTBiBARtcMMrpm+Byh1JYSgYgcUn4+XHYZzJ8PV1wBU6dCo0Z+RyXhpkQgIge1eTMMHQrvvQcPPwzXX6/RQ+srJQIROcDrr8NFF0FpqTeM9KBBfkckkaRWQyJSwbRp8NOfQuvW8PbbSgLxQE8EInGs/HhBPVun8O9H05k+NcDPfubNKta8ud8RSjQoEYjEqcrjBVlJEmUFadx8y0LuuydAQC1B44aKhkTiVHBTkMyt348XVHZUPo1+lMmAq4JKAnFGiUAkTq3eduB4QUVO4wXFIyUCkTj01Vfw76dTYK/GCxIlApG488Yb0Ls3vP33dLo0TqNpgsYLincxUVlsZkcDfwN6AA640jm3wtegROqZkhL44x+9oSI6dYIVywP0OVXjBUmMJALgr8C/nXMXmVlDoInfAYnUJ59+Cr/4Bbz5pjdkxGOPQbNmoPGCBGIgEZhZc+BMYAyAc64IKPIzJpH6ZN48GDvWeyJ49lkYNcrviCTWxEIdQScgD3jazLLN7G9mllR5IzMbZ2ZZZpaVl5cX/ShF6pjCQhg/3hsq4sQTITtbSUAOLhYSwVFAH2Cqcy4FKABurbyRc26Gcy7VOZfaunXraMcoUqesXQt9+8KMGXDLLV6RUJcufkclscr3oiEgF8h1zmWG3s/lIIlARA6u/DARycemsGVxOpNuDnD00bBokTdukMjh+J4InHPbzewzMzvJObcROBtY73dcInVB5WEiGpQmUfpJGkN+spBnZgb44Q/9jlDqglgoGgK4DnjOzNYCycA9/oYjUjdUHiaiNJBPw86ZXPvnoJKAVFtMJALnXE6o/L+Xc+4C59zXfsckUhes2JJN/t6Kw0QUU8DaHTn+BCR1UkwkAhGpuZdfhml3pECRhomQ2lEiEKljtm+Hiy+GCy6A9gXp9G2XRtOGGiZCjpzvlcUiUj3OwVNPwc03w+7dcPfdMGlSgAYBDRMhtaNEIFIHbNoE48bB0qVw5ple/4CTTtq3VsNESO2oaEgkhhUXw/33Q8+esHq1lwCWLi2fBERqT08EIjFq1SpvjKCcHLjwQpgyBdq18zsqqY+UCERiQPnewd1bpvDWrHQe/nOANm28QeMuvNDvCKU+UyIQ8Vnl3sEUJ+F2pDF27EIevN8bKkIkklRHIOKz4KYgKz/7vnewS8incZdMht4UVBKQqFAiEPHRnj3w1xeyvSeB8svLNIm8RI8SgYgPnIN//AO6dYP/PJtCwKl3sPhHiUAkylasgDPOgJEjoUULWDw1nYFd1DtY/KPKYpEo2bIFbr0V5syBtm29XsKXXw6BQICzytQ7WPyjRCASYd9+C/feC3/5CzRoALffDpMmQdOm328TaKDeweIfJQKRCCkpgSefhN/9DvLy4LLL4J57oEMHvyMTqUiJQCQMyncIS2mbQmBzOrfcHODdd+HHP4bXXoPUVL+jFDk4JQKRWjrUdJGddi9k3rwAw4aBmd9RihxazLQaMrOAmWWb2QK/YxGpieCmICs+O3C6yD+9HOTCC5UEJPbFTCIAfgVs8DsIkZrYuBEmP5ZNYdGB00W+tzPHn6BEaigmEoGZdQDOBf7mdywi1bF+PVx6KZx8Mmz4bwoJpg5hUnfFRCIAHgZuAcoOtYGZjTOzLDPLysvLi1pgIuWtWweXXAI9esArr3jNQLcsTufMzuoQJnWX75XFZpYB7HDOrTKzgYfazjk3A5gBkJqa6qITnYgnJwfuugtefBGaNYPbboMbboBWrQACLBylDmFSd/meCIB+wPlm9jMgEfiBmc12zo3yOS6JM5WbgKZ3SScnO8Bdd8HLL8MPfuD1CZg4EVq2rLivOoRJXeZ7InDO3QbcBhB6IrhZSUCirXIT0MRAEk2+TmPnXxZydPMAd94J11+PhoWWesn3RCASC4KbgmTmZpJfnA/A7tJ8didmMurOII9en0Hz5j4HKBJBsVJZDIBz7nXnnJ6tJap274YnX80mv1ITUGtYwEkDcpQEpN6LqUQgEk25ufCb38Bxx8H8aSlYiZqASnxSIpC44hwsX+41Ae3YEe6/3xsLaPHUdM7qqiagEp9URyBxYc8ebx6ARx6B1au9St8bboAJE6BTJwDNCSDxS4lA6pXKTUB7N0lnxvQA06d7Q0F37w7TpsGoUZBUsSRITUAlbikRSL2xvwlobiYFxQUEypIo+SQNnl3I+RkBrr8efvITDQInUpkSgdQbc1YHWbY5kyK8JqAlDfJJ6JjJ40uDjD1Td/kih6LKYqnTiou9Xr8XXgiX3ZJNUVnFJqAlVsB2cvwJTqSOUCKQOiknxxvqoX17uOACeOstOL9vCk0S1ARUpKaUCKTO+OILbwL43r0hJQWmToUBA2DBAq9PwNx70/l/x6sJqEhNqY5AYk75lj89WqVQvCGdZ58JEAxCaSmcdho89hiMGFF58DeNAipyJJQIJKaUlpUyaPZgVnyaye6SAihOgtw02mYv5OabA4we7U0GcyhqAipSc0oEEhNKS71y/gdfCbK0USYuIR8MaJhPYpdMpt4WZGg3/biLRILqCMQ3RUWwcCGMHw/t2sGZZ8Jrq7NxR1Vs+bO3rIB1O3L8CVIkDuiJQMLuYBO87CunLyyERYtg3jz417/gm2+8Hr7nngvDhwMnpXDVa0nkF+XvP55a/ohElhKBhFXlCV6SGiZxaps0xjZayPyXvArfwkJo0QKGDfPa///0p5CYuG//dGasSauwv1r+iESWEoGEVXBTkMytmfvv6POL8vnfR5n8b26QY7/NYPRo78d/wABISDhw/0ADtfwRiTYlAgmLwkJYtgzueT2b/EYFXkXvPgkFXP27HKZdmkGDatRKqeWPSHT5ngjM7DhgFtAGcMAM59xf/Y1KquIcrFvnVfYuWuQlgb174ahTUghcmERp4Psy/qaNkjj/tORqJQERiT7fEwFQAtzknFttZs2AVWa22Dm33u/A4tnBKny/zAuweLH3w794MWzf7m17yineuP6DBkG//ukMe1Fl/CJ1ie+JwDm3DdgWev2dmW0A2gNKBD7ZV+G7MjeTwuICjnJJNMpLI3/qQnABWrXyKngHDfL+tm9ffm+V8YvUNeac8zuG/cysI/AG0MM5922ldeOAcQDHH3/8qZ988kn0A6znduzwpnGctXIBrzQcSdlR3xfvNChpyqjE57l+SAYpKaiYR6QOMrNVzrnUyst9fyLYx8yaAvOAiZWTAIBzbgYwAyA1NTV2slcMO1x7fudg40Z4803vx3/5cvjwQ2+/wE+yKftxxU5d7qgCTvxxDqeeqgpckfomJhKBmSXgJYHnnHMv+h1PfVC5PX+ThCS6NU3jot0LeWt5gLfegp07vW2POQb694err4Z+/eCL5ilc/oo6dYnEC98TgZkZ8CSwwTn3Z7/jqQ9KS2H60iBvbslkr/N+zAuK81m1PZNV84J0JYPzz/d+/Pv1g65dK07fWFqWTtoqVfiKxAvfEwHQD7gMWGdmOaFlv3HOveZfSLHjcMU74DXZfO89WL0asrO9v2vXQmFqNgwsqDCalDUq4JaHc7gv/fDFO+rUJRJffE8Ezrk3qdj9SEIOVbxzuS0kJztAdraXBIqLve2bNfMmbLn6amjQLYXpXyZRWFqxeKf/j5KrdW516hKJH74ngnhQ1V19ZXv2wAcfwOx3gizL/X4y9vLFO62+yqBPHxgyxPvx79MHOnf+vjVPaVk6a2ereEdEqqZEEGEHG4QtrX0aC0ct5NtvAmzYABs2wPvvs//15s1eqx7OPHjxzk0P5fBARkaFcv3KVLwjItWlRFANNb2jL++VDUFWfJZJYcn3g7At/TCTY04P8s073xe7NGrkVdqmpsKoUd4sXHktUrhtVRIFxRWLdwaclHzYJLCPindEpDqUCKpwuDv6QIMAJSXw2WewZYt3J7/v777XW7tkw4CKd/VlgQJ+1C+HkT/PoFs370e/Y0cIBCqfO52Xv1DxjohEVtwkgiO5q3cO/plz4B39/z7KJPnnQb7LyiA312uuuU+DBtChA3Tq5A2/sLdjCi9a0v5mnOANwnbntclkdD18zCreEZFoiItEcMBdfUISKW3SmHLaQnZsD/D557Btm/ev8us9px1YTl9CAd8k5nBm/ww6dvR+9Dt18u7qjzuu4jj7pWXp7KhFpa2Kd0Qk0uIiERwwWUpxPss+yiT5gSB88P0PbLNm0LatN39uWpr39+vWKTy3t9IdfWISj99e9R096K5eRGJfXCSC7G3ZFBRVHDuHRgX8/Poc/q9nBm3begmgadMD9y0tS+eTWjbD1F29iMSyuEgEKW1TSGpYceycpg2TuOynyZypcnoRiXNxkQjSu6ST1l7l9CIiBxMXiUB39SIihxYXiQB0Vy8iciiaZ0pEJM4pEYiIxDklAhGROKdEICIS55QIRETinDnn/I6hxswsD/jkCHdvBXwZxnDCRXHVjOKqGcVVM7EaF9QuthOcc60rL6yTiaA2zCzLOZfqdxyVKa6aUVw1o7hqJlbjgsjEpqIhEZE4p0QgIhLn4jERzPA7gENQXDWjuGpGcdVMrMYFEYgt7uoIRESkonh8IhARkXKUCERE4ly9SgRmNsTMNprZJjO79SDrG5nZnND6TDPrWG7dbaHlG81scJTjutHM1pvZWjNbYmYnlFtXamY5oX+vRDmuMWaWV+78Y8utG21mH4b+jY5yXH8pF9MHZrar3LqIXC8ze8rMdpjZu4dYb2b2SCjmtWbWp9y6SF6rquL6RSiedWb2lpn1LrduS2h5jpllRTmugWb2TbnP6vZy6w77+Uc4rknlYno39H1qGVoXyet1nJktDf0OvGdmvzrINpH7jjnn6sU/IAB8BHQGGgJrgO6VtpkATAu9HgHMCb3uHtq+EdApdJxAFOM6C2gSen3tvrhC7/N9vF5jgEcPsm9L4OPQ3xah1y2iFVel7a8DnorC9ToT6AO8e4j1PwOCgAGnA5mRvlbVjOuMfecD0vfFFXq/BWjl0/UaCCyo7ecf7rgqbXse8N8oXa+2QJ/Q62bABwf5/zFi37H69ERwGrDJOfexc64I+AcwtNI2Q4FnQq/nAmebmYWW/8M5t9c5txnYFDpeVOJyzi11zhWG3q4EOoTp3LWK6zAGA4udc185574GFgNDfIprJPB8mM59SM65N4CvDrPJUGCW86wEjjaztkT2WlUZl3PurdB5IXrfrepcr0Opzfcy3HFF5bsF4Jzb5pxbHXr9HbABaF9ps4h9x+pTImgPfFbufS4HXsj92zjnSoBvgGOquW8k4yrvKrysv0+imWWZ2UozuyBMMdUkruGhx9C5ZnZcDfeNZFyEitA6Af8ttzhS16sqh4o7kteqpip/txywyMxWmdk4H+L5f2a2xsyCZnZKaFlMXC8za4L3Yzqv3OKoXC/ziqxTgMxKqyL2HYubGcrqAjMbBaQCA8otPsE5t9XMOgP/NbN1zrmPohTSv4DnnXN7zWw83tPUT6J07uoYAcx1zpWWW+bn9YpZZnYWXiLoX25x/9C1+iGw2MzeD90xR8NqvM8q38x+BswHTozSuavjPGC5c67800PEr5eZNcVLPhOdc9+G89iHU5+eCLYCx5V73yG07KDbmNlRQHNgZzX3jWRcmNk5wG+B851ze/ctd85tDf39GHgd704hKnE553aWi+VvwKnV3TeScZUzgkqP7hG8XlU5VNyRvFbVYma98D6/oc65nfuWl7tWO4CXCF9xaJWcc9865/JDr18DEsysFTFwvUIO992KyPUyswS8JPCcc+7Fg2wSue9YJCo+/PiH93TzMV5Rwb5KplMqbfN/VKwsfiH0+hQqVhZ/TPgqi6sTVwpeBdmJlZa3ABqFXrcCPiRMFWfVjKttudfDgJXu+8qpzaH4WoRet4xWXKHtuuFV3lk0rlfomB05dOXnuVSsyHs70teqmnEdj1fndUal5UlAs3Kv3wKGRDGuY/d9dng/qJ+Grl21Pv9IxRVa3xyvHiEpWtcr9N8+C3j4MNtE7DsWtosbC//watU/wPtR/W1o2R/w7rIBEoF/hv7HeBvoXG7f34b22wikRzmu/wBfADmhf6+Elp8BrAv9z7AOuCrKcd0LvBc6/1KgW7l9rwxdx03AFdGMK/T+DuC+SvtF7Hrh3R1uA4rxymCvAq4BrgmtN+CxUMzrgNQoXauq4vob8HW571ZWaHnn0HVaE/qMfxvluH5Z7ru1knKJ6mCff7TiCm0zBq/xSPn9In29+uPVQawt91n9LFrfMQ0xISIS5+pTHYGIiBwBJQIRkTinRCAiEueUCERE4pwSgYhInFMiEBGJc0oEIiJxTolAJEzMrIOZXeJ3HCI1pUQgEj5n4411L1KnqGexSBiYWX/gZWAX8B1wofMGvhOJeUoEImFiZv8GbnbOHXQaRJFYpaIhkfA5CXjf7yBEakqJQCQMQmPpf+O8me9E6hQlApHw6Ah87ncQIkdCiUAkPN4HWpnZu2Z2ht/BiNSEKotFROKcnghEROKcEoGISJxTIhARiXNKBCIicU6JQEQkzikRiIjEOSUCEZE49/8BO0CCGV4R9BMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# gráfico aqui (você pode explorar ainda elementos do gráfico para mudar a configuração dos pontos, cores das linhas, etc.)\n",
    "plt.figure()\n",
    "plt.title('')\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel(f'$y(t)$')\n",
    "plt.plot(t, yexato, 'b', label='Sol. exata')\n",
    "plt.plot(t, y_solver, 'go', markersize=5, label='solver odeint')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notou que não tem mais diferença entre a solução exata e a solução via ODEINT ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 390
    },
    "id": "PQJd2t2fagmU",
    "outputId": "85dd7f8f-9ec6-47bf-a972-b3d14209666b"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1966817/3854516553.py:19: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
      "  fig.show()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAF1CAYAAADRBwbsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABedElEQVR4nO3debzUY//H8dd1pk7LnBBFKSpSiMwkjl22TsO57cQtssRtX37itkf2fe1GloosuZUtxpFukTCUM9ZKIZQiUZqpszTn+v1xneqc1rPOd5b38/Hocc585zvf+TTn6OPzva7rcxlrLSIiIiIiIiKpLsfrAERERERERERqQgWsiIiIiIiIpAUVsCIiIiIiIpIWVMCKiIiIiIhIWlABKyIiIiIiImlBBayIiIiIiIikhSZeB1AXbdq0sZ07d/Y6DBERyRDTpk37w1rb1us40plys4iINKT15ea0LGA7d+7M1KlTvQ5DREQyhDHmJ69jSHfKzSIi0pDWl5s1hVhERERERETSggpYERERERERSQsqYEVERERERCQtpOUa2HUpLy9n7ty5lJSUeB2KNIDmzZvTsWNHmjZt6nUoIiJSR8rNqUf5VUTSXcYUsHPnzqVVq1Z07twZY4zX4Ug9WGtZtGgRc+fOpUuXLl6HIyIidaTcnFqUX0UkE2TMFOKSkhK22GILJcgMYIxhiy220B17EZE0p9ycWpRfRSQTZEwBCyhBZhD9LEVEMoP+PU8t+nmISLrLqALWa7feeis9evSgZ8+eBAIBIpHIBs/v06fPRvfMu/HGG+nQoQOBQGDVn8WLF9c6tmg0yltvvVXr14mIiKQz5WYRkcySMWtgvfbxxx8zfvx4Pv/8c5o1a8Yff/xBWVlZg1z7sssuY/DgwfW6RjQaZerUqRx++OENEpOIiEiqU24WEck8GoFtIPPnz6dNmzY0a9YMgDZt2rD11lsDMHHiRILBILvuuitnnnkmpaWl9X6/+++/nzPPPBOAr776il122YVly5bx6aefsvfeexMMBtlnn32YOXMmZWVl3HDDDYwZM4ZAIMCYMWPWeZ6IiEgmUW4WEUmOREWC8d+N5+b3b2b8d+NJVCQa7b0ycgT20kshGm3YawYC8MAD63++b9++DB06lG7dunHooYfSv39/DjzwQEpKSjj99NOZOHEi3bp147TTTuPRRx/l0ksvrfF733///YwePRqA1q1b895773HJJZfQp08fXnnlFW699VYef/xxWrZsyY477sjkyZNp0qQJ7777Ltdccw1jx45l6NChTJ06lUceeQSAv//+e53niYiINAblZuVmEclMiYoEBaMLiMyLEC+L48/1k98hn6IBRfhyfA3+fhlZwHohLy+PadOmMXnyZN577z369+/PHXfcQTAYpEuXLnTr1g2AgQMHMmzYsFolyXVNU8rJyWHkyJH07NmTf/3rX+y7774ALFmyhIEDBzJr1iyMMZSXl6/zmjU9T0QkFf3+O0ybBv36gXrSyPooN4uINL7w7DCReRFiZTEAYmUxIvMihGeHKexW2ODvl7QC1hjzNFAI/G6t3aXK8YuAC4AE8Ka19sr6vteG7sY2Jp/PR58+fejTpw+77roro0aNIhgMNtr7zZo1i7y8PH799ddVx66//noOOuggXnnlFebMmUOfPn3W+dqanicikoruvRfuvhtmzYLtt/c6mvSl3NzwlJtFJNsUzy8mXhavdixeFie6INooBWwy18COBPpVPWCMOQg4CtjNWtsDuCeJ8TSomTNnMmvWrFWPo9EonTp1onv37syZM4fZs2cD8Oyzz3LggQfW+/2WLFnCxRdfzAcffMCiRYt4+eWXVx3v0KEDACNHjlx1fqtWrVi6dGm116/rPBGRVPfHHzBsGJx0korXBjAS5WblZhGRegi2D9Kyib/aMX+un0C7QKO8X9IKWGvtB8Cfaxw+D7jDWltaec7vyYqnocViMQYOHMjOO+9Mz549+fbbb7nxxhtp3rw5I0aM4IQTTmDXXXclJyeHc889d63XDxo0aL1t+++///5qrfrnzJnDZZddxgUXXEC3bt146qmnuOqqq/j999+58sorufrqqwkGg6xYsWLVNQ466CC+/fbbVY0i1neeiEiqu+8+WLYMrrvO60jSn3KzcrOISH2FuobYoiQfSvMwGPJy88jvkE+oa6hR3s9Yaxvlwut8M2M6A+NXTlMyxkSB13B3f0uAwdbazzZ2nd69e9s1E8r06dPZaaedGjpk8ZB+piKypkWLoHNnOPxwGDOm4a5rjJlmre3dcFdMH8rN2Uc/FxFpSAsXQqcuCfY8Jcyhp0QJtAsQ6hqqdwOn9eVmr5s4NQE2B/YC9gBeMsZsZ9dRVRtjzgHOAdh2222TGqSIiKSGBx6AWAyuv97rSDKacrOIiNTYffdByTIfj11WyI47Nvya1zV5vQ/sXGCcdT4FKoA26zrRWjvcWtvbWtu7bdu2SQ1SRES899df8NBDcNxxsMsuGz9f6ky5WUREamTRInjkEejfH3bcMTnv6XUB+ypwEIAxphuQC/zhZUAiIpKaHnwQ/v5bo69J8CrKzSIiUgP33w/xeHL7UiRzG50XgD5AG2PMXGAI8DTwtDHma6AMGLiuKUoiIpLdFi9204ePPhp2283jYDKIcrOIiNTVn3+6mVHHHw89eiTvfZNWwFprT17PUwOSFYOIiKSnhx+GJUvghhu8jiSzKDeLiEhdPfggLF2a/JlRXk8hzlqPP/44f/31l9dhiIikvL//dlOU/vEPCAa9jkZEREQWL3YF7LHHwq67Jve9VcA2IGMMAwasvmm9YsUK2rZtS2Fh9W5cQ4cOpXXr1rRu3RpwG6u/9dZbtX6/zp0788cfdVuWNHLkSC688MIanz9p0iQ23XTTanvevfvuuxt8zemnn75qE3cRkbp65BHXwEmjr1IXt956Kz169KBnz54EAgEikcgGz+/Tp896935tbHPmzGGXjXQo+/XXXzn++OM3eq3bbrutocISEVnLQw+5mVFe9KXwehsd7yQSEA5DcbG7pR8Kga9+exX5/X6+/vprli9fTosWLZgwYQIdOnRY67wb1vi/sGg0ytSpUzn88MPr9f6Nbf/992f8+PGNdv0VK1bQpEn2/kqKyNqWLoV773X7vvbOyl1as0uiIkF4dpji+cUE2wfrvY/gxx9/zPjx4/n8889p1qwZf/zxB2VlZQ0Y8folEgl89fz/inXZeuuta3Rz+LbbbuOaa65p8PcXEVmyxM2MOuooCASS//7ZOQKbSEBBAZx8MgwZ4r4WFLjj9XT44Yfz5ptvAvDCCy9w8smrlxfF43HOPPNM9txzT4LBIK+99hplZWXccMMNjBkzhkAgwJgxY/jzzz85+uij6dmzJ3vttRdffvklAIsWLaJv37706NGDQYMGUbWnxn333ccuu+zCLrvswgMPPLDO2EaMGEG3bt3Yc889mTJlyqrjCxcu5LjjjmOPPfZgjz32qPbcxqx5t/iee+7hxhtvXOu8adOmceCBB7L77rtTUFDA/PnzAXen+9JLL6V37948+OCDNX5fEckO//mPaxIxZIjXkUhjS1QkKBhdwMljT2bIpCGcPPZkCkYXkKioe26eP38+bdq0oVmzZgC0adOGrbfeGoCJEycSDAbZddddOfPMMyktLa322scee4wrrrhi1eOqM5dGjx7NnnvuSSAQ4F//+heJyv9/yMvL4/LLL2e33Xbj448/rna9aDTKXnvtRc+ePTnmmGNWLSOaNm0au+22G7vtthvDhg1b/XkkElxxxRXsscce9OzZk8cffxyonndHjhzJscceS79+/dhhhx248sorAbjqqqtYvnw5gUCAU045pc6fn4jIujz8sJtC7NXMqOwsYMNhiEQgFgNr3ddIxB2vp5NOOokXX3yRkpISvvzyS/Lz81c9d+utt3LwwQfz6aef8t5773HFFVdQXl7O0KFD6d+/P9FolP79+zNkyBCCwSBffvklt912G6eddhoAN910E/vttx/ffPMNxxxzDD///DPgkt+IESOIRCJ88sknPPHEExQXF1eLa/78+QwZMoQpU6bw4Ycf8u2336567pJLLuGyyy7js88+Y+zYsQwaNGidf7fJkydXm0L8/fff1+gzKS8v56KLLuLll19m2rRpnHnmmVx77bWrni8rK2Pq1KlcfvnlNfuQRSQrxGJwzz3Qrx/suafX0UhjC88OE5kXIVYWw2KJlcWIzIsQnl333Ny3b19++eUXunXrxvnnn8/7778PQElJCaeffjpjxozhq6++YsWKFTz66KPVXnvcccfxyiuvrHo8ZswYTjrpJKZPn86YMWOYMmUK0WgUn8/Hc889B7gb1fn5+XzxxRfst99+1a532mmnceedd/Lll1+y6667ctNNNwFwxhln8PDDD/PFF19UO/+pp55i00035bPPPuOzzz7jiSee4Mcff1zr7xiNRlf9PcaMGcMvv/zCHXfcQYsWLYhGo6tiExFpCEuXwn33QWEh9OrlTQzZOV+zuNhtWFRVPA7RqPtp1EPPnj2ZM2cOL7zwwlpTgt955x1ef/117rnnHsAl0JVFaFUffvghY8eOBeDggw9m0aJF/P3333zwwQeMGzcOgCOOOGLVGtoPP/yQY445Br/fD8Cxxx7L5MmTCVbpdhKJROjTpw8rN5rv378/3333HQDvvvtutYL277//JhaLkZeXVy2udU0hnjNnzkY/k5kzZ/L1119z2GGHAe6ucvv27Vc9379//41eQ0Syz2OPwR9/aO1rtiieX0y8rHpujpfFiS6IUtitbrk5Ly+PadOmMXnyZN577z369+/PHXfcQTAYpEuXLnTr1g2AgQMHMmzYMC699NJVr23bti3bbbcdn3zyCTvssAMzZsxg3333ZdiwYUybNo099tgDgOXLl7PlllsC4PP5OO6449aKY8mSJSxevJgDDzxw1fudcMIJLF68mMWLF3PAAQcAcOqppxKuvJn+zjvv8OWXX66aLrxkyRJmzZq1KuaVDjnkEDbddFMAdt55Z3766Se22WabOn1eIiIbkwp9KbKzgA0Gwe93t/dX8vsbbBL3kUceyeDBg5k0aRKLFi1addxay9ixY+nevXu18zfWUKKxVVRU8Mknn9C8efNav7ZJkyZUVFSselxSUrLWOdZaevTosdZ0qpVWFt4iIistWwZ33w2HHQZ77+11NJIMwfZB/Ll+YmWrc7M/10+gXaBe1/X5fPTp04c+ffqw6667MmrUqGo3eDfkpJNO4qWXXmLHHXfkmGOOwRiDtZaBAwdy++23r3V+8+bNG2zdq7WWhx9+mIKCgmrH17xxvHJ6NLi/64oVKxrk/UVE1hSLub4UoRBU3sPzRHZOIQ6FID8f8vLAGPc1P98dbwBnnnkmQ4YMYdc1ekoXFBTw8MMPr1q7unKab6tWrVi6dOmq8/bff/9VU34mTZpEmzZt2GSTTTjggAN4/vnnAQiHw6vWz+y///68+uqrLFu2jHg8ziuvvML+++9f7b3z8/N5//33WbRoEeXl5fz3v/9d9Vzfvn15+OGHVz2ORqM1/rtutdVW/P777yxatIjS0tJ1Nnnq3r07CxcuXFXAlpeX880339T4PUQk+zz+OPz+u0Zfs0moa4j8Dvnk5eZhMOTl5pHfIZ9Q17rn5pkzZzJr1qxVj6PRKJ06daJ79+7MmTOH2bNnA/Dss8+uGh2t6phjjuG1117jhRde4KSTTgLciOfLL7/M77//DsCff/7JTz/9tME4Nt10U1q3bs3kyZOrvd9mm23GZpttxocffghQbbpvQUEBjz76KOXl5QB89913xNecPbYBTZs2XfVaEZGG8J//wKJF3ufm7ByB9fmgqMiteY1G3chrA3QhXqljx45cfPHFax2//vrrufTSS+nZsycVFRV06dKF8ePHc9BBB3HHHXcQCAS4+uqrufHGGznzzDPp2bMnLVu2ZNSoUQAMGTKEk08+mR49erDPPvuw7bbbAtCrVy9OP/109qxcJDZo0KC17i63b9+eG2+8kb333pvNNtuMQJXR5oceeogLLriAnj17smLFCg444AAee+yxteJfuQZ2peuuu47jjz+eG264gT333JMOHTqw4447rvW63NxcXn75ZS6++GKWLFnCihUruPTSS+nRo0etP1sRyXzLl8Ndd8HBB8Maywglg/lyfBQNKCI8O0x0QZRAu0C9uxDHYjEuuugiFi9eTJMmTejatSvDhw+nefPmjBgxghNOOIEVK1awxx57cO655671+tatW7PTTjvx7bffrsqxO++8M7fccgt9+/aloqKCpk2bMmzYMDp16rTBWEaNGsW5557LsmXL2G677RgxYgTgGiyeeeaZGGPo27fvqvMHDRrEnDlz6NWrF9Za2rZty6uvvlrjv/s555xDz5496dWrl9bBiki9xeOuL0XfvrDXXt7GYqp2sk0XvXv3tmvu0TZ9+nR22mknjyKSxqCfqUh2eughuOQSmDQJ1jEo1iiMMdOstdqopx6Um9OHfi4iUlv33guDB8OUKbDPPsl5z/Xl5uycQiwiIimppATuvNMVrskqXkVERGT9VvalOOSQ5BWvG5KdU4hFRCQlPfUU/PorPPus15GIiIgIwPDh8NtvUKWFjqc0AisiIimhtBTuuMOtez3oIK+jERERkeXL3cyoPn1gjR6xnsmoEVhrLcYYr8OQBpCOa7NFpH5GjIC5c+Hpp12DeMkMys2pRflVRGrjySdhwQJ44QWvI1ktY0ZgmzdvzqJFi/QPcwaw1rJo0aI67UsrIumprAxuv93t+XrooV5HIw1FuTm1KL+KSG2UlLiZUfvvn1p9KTJmBLZjx47MnTuXhQsXeh2KNIDmzZvTsWNHr8MQkSQZNQp+/tnt/6rBusyh3Jx6lF9FpKaeftr1pXjmmdTKzRlTwDZt2pQuXbp4HYaIiNRSeTncdhvsuScUFHgdjTQk5WYRkfRUWupmRu2zj9uXPZVkTAErIiLp6dlnYc4ceOSR1LrDKyIikq1GjnR9KZ56KvVyc8asgRURkfRTXg633gq77w6HH+51NCIiIlJW5mZG5efDYYd5Hc3aNAIrIiKeef55+OEHeO211LvDKyIiko2eecb1pXjssdTMzRqBFRERT6xYAbfcAoEA/OMfXkcjIiIiK/tS9O4N/fp5Hc26aQRWREQ88eKLMHs2jBuXmnd4RUREss3o0fDjj/DQQ6mbmzUCKyIiSZdIuNHXnj3hqKO8jkZERERWrHB9KXr1giOO8Dqa9dMIrIiIJN1LL8HMmfDyy5CjW6kiIiKee/55+P57ePXV1B19BY3AiohIkiUScPPNsMsucMwxXkcjIiIiK/tS7LYbHHmk19FsmEZgRUQkqcaOhenTYcwYjb6KiIikgjFjYNYsl6NTefQVVMCKiEgSVVTA0KGw005w3HFeRyMiIpLdEhUJxn8X5pJxxWx7WJB/HBkCfF6HtUEqYEVEJGleeQW++cats/HVNT8mEhAOQ3ExBIMQCtXjYiIiItkpUZGgYHQBH86JULpLnLjPT+j5fIoGFOHLSd28qgJWRESSIpGAm26C7t3hxBPrcZGCAohEIB4Hvx/y86GoSEWsiIhILYRnh4nMjVBqY5ADJTZGZF6E8Owwhd0KvQ5vvbT6SEREkuLZZ+Grr1wRW+daMxx2xWssBta6r5GIOy4iIiI1Vjy/mFh5vNqxeFmc6IKoNwHVUNIKWGPM08aY340xX6/jucuNMdYY0yZZ8YiISPIsWwbXXQd77lmP0Vdw04bj1ZMt8ThEo/UJL2spN4uIZK8dNwtCmb/aMX+un0C7gDcB1VAyR2BHAv3WPGiM2QboC/ycxFhERCSJ7r8f5s2De+6pZ3fDYNBNG67K74dAoD7hZbORKDeLiGSlr14Jwdx8WvjyMBjycvPI75BPqGvI69A2KGlrYK21HxhjOq/jqfuBK4HXkhWLiIgkz++/w513wtFHw/771/NioZBb87rmGthQaifbVKXcLCKSnebPh3vv9nHc4UWcfmKY6IIogXYBQl1DKd3ACTxu4mSMOQqYZ639wqT6hkMiIlInN93kphDfcUcDXMzncw2bwmE3bTgQUBfiBqbcLCKS+YYMgfJyuON2H127FqZ006Y1eVbAGmNaAtfgpijV5PxzgHMAtt1220aMTEREGsqMGfD443Duua77cIPw+aCw0P2RBqXcLCKS+b75Bp56Ci66CLp29Tqa2vOyC/H2QBfgC2PMHKAj8Lkxpt26TrbWDrfW9rbW9m7btm0SwxQRkbq66ipo2dLd6ZW0oNwsIpLhrrwSWrWC66/3OpK68WwE1lr7FbDlyseVibK3tfYPr2ISEZGG88EH8NprcOutoNomPSg3i4hktokT4a234K67YIstvI6mbpK5jc4LwMdAd2PMXGPMWcl6bxERSa6KChg8GDp2hEsv9ToaWR/lZhGR7FFRAVdcAZ06uenD6SqZXYhP3sjznZMUioiINLKXXoLPPoORI90UYklNys0iItnjuefcduqjR0Pz5l5HU3deroEVEZEMVFoKV18Nu+0GAwZ4HY2IiIgsXw7XXgu77w4nb/DWZerzdBsdERHJPI88AnPmwIQJ2t1GREQkFTz4IPzyCzzzDOSk+RBmmocvIiKp5M8/4ZZboF8/OPRQr6MRERGRhQvh9tvd7nN9+ngdTf2pgBURkQZzyy3w99+uu6GIiIh4b+hQiMczJzergBURkQbxww9u+vAZZ8Cuu3odjYiIiHz3HTz2GAwaBDvt5HU0DUMFrIiINIhrroGmTd2dXhEREfHe1Ve7jsM33uh1JA1HBayIiNRbJAJjxsDll8PWW3sdjYiIiHz4IYwbB1deCe3aeR1Nw1EBKyIi9WItDB4MW23lNkgXERERb1nrcnL79vB//+d1NA1L2+iIiEi9vPaau8v72GPQqpXX0YiIiMjLL8Mnn8BTT4Hf73U0DUsjsCIiUmfl5W5q0k47wVlneR2NiIiIlJbCVVe5hooDB3odTcPTCKyIiNTZ8OEwaxa88QY0UUYRERHx3KOPup0BwmHw+byOpuFpBFZEROpkyRLX1bBPHzjiCK+jERERkb/+crsBHHYYFBR4HU3jUAErIiJ1cued8McfcM89YIzX0YiIiMhtt8HixXD33Zmbm1XAiohIrf3yC9x/P5xyCuy+u9fRiIiIyJw58NBDcNppsNtuXkfTeFTAiohIrV13nWvRf+utXkciIiIiANdcAzk5cMstXkfSuFTAiohIrUSj8OyzcPHF0KmT19GIiIjI1Knwwgtuz9eOHb2OpnGpgBURkRqzFgYPhtat3Z1eERER8dbK3Ny2Lfz7315H0/i06YGIiNTY22/DxInwwAOw2WZeRyMiIiJvvAHvvw/DhsEmm3gdTePTCKyIiNRIIgFXXgnbbw/nned1NCIiIrJihRt17dYNzj7b62iSQyOwIiJSIyNHwtdfw3//C7m5XkcjIiIiTz4JM2bAq69C06ZeR5McKmBFRGSj4nG4/nrYe2847rg6XiSRgHAYioshGIRQCHy+Bo1TREQk0yUqEoRnh/nkp2IeeTrIfvuHOPLI7MmnKmBFRGSj7r0X5s+Hl1+u48boiQQUFEAk4qphvx/y86GoSEWsiIhIDSUqEhSMLiAyL0KsNA6H+inZOp8KW4TPZEc+1RpYERHZoAUL4K673MjrPvvU8SLhsCteYzHXLjEWc4/D4QaNVUREJJOFZ4dd8VoWA2OhWYwZSyOEZ2dPPlUBKyIiGzR4MJSVwe231+MixcVu5LWqeNxtKisiIiI1Ujy/mHhZ9XwaL4sTXRD1JiAPqIAVEZH1mjgRnnsOrroKdtihHhcKBt204ar8fggE6hOeiIhIVgm2D9Isp3o+9ef6CbQLeBOQB1TAiojIOpWUuO1ytt8err66nhcLhdya17w8t4g2L889DoUaJFYREZFs0KdDCOblY8rzMBjycvPI75BPqGv25FM1cRIRkXW6806YNcv1WWrRop4X8/nchcJhN204EFAXYhERkVq643YfJcOLuPWlMBVbRgm0CxDqGsKXkz35VAWsiIisZdYsuO02OOkk6Nu3gS7q80FhofsjIiIitTJ9umuqeOoAH9ccXwhkZz7VFGIREanGWjj/fGjeHO67z+toRERExFq3rMfvh3vu8Toab2kEVkREqnnxRXj3XXjkEWjf3utoRERE5Nln4f334fHHYcstvY7GW0kbgTXGPG2M+d0Y83WVY3cbY2YYY740xrxijNksWfGIiMjaFi+Gyy6D3r3h3HO9jkYam3KziEjq+/NPuPxy2GsvGDTI62i8l8wpxCOBfmscmwDsYq3tCXwH1LfPpYiI1MO118LChe4Or/orZYWRKDeLiKS0q66Cv/5yuTlHC0CTV8Baaz8A/lzj2DvW2hWVDz8BOiYrHhERqe7TT+HRR+Gii6BXL6+jkWRQbhYRSW1TpsATT7jZUT17eh1NakilGv5MIOx1ECIi2WjFCvjXv9ya16FDvY5GUohys4iIR8rL3XKebbaBIUO8jiZ1pEQTJ2PMtcAK4LkNnHMOcA7Atttum6TIRESywyOPuO1Z//tf2GQTr6ORVKDcLCLirQcegK+/hldfhbw8r6NJHZ6PwBpjTsdtYnSKtdau7zxr7XBrbW9rbe+2bdsmLT4RkUw3dy5cfz2EQnDccV5HI6lAuVlExFs//QQ33ghHHglHHeV1NKnF0xFYY0w/4ErgQGvtMi9jERHJVpde6qYQDxsGxngdjXhNuVlExHsXX+y+PvSQt3GkomRuo/MC8DHQ3Rgz1xhzFvAI0AqYYIyJGmMeS1Y8IiICb74JY8fCDTdAly5eRyPJptwsIpJ6XnsNXn/djcB26uR1NKnHbGBmUMrq3bu3nTp1qtdhiIiktWXLoEcPaNHCrX/NzfU6Iu8YY6ZZa3t7HUc6U24WEam/WAx23hk22wymTYOmTb2OyDvry80p0cRJRESS7+abYc4ceP/97C5eRUREUsWNN8Ivv8CLL2Z38bohnjdxEhGR5PvmG7jnHjj9dDjgAK+jERERkS++cJ2Hzz4b9tnH62hSlwpYEZEsU1Hh9pXbZBO4+26voxEREZGVubl1a7jjDq+jSW2aQiwikmVGjYIPP4SnnoI2bbyORkRERJ58Ej75xOXozTf3OprUphFYEZEs8scfcMUVsN9+bvqwiIiIeOv33+Hf/4Y+feDUU72OJvWpgBURySL//jcsWQKPPgo5ygAiIiKeGzwY4nGXm7Uf+8bpf19ERLLE5Mnw9NNw+eWwyy5eRyMiIiL/+x88+6y7wbzjjl5Hkx5UwIqIZIGyMtcconNnuOEGr6MRERGR0lI47zzYbju45hqvo0kfauIkIpIF7rsPvv0Wxo+Hli29jkZERETuugu++w7CYWjRwuto0odGYEVEMtyPP8LQoXDssXDEEV5HIyIiIrNnw623woknQr9+XkeTXjQCKyKSwayFCy8Enw8efLAOF0gk3K3h4mIIBiEUchcTERGROrEWLrgAcnPh/vu9jib9qIAVEclgr7wCb73lphB37FjLFycSUFAAkYhrj+j3Q34+FBWpiBUREamlREWC8Owwz0wo5p05QR64JcTWWyuf1pYKWBGRDLV0KVx8MQQCcNFFdbhAOOyK11jMPY7F3ONwGAoLGzJUERGRjJaoSFAwuoBP5kaIl8bJOdHP65vmc2FFEb4cFbG1oTWwIiIZ6vrr4ddf4bHHoEldblcWF7uR16ricYhGGyI8ERGRrBGeHSYyL0K8PAY5loomMT79NUJ4dtjr0NKOClgRkQw0aRI89BCcf76b9VsnwaCbNlyV3++GdEVERKTGiucXEy+rflM4XhYnuiDqTUBpTAWsiEiGWbIEBg6Erl3hzjvrcaFQyFW/eXlgjPuan++Oi4iISI3t0CoI5dVvCvtz/QTaBbwJKI1pDayISIa55BKYNw+mTFl7ALVWfD7XsCkcdtOGAwF1IRYREamD1+4NYZvm06JrhJJEHH+un/wO+YS66qZwbamAFRHJIOPGwahRbv1rnacOV+XzuYZNatokIiJSJ2PGwIvP+xhyYxG9TwwTXRAl0C5AqGtIDZzqQAWsiEiGWLAAzjkHdt/dFbAiIiLirXnz4LzzYM894bprfTRpUkhhN90Urg+tgRURyQDWwqBBrknws89C06ZeRyQiIpLdrIWzzoKSEpeb67QjgKxFH6OISAZ48kl480148EHYaSevoxEREZFHH3WtJIYNg27dvI4mc2gEVkQkzX3/PVx2GRxyCFx4odfRiIiIyHffweDBUFDgphBLw1EBKyKSxhIJOO00Ny1pxAjI0b/qIiIinlqxAk49FZo3h6efdjvRScPRFGIRkTR2113w0Ufw3HOwzTZeRyMiIiK33Qaffuq6D2+9tdfRZB7dqxcRSVPRKAwZAieeCCef7HU0IiIiMnUqDB0K//yny8/S8FTAioikoZISGDAA2rRxTSI0PUlERMRby5a5qcPt2sEjj3gdTebSFGIRkTR07bXwzTfw9tuw+eZeRyMiIiJXXQUzZsCECdC6tdfRZC6NwIqIpJlJk+D+++H88113QxEREfHWhAnw8MNw8cVw6KFeR5PZVMCKiKSRJUtg4EDo2tU1cBIRERFv/fUXnHEG7Lgj3HGH19FkPk0hFhFJI5dcAvPmwZQp4Pd7HY2IiIhceCH89hu8+iq0aOF1NJlPI7AiImli3DgYNcqtf83P9zoaERERGTMGnn8ebrgBevf2OprskLQC1hjztDHmd2PM11WObW6MmWCMmVX5VcudRUTWYcECOOcclxyvu87raCRTKDeLiNTdvHlw3nnupvLVV3sdTfZI5gjsSKDfGseuAiZaa3cAJlY+FhGRKqyFQYMgHodnn4WmTb2OSDLISJSbRURqzVo46yy3rd0zz0ATLcxMmqQVsNbaD4A/1zh8FDCq8vtRwNHJikdEJF08+SS8+SbceadrECHSUJSbRUTq5tFHoagI7rkHunXzOprs4vUa2K2stfMrv18AbLW+E40x5xhjphpjpi5cuDA50YmIeOz77+Gyy1xL/gsv9DoayRLKzSIiGzBzJgwe7LayO+88r6PJPl4XsKtYay1gN/D8cGttb2tt77Zt2yYxMhERbyQScNppbsrwiBGQkzL/Yku2UG4WEaluxQqXm1u0gKefBmO8jij7eD1b+zdjTHtr7XxjTHvgd4/jERFJGXfdBR99BM89Bx071vLFiQSEw1BcDMEghELg8zVKnJJxlJtFRNbjttvg009d9+Gtt/Y6muzkdQH7OjAQuKPy62vehiMikhqKi2HIEDjxRDj55Fq+OJFw85oiEdf5ye93LRKLilTESk0oN4uIrMNnn8HQofDPf7r8LN5I5jY6LwAfA92NMXONMWfhkuNhxphZwKGVj0VEslpJCQwYAG3auCYRtZ6eFA674jUWc20SYzH3OBxulHglfSk3i4hsXKIiwdivxhO642Y223M8Dz6U8DqkrJa0EVhr7frGEA5JVgwiIung2mvh22/h7bdh883rcIHiYjfyWlU8DtEoFBY2RIiSIZSbRUQ2LFGRoGB0AR/8EKF8lzgtAn5OejOfogFF+HI0q8kLagkiIpJC3nkH7r8fzj/fzQKuk2DQTRuuyu+HQKC+4YmIiGSV8OwwU36KUG5ikGNZXhEjMi9CeLZmNXlFBayISIr46Se33nWXXeDuu+txoVDIrXnNy3Pzj/Py3ONQqMFiFRERyQYTvy2mZEX1WU3xsjjRBVFvAhLPmziJiAhu3evxx7v2/GPHQsuW9biYz+caNoXDbtpwIKAuxCIiIrW0fDm89lgQdvdDbmzVcX+un0C7gHeBZTkVsCIiKeCSS2DqVHj1Vdhhhwa4oM/n1rtqzauIiEitWQsXXAA/vhNit775fF8aIV4Wx5/rJ79DPqGumtXkFRWwIiIeGzkShg+Hq66Co47yOhoRERF58kkYMQKuu87HjZcWEZ4dJrogSqBdgFDXkBo4eUgFrIiIh6JROO88OPhguPlmr6MRERGRqVPhwguhb1+48Ubw5fgo7FZIYTfNakoFauIkIuKRv/6CY4+FLbaAF16AJrqlKCIi4qlFi1xPinbt4Lnn1D4iFel/l0REPFBRAaeeCnPnwgcfwJZbeh2RiIhIdksk4JRTYP58+PBDaNPG64hkXVTAioh44Lbb4M034ZFHYK+9vI5GREREhg51Tfwfewz22MPraGR9NIVYRCTJ3nkHbrjB3eU9/3yvoxEREZG33nIF7MCBcM45XkcjG6ICVkQkiX76Cf75T+jRAx5/HIzxOiIREZHs9uOPMGAA7LYb/Oc/ys2pTgWsiEiSlJbCCSdAeTmMGwd+v9cRiYiIZLfly+G441xvirFjoWVLryOSjdEaWBGRJLnkEvjsM3jlFdhhB6+jERERkQsvhOJieOMN2H57r6ORmtAIrIhIEowa5aYM//vfcPTRXkcjIiIiTz4JTz8N114LhdriNW2ogBURaWTRKJx7Lhx0ENxyi9fRiIiIyLRpbvT1sMPgppu8jkZqo9YFrDHGb4zRlr4iIjXw119ubc0WW8CLL0ITLdyQRqDcLCJSc4sWudy85Zbw/PPg07+eaWWj/ytljMkBTgJOAfYASoFmxpg/gDeBx621sxs1ShGRNFRRAaedBr/8Au+/7xKlSENQbhYRqZtEwnUcnj8fJk+GNm28jkhqqyYjsO8B2wNXAe2stdtYa7cE9gM+Ae40xgxoxBhFRNLS7bfD+PFw332w995eRyMZRrlZRKQObr4Z3n4bHnwQ9tzT62ikLmoyme1Qa225MeZza22vlQettX8aY+ZZa48zxjRtxBhFRNLOhAlw/fVw8slwwQVeRyMZSLlZRKSWwmEYOhROPRX+9S+vo5G6qkkBe4wxphfQyhizEzDTWltR+dxwoKe1trzRIhQRSTM//+wK1513hieeqOWG6ImEy7DFxRAMQiikxTmyLsrNIiK18OOPcMopsOuu8NhjtczNklJqUsBOAZoDg4D7gO7GmMXAr8DyxgtNRCT9lJbC8cdDWRmMGwd+fy1enEhAQQFEIhCPuxfn50NRkYpYWZNys4hIDZWUuNxcUQFjx0LLll5HJPWx0QLWWjsPeMYY8721dgqAMWYLoDMwo3HDExFJL5deCp995hJkt261fHE47IrXWMw9jsXc43BYG9RJNcrNIiIblqhIEJ4dpnh+MZP/G+Tz4hCvveqja1evI5P6qkkXYmOdKSuPWWsXAYvWPKeRYhQRSQvPPOOmJV1xBRx7bB0uUFzsRl6risfdRrIqYKUK5WYRkfVLVCQoGF1AZF6EeFkcu4mfTtflc0RhEaAZTemuRl2IjTEXGWO2rXrQGJNrjDnYGDMKGNg44YmIpIcvvnANIfr0gdtuq+NFgsG15xz7/RAI1DM6yUDKzSIi6xGeHSYyL0KsLIbFQrMYfzSPEJ4d9jo0aQA1KWD7AQngBWPMr8aYb40xPwCzgJOBB6y1IxsxRhGRlDZ/Phx5JGy+Obz4IjSpSXeBdQmF3JrXvDzXXSIvzz0OhRo0XskIys0iIutRPL+YeFn1GU3LyuNEF0S9CUgaVE3WwJYA/wH+U9mSvw2w3Fq7uJFjExFJebGYm927aBG8/z5stVU9LubzuYZN4bCbNhwIqAuxrJNys4jI+nXfNIgp92ObxlYd8+f6CbQLeBeUNJgajxMYYyYAg621XzRiPCIiaWPFCujf39Wab7wBu+/eABf1+VxFrDWvUgPKzSIi1ZWXw5P/DlHRLp/mXSOUVsTx5/rJ75BPqKtmNGWC2kx0+zfwgDFmDnCNtXZ+44QkIpL6rIULLoC33oLHH4fDD/c6IslSys0iIpWshXPPhQnv+Bj+RBHtDwgTXRAl0C5AqGsIX45mNGWCGhew1trPgYOMMccBbxtjxgF3WWu135yIZJ0774Thw+Gqq+Ccc7yORrKVcrOIyGq33AJPPw3XXw9nD/IBhRR204ymTFOTJk6rGGMMMBN4FLgImGWMObUxAhMRSVXPPw9XXw0nnwy33up1NJLtlJtFRGDUKLjhBjj1VLjpJq+jkcZU4wLWGDMFmAfcD3QATgf6AHsaY4bXJwhjzGXGmG+MMV8bY14wxjSvz/VERBrL++/DGWfAgQfCiBGQU6vbgCINS7lZRAQmToRBg+Cgg+DJJ10jf8lctVkDew7w7To2Rb/IGDO9rgEYYzoAFwM7W2uXG2NeAk4CRtb1miIijWH6dDj6aNh+e3jlFWjWzOuIRJSbRSS7ff01HHssdO8O48ZBbq7XEUljq80a2G828PQRDRBHC2NMOdAS+LWe1xMRaVALFrgdbZo1c42bWrf2OiIR5WYRyW7z5rnc7Pe73LzZZl5HJMnQIJPfrLU/1OO184B7gJ+B+cASa+07a55njDnHGDPVGDN14cKFdQ9WRKSWYjE44ghYuBDGj4fOnb2OSGTjlJtFJJMtXepy8+LF8OabsO22XkckyeL56i1jTGvgKKALsDXgN8YMWPM8a+1wa21va23vtm3bJjtMEclSK1bASSe5vV7HjIHevb2OSKTxKTeLSCorL4cTTnDTh//7XwgGvY5IksnzAhY4FPjRWrvQWlsOjAP28TgmERGshYsvdnd2hw2DQnXil+yh3CwiKclaOP98KCqCRx+Ffv28jkiSLRUK2J+BvYwxLSu3AjgEqHPjCRGRhnL33S45Xnml2xhdJIsoN4tISrrtNtdp+Npr4eyzvY5GvOB5AWutjQAvA58DX+FiqlfrfxGR+nrxRfj3v6F/f7j9dq+jEUku5WYRSUWjR8N118GAAXDzzV5HI16pzTY6jcZaOwQY4nUcIiIAH3wAAwfC/vvDyJHa61Wyk3KziKSS//0PzjwT+vSBp57SXq/ZLCUKWBGRVDFjhtvrtUsXePVVaN68Bi9KJCAchuJi10kiFAKfr5EjFRERyQ7ffOP2et1hB+31KipgRURW+e03V3s2ber2k9t88xq8KJGAggKIRCAed5vR5ee77hIqYkVEROrl119dbm7RQvuwi6OJcSIiuNqzsNAVsePHw3bb1fCF4bArXmMx1xoxFnOPw+FGjVdERCTTLV3qcvOff7odATp18joiSQUqYEUk663c6/Xzz13zpj32qMWLi4td9VtVPO42jhUREZE6WbECTjwRvvzS7fXaq5fXEUmqUAErIlnNWrjkEjfq+tBDcOSRtbxAMOimDVfl90Mg0FAhioiIZI1ERYI3Zo5nr3/fzNs/jOeRYQlCIa+jklSiNbAiktXuvRf+8x8YPBguuKAOFwiF3JrXNdfAKtuKiIjUSqIiQcHoAib/GKEsL07Tk/283CKfsyuK8OWor4Q4KmBFJGs9+SRccQWccALceWcdL+LzuYZN4bCbNhwIqAuxiIhIHYRnh/nwxwhlxCAHyokRmRchPDtMYbdCr8OTFKECVkSy0ogRcM450K8fPPNMPfd69flcl4lCJVcREZG6emJ8MaUV8WqLHONlcaILoipgZRWtgRWRrPPss3DWWXDooW4/uRrt9SoiIiKN5umn4fXHg/hs9b4S/lw/gXYBb4KSlKQCVkSyyvPPw+mnw0EHwauvun3lRERExDvPPAODBsGhnUIcsF0+ebl5GAx5uXnkd8gn1FV9JWQ1TSEWkawxZgyceioccAC8/jq0bOl1RCIiItntuefcjeWDD4bXX/OR26yI8Oww0QVRAu0ChLqG1MBJqlEBKyJZ4eWX4ZRTYN994Y031t75RkRERJLrxRfhtNPgwAPdjWU3K8pHYbdCrXmV9dIUYhHJeK+8Aief7Ha3efNNyMvzOiIREZHs9t//woAB7sby+PGaFSU1pwJWRDLa66/DiSdC795up5tWrbyOSEREJLuNG+duLO+1l7uxrFlRUhsqYEUkY735Jhx/PASD8PbbsMkmXkckIiKS3V57Dfr3hz32gLfe0o1lqT0VsCKSkd5+G449Fnr2hHfegU039ToiERGR7DZ+PJxwAvTqpRvLUncqYEUk40yYAEcfDTvv7IrXzTbzOiIREZHsFg7Dcce5G8tFRbqxLHWnAlZEMsr//gdHHgndu8O778Lmm3sdkYiISHZ75x045hjo0UM3lqX+VMCKSMZ4/30oLISuXV3xusUWXkckIiKS3SZOhKOOcjeWJ0zQjWWpP+0DKyIZYfJkOOII6NLFJcu2bWvwokTCzWkqLnadnkIh8GmzdBERkYYwaRL84x+6sSwNSwWsiKS9KVPg8MOhY0dXvG65ZQ1elEhAQQFEIhCPux7++fluYY6KWBERkXr54IM63FgWqQFNIRaRtPbJJ27gtH17t/61XbsavjAcdsVrLAbWuq+RiDsuIiIidbbyxvK227rcXKMbyyI1pBFYEUlbn33mBlG33BLeew+23roWLy4udiOvVcXjEI26hbQiIiJSI4mKBOHZYYrnF9NiSZAbB4To0MHH//4HW23ldXSSaVTAikhamjYN+vZ162neew86dKjlBYJBN204Flt9zO+HQKAhwxQREcloiYoEBaMLiMyLEC+LY8v8tPhnPhOuLaJ9ey3JkYanKcQiknY+/xwOO8ztIffee7DNNnW4SCjk1rzm5YEx7mt+vjsuIiIiNRKeHSYyL0KsLIbFQm4Ms02EL5drSY40Do3AikhamTjR7SXXurUrXjt1quOFfD7XsCkcdtOGAwF1IRYREaml4vnFxMuqL8lZviJOdEGUwm5akiMNTwWsiKSNF1+E005ze8mFw67rcL34fG69q9a8ioiI1MniGUFsqR+arV6S48/1E2gX8C4oyWiaQiwiaeGBB+Dkk2GvvVxr/noXryIiIlJn1sI998B954XYLJ6Pv2keBkNebh75HfIJddWSHGkcGoEVkZRWUQFXXQV33w3HHQejR0Pz5l5HJSIikr0qKmDwYLj/fjjxRB8j7ijif7+EiS6IEmgXINQ1hC9HS3KkcaiAFZGUVVYGZ53litYLLoAHH9QSVRERES+VlsLpp7tlPRdf7IrYnBwfhd0KteZVkiIlphAbYzYzxrxsjJlhjJlujNnb65hExFtLl8I//uGK11tvhYcfVvEqkkzKzSKypr//hsMPd8XrnXe65T05KVFNSDZJlRHYB4G3rbXHG2NygZZeByQi3vntNzjiCNcc+Omn4YwzvI5IJCspN4vIKvPnu+L1669h1CjXVFHEC54XsMaYTYEDgNMBrLVlQJmXMYmId77/HgoK4Ndf4bXXXCErIsml3CwiVX33ncvNCxfCG29Av35eRyTZLBUG/bsAC4ERxphiY8yTxhj/micZY84xxkw1xkxduHBh8qMUkUY3bRrssw8sXuz2eFXxKuIZ5WYRASASgX33hXjc5WYVr+K1VChgmwC9gEettUEgDly15knW2uHW2t7W2t5t27ZNdowi0sjeeQcOPBBatoQpUyA/3+uIRLKacrOIEA7DwQdDq1YuN++xh9cRiaRGATsXmGutjVQ+fhmXNEUkSzz3nBtt7doVPvoIunffwMmJBIwfDzff7L4mEkmLUySLKDeLZLlRo1wzxe7dXW7eYQevIxJxPF8Da61dYIz5xRjT3Vo7EzgE+NbruEQkOe691+0ld/DBMG4cbLrpBk5OJNwinEjEzWXy+91QbVGRWhSLNCDlZpHsZa3rMHz11XDooS43t2rldVQiq6XCCCzARcBzxpgvgQBwm7fhiEhjq6iA//s/V7z27w9vvbWR4hXcXKZIBGIxl2FjMfc4HE5KzCJZRrlZJMskEnDJJa54/ec/4c03VbxK6vF8BBbAWhsFensdh4gkR2mp2xrnhRdcorzvvhruI1dc7EZeq4rH3X47hdo8XaQhKTeLZJeSErc1zn//C5dfDnfdpT1eJTWlRAErItnj77/h2GNh4kSXHAcPBmNq+OJg0E0bjsVWH/P7IRBojFBFREQyWqIiQXh2mI/nFPP68CBfjwtxzz0+Lr/c68hE1k8FrIgkzYIFEAq5TdCfeQZOPbWWFwiF3JrXNdfAhkKNEq+IiEimSlQkKBhdwCe/RIiXx6Gbnx535HPpZUWA+kpI6lIBKyJJ8cUXcPTRbhP08eNdL6Za8/lcw6Zw2E0bDgRc8aoGTiIiIrUSnh3mo58jLE/EwADNYvyUiBCeHaawm5blSOrSzHYRaXTPPQd77w3l5TBpUh2L15V8Prfe9brr3FcVryIiIrU2IlzM8vLqfSXiZXGiC6LeBCRSQypgRaTRlJfDpZfCgAFu8/Np06C3WsKIiIh4pqwMLrwQxv0nSE6Fv9pz/lw/gXYBbwITqSEVsCLSKH77DQ47DB580BWx774LW23ldVQiIiLZa8ECOOQQGDYMLvtHiAO75pOXm4fBkJebR36HfEJd1VdCUpvWwIpIg/v0U9dp+M8/YfRoOOUUryMSERHJbh9/DMcdB0uWuG3sTjrJR6KiiPDsMNEFUQLtAoS6hvDlaGmOpDYVsCLSoJ58Ei64ADp0gI8+0g43IiIiXrIWHn8cLr4YttkG3n4bevZ0z/lyfBR2K1TTJkkrmkIsIg2itBT+9S84+2zo0wemTlXxKiIi4qWSEhg0CM47Dw491OXmlcWrSLpSASsi9TZvnitahw+Hq6+Gt96CzTf3OioREZHs9csvcMAB8PTTrnH/G29A69ZeRyVSf5pCLCL1MnkynHACxOMwdqxb+7pOiYTbv7W4GIJB7d8qIiLSSCZNghNPdCOwr7zi9mEXyRQqYEWkTqyFhx+Gyy+H7baD//0Pdt55PScnEm7z10jEVbp+P+TnQ1GRilgREZEGYi088ABccQXssIMrXnfc0euoRBqWphCLSK0tWwYDB8Ill8Dhh7uuw+stXsGNvEYiEIu57BqLucfhcNJiFhERyWTxuOv6/3//B0cd5XKzilfJRCpgRaRW5syBffd12+PcfLO7u7vppht5UXGxy6xVxeMQjTZSlCIiItnj++9hn33gxRfhttvg5ZehVSuvoxJpHJpCLCI1NmECnHQSVFTA+PFu9LVGgkE3bTgWW33M71ebYhERkVpKVCQIzw5TPL+YYPsg5vsQA/7pwxg3samgwOsIRRqXClgR2Shr4a674JproEcPGDcOunatxQVCIbfmdc01sKFQo8UsIiKSaRIVCQpGFxCZFyFeFqcpfsp+yKdnpyJefcVHly5eRyjS+FTAisgGLV0KZ57ppiP17w9PPeXqz1rx+VzDpnDYTRsOBNSFWEREpJbCs8NE5kWIlbkZTWXEaNI5wg2Dw3TpUuhxdCLJoQJWRNZryhQ47TS37vWee1xjCGPqeDGfDwoL3R8RERGpteL5xcTLqveUSOTEmf5XFFB+leygJk4ispayMrj2WrcBurXw/vtuu5w6F68iIiJSL6WlEH07iC2tPg3Kn+sn0C7gTVAiHlABKyLVTJ8Oe+/tuhiefjp88QXst5/XUYmIiGSvr792rSPG3Rlia/LxN83DYMjLzSO/Qz6hruopIdlDU4hFBHCdhR95BP79b8jLc9vjHH2011GJiIhkr4oKePBBuPpqt2XdG6/7CB1eRHh2mOiCKIF2AUJdQ/hy1FNCsocKWBFh3jw44wy3Tc4RR7hGTVtt5XVUIiIi2euXX9xMqP/9D448Ep54ArbcEsBHYbdCCrtpzatkJ00hFslyY8bArru6hk2PPw5vvKHiVURExEvPP+9ycyQCTz4Jr766sngVEY3AimSpxYvhwgvhuefcuppnn4UddljjpETCbX1TXAzBoLa+ERERaUR//QXnnw8vvgj77APPPAPbb+91VCKpRQWsSBZ67z0YOBB+/RWGDnVra5qs+a9BIgEFBe72bzzuNn/Nz3f7uaqIFRERaVDvvuumDP/2G9x6q+tJoXQrsjZNIRbJIiUlbjucgw+GFi3g44/h+uvXUbyCG3mNRCAWc3vpxGLucTic9LhFREQy1fLlcOmlcNhhsMkmLtVec42KV5H1UQErkiW++AL22APuuw8uuMDNCt5jjw28oLjYjbxWFY9DNNqYYYqIiGSN4mLo3dt1Gr74Ypg2DXr18joqkdSmKcQiGS6RcEXrddfB5pu7AdR+/WrwwmDQTRuOxVYf8/shEGisUEVERDJWoiJBeHaY4vnF9NwqyNfjQtx0o4+2beGdd9wIrIhsnApYkQw2Z45b6/rBB3Dssa7LcJs2NXxxKOTWvK65BjakzdJFRERqI1GRoGB0AZF5EeJlccwKPxU/53P8MUU8/qiPzTf3OkKR9KECViQDWQujRrnpSAAjR8Jpp4ExtbiIz+caNoXDbtpwIKAuxCIiInUQnh0mMjdCrNzNarJNYjTbPsJp/cNsvrn2cxWpjZRZA2uM8Rljio0x472ORSSdffstHHQQnHGGqzm//NKNwtaqeF3J54PCQjf/uLBQxatIllFuFmkYb39RTKysel+JMhvni9+i3gQkksZSpoAFLgGmex2ESLqKx+Gqq2C33VzROnw4TJoEnTt7HZmIpDHlZpF6iMXgiivg0RuCUO6v9pw/10+gXcCbwETSWEoUsMaYjsARwJNexyKSbqyFV1+FnXeGO+90U4VnzoSzz4aclPgvXETSkXKzSN1ZC2PHwk47wT33wBn7h9h/u3zycvMwGPJy88jvkE+oq/pKiNRWqqyBfQC4Emi1vhOMMecA5wBsu+22yYlKJMX9+KNb5zp+POy6K3z4Iey7r9dRiUiGeADlZpFa+/57uOgi10Jit93gpZdg7719JCqKCM8OE10QJdAuQKhrCF+OluaI1Jbn4zPGmELgd2vttA2dZ60dbq3tba3t3bZt2yRFJ5KaSkvh1lvdqOukSXDvvW7vOBWvItIQlJtFaq+0FG6+GXbZxd1QfuABmDoV9t7bPe/L8VHYrZDrDriOwm6FKl5F6igVRmD3BY40xhwONAc2McaMttYO8DgukZQ0cSJccIGbJnz88XD//dCx4xonJRLu1m9xsdvPVd2DRaR2lJtFamHCBJebZ82C/v3d/utbb+11VCKZyfMC1lp7NXA1gDGmDzBYCVJkbfPnw+WXwwsvwPbbu/q0X791nJhIQEHB2vu3FhWpiBWRGlFuFqmZX3+Fyy5z04R32AHeeQcOO8zrqEQym+dTiEVkw1asgIcfhh13hHHjYMgQ+Prr9RSv4CrbSMS1PrTWfY1E3HERERGptxUr3BThHXeE1193U4e/+krFq0gyeD4CW5W1dhIwyeMwRFJGJALnnedmAvftC4884u7wblBxsRt5rSoeh2jU7eUqIlILys0i1X38scvNX3zhVug88ghst53XUYlkj5QqYEXE+fNPuOYat5dr+/ZuatLxx4MxNXhxMOimDcdiq4/5/RAINFa4IiIiGSlRkSA8O0zx/GK29weZ+HiIp5/00bGj2ybnmGNqmJtFpMGogBVJIRUV8MwzbtPzv/5y62puvBFarXcTi3UIhdya1zXXwIa015yIiEhNJSoSFIwuIDI3Qrw8ji3zQ3k+l19RxI03+MjL8zpCkeykAlYkBVjrOhhefTV8/jnssw88+ij07FmHi/l8rmFTOOymDQcC6kIsIiJSS2/NCjPlpwglFZUzmnJjtOgaoc+JYfLytCRHxCtq4iTisUgEDjnENQ7+80949lmYPLmOxetKPp9b73rdde6rilcREZEa+/hjuPi2YkpWVO8pUZKIE10Q9SYoEQFUwIp4Zvp0OPZY2Gsv11X4oYdgxgwYMABy9F+miIhI0n3zDRx9tJsJtXhGkGY5/mrP+3P9BNoFPIlNRBz9b7JIkv3yC5x1FuyyC7z7LgwdCt9/DxddBM2aeR2diIhI9vnpJzjjDDf76b334JZb4Kd3Q+zXJZ+83DwMhrzcPPI75BPqqp4SIl7SGliRJFm0CG67DYYNc2teL7nEdRpu06byhETCrVstLnadhLVuVUREpFEtXOhy83/+47oJ/9//wVVXwRZbAPgoGlBEeHaY6IIogXYBQl1D+HKUm0W8pAJWpJHFYm6z87vvdt8PHOg6C2+7bZWTEgm3CHbNzsFFRSpiRUREGtjSpXD//XDPPS7tnnEGDBkC22xT/Txfjo/CboUUdlPTJpFUoSnEIo2krMyNtnbtCtdfDwcfDF9+CU8/vUbxCm7kNRJxFa617msk4o6LiIhIgygtdT0ntt/eFayHHebWvT755NrFq4ikJhWwIg2sogKeew523BEuvBC6d4ePPoJXXoEePdbzouJidwu4qnjcbYMjIiIitZKoSDD+u/Hc/P7NjP9uPGXlCZ591uXmSy5xfSgiERg71h0TkfShKcQiDcRaN2B69dVupDUQcI8LCty6mg0KBt204Vhs9TG/311EREREaixRkaBgdAGReRHiZXGa5fjxzc8n/lgRvYI+hg+HQw+tQW4WkZSkEViRerLW7dt64IFwxBGuBn3+eZg2Dfr1q2GCDIXcmte8PPeCvDz3OKROhyIiIrURnh0mMjdCrCyGxVJSEWPZ5hGufCLMZ5+5acMqXkXSl0ZgReqoogJefx3uustteN6unetieNZZkJtby4v5fK5hUzjspg0HAupCLCIiUkuJBDw3sZhYWRyqFqlN47TaIUpOjpoxiaQ7FbAitVRaCs8+6zoXzpwJXbrAI4+4DoYtW9bjwj4fFBa6PyIiIlJjJSXwzDMuN88yQcwJfmzT1cty/Ll+Au0C3gUoIg1GU4hFamjxYrjjDujcGc4+2xWrL74I301PcEGn8bS892YYP97d/hUREZFG99dfbh/Xzp3hX/+CTTaBF28OcdAO+eTl5mEw5OXmkd8hn1BXLcsRyQQagRXZiLlz4cEH4fHH3b5xffvC6NFuWxxTof1bRUREku2XX9w+rsOHu/Tbrx9ceSX06QPG+Di+oojw7DDRBVEC7QKEuobw5Sgvi2QCFbAi6/HNN24q0nPPufWuJ54IV1zhGgavUnX/Vqi+f6umAouIiDSor7+Gu+92zRKthZNPhsGDYbfdqp/ny/FR2K2Qwm7KxSKZRlOIRapY2VH4H/9we8SNGQPnnguzZrlkWa14Be3fKiIi0sishfffd53+d93V7d164YXw/feuJ8WaxauIZDaNwIrgRlhfe811FP7kE2jTBm66Cc4/332/Xtq/VUREpMEkKhKEZ4cpnl9Mz62ClH4d4t67fXz6KbRtCzff7HLz5pt7HamIeEUFrGS1kpLVHYW/+851FB42DE4/vYYdhVfu37rmGljt3yoiIlIriYoEBaML+GRuhHh5HFPux/6Sz3Z/FvHooz4GDoQWLbyOUkS8pgJWstKPP8JTT8GTT8Jvv0GvXm668LHHQpPa/Feh/VtFREQaxIgPw3zwQ4Ry42Y12aYxmm0f4b6rwxy1o9ayioijAlayRnk5vP6661g4YQIY42rNyy6r7ChsNn6NddL+rSIiInVSVuaW8DzxBEwoLYY+caiSj8tsnK9+j6qAFZFVVMBKxvv+ezfSOmKEG23t2BGGDIEzBybY5uswfFQMy4MaORUREUmSWbNW5+aFC2HbbeGUs4K84vOzbMXqvhL+XD+BdgHvAhWRlKMCVjJSWRm8+qobbZ04cfUg6TnnuG1bfWj/VhERkWQqLYVXXnG5+b33XLo98kiXmw87DDAhFozOJzIvQrwsjj/XT36HfEJd1VdCRFZTASsZ5bvv3DSkUaPcHd1OnVzHwjPOgA4dqpw4Xvu3ioiIJMOMGatz86JFrmHirbe63Ny+fdUzfRQNKCI8O0x0QZRAuwChriF8ObqxLCKrqYCVtFdSsvqO7qRJrgnTyju6hx66ngHVDe3fqgJWRESkxqpufRNsHyTUNUR5mY+xY11u/uADl5uPPtrl5kMOgZycdV/Ll+OjsFshhd2Ui0Vk3VTAStqaPt3d0X3mGXdHd7vt4Pbb3RY47dpt5MXav1VERKTeVm59s3Labwufny1K8ln6aBGL//Sx/fZw550wcCBstZXX0YpIJlABK2nlr7/caOvIkTB5MjRtCsccA2ef7ToJr++O7lq0f6uIiEi9hWeHK/dtdTeElyViLKuIsP8xYW78ZyF9+tQiN4uI1IAKWEl5f//tWuyPGQPvFiU4dEWYYzcv5pLTg+x/W4gt29dhbYz2bxUREamzJUtcs8RbJxcT7xCHKkWqaRan72lRDj5A04BFpOGpgJWUFIvBG2/ASy+5GrO0FDpvk+CrrQvYbmGEnL/imJf98Es9Ogdr/1YREZEaW7rU5eYxY+Dtt13H/7b7Bmna0U852vpGRJLD80kdxphtjDHvGWO+NcZ8Y4y5xOuYxBvLlsHLL8MJJ8CWW8I//wmffgrnngsffQQ/DAuzw58RfMtjGGurdw4WEZEGo9wsK8Xj7mby8ce73HzKKTBtGlxwAXzyCfz6fogDtssnLzcPgyEvN09b34hIo0qFEdgVwOXW2s+NMa2AacaYCdbab70OTBpfSYm7i/vSS/D66y5RbrklnHkm9O8P++5bZe3MzeocLCKSJMrNWWBd3YN9OT5KSty94TFj3IjrsmWuAdOgQS4377NP1XWt2vpGRJLL8wLWWjsfmF/5/VJjzHSgA6AkmaHKymDCBJcYX3vNrXHdYgt3V7d/fzjwwPXMCFbnYBGRpFBuznxrdg9u2dTP9s3z2aW4iDde87F0KbRpA6edBieeCAccsP7VOtr6RkSSyfMCtipjTGcgCETW8dw5wDkA2267bXIDk3orLXV7tL70kusi/NdfsNlmcNxxrmg9+GDXUXiD1DlYRCTplJsz05rdg+PlMb6MRfj+2zAnnVhI//5w0EFu/1YRkVSSMv8sGWPygLHApdbav9d83lo7HBgO0Lt3b5vk8KQO5sxxU5DCYXjv3QQHLg+zV24x1+4TpPulIfqGfOTm1uKC6hwsIpJUys2Z58cfXRp9MFpMvP3a3YOvuCfKkIM0kioiqSslClhjTFNcgnzOWjvO63ikbkpL4cMP4a23XHKcPt0d365Tgs82L6Drogi+0jhmqh8ezofCIqCWxac6B4uIJIVyc2YoLXX7pq/MzTNmuONb7R+kafu1uwfv3iHgTaAiIjXkeQFrjDHAU8B0a+19XscjtfPzz6tHWd99183szc11a2XOPtsNkHafFcb8MwIllUmyavdgFaIiIilHuTm9/fTT6tw8ceLq3NynD/zrX3D44bDd9iH6PZe/ag2sP9ev7sEikhY8L2CBfYFTga+MMdHKY9dYa9/yLiRZn7IymDLFJcW33oJvvnHHO3WCU091SfGggyAvr8qL/qvuwSIiaUa5OUWtq3NwYoWv2gyobytbbXXu7JowrczNfn/VK6l7sIikJ88LWGvth4DxOg5Zh0QCwmGWTCrmk9IgT8wN8c5E15mwaVPYf3844ww3yrrTTmDW91NU92ARkbSi3Jya1uwcnGv8tFqSz/Iniogv9dG0qevkf9ZZLjfvuOMGcjPqHiwi6cnzAlZSi7Vu6tGH7yfY/ZoCtp0fIc/G2Qc//mb5tDm1iH5H+DjkEGjVqoYXVfdgERGROrPWNUZ8uCjM+/MjrMhxN4RLbYyyFhEKzghz3iGFHHzwGjOgREQykArYLFdR4ZotTZ7s/nzwAcydC0cQ5mgi+CubO7Qixr5NI+x3VB3Wrap7sIiISI1VVLhpwFVz87x5wAHF0GeNJTlN4+x7XJQjD9AoqohkBxWwWaa8HIqLVyfFDz+ERYvcc+3bu2nB++8Px00vxv9oHKpsimDqs25V3YNFRCSLrWvt6sr1pitz8wcfrM7Nf/7pXrf11qtzc0XXINdM8xMrr945ONAu4MHfSETEGypgM1Hl2lWKiyndOcjHm4X4YIqPyZPh449X91Pafns48sjViXH77auslRkfhGe0blVERKS+1ly72rKpn27+fP6xuIgpk318/DEsW+bO3WEHOPro1bl5u+1W5+ZERYjXflPnYBHJbipgM0hFBfwwK0He8QW0nhmhaXmcMvysIJ+hFNGjp48zzlidFNu338DFtG5VRESk3ioq4KkPwkyZE6HEupvC8fIYxb9FKB4bZrcWhZx1lsvL++234dzsy1HnYBERFbBpKpGA776Dzz+HadPc1+Ji2P/vMC8QoVmVtat9mkdYPCJM3km1mL6rdasiIiK1kkjAzJkuJ6/Mz8XFsDRYuXY1Z/W5plmcax6Ockvf2i2tUedgEcl2KmBTRZVpvwSD1YrFFStcM4eqCTEaXT3dqHlzV18OGABnzism7/Xqa1eblMbJmx0F6tB8SetWRUQkS21s3er06atvIn/+efXc3KIF7Lab24fVt1OQJ/7yszxRfe3qXp0Dyf9LiYikORWwqSCRgIICiESw8TgVzf0s6JTPrQcUMS3q48svoaTEner3u/p20CDYfXfo1cvt89Zk5U9yfBAmau2qiIhIfay5brW5z08nXz77/1RE9HOXm0tL3bl5eS43n3326tzcvfvq3JyoCPHNaK1dFRFpCCpgPRKLuWlGM2ZAxethjp8UoUUihgF8y2NsMiPCop/D+PMLueAClwx33x26dt3ILF6tXRUREamTWMzl5Rkz4I2ZYSaZCAmfuyG8PBFjxrIIv3waJr91IRdd5HJzr16u8VJOzvqvq7WrIiINRwVsfW1g6q+18NtvborRjBmrv86YAb/8svoSN5hicm31fd3yTJwXr4pirteeqyIiIrWxoam/1sKCBdVz8srv585dfQ3Tpxh7YPXcbJrF+ff9Ua4/sPZLa7R2VUSkYaiArY/Kqb+2crQz0czP3A75DN23iG9n+pgxA5YsWX16Xp6b7nvgge7rTju5rzvMDOI7rfq0X+P3QzBQt7i0dlVERLLUqqm/cyPEy+M0y/HTweazzw9FzJzucvPff68+f2VuPugg93Vlfp5REeS01/3EyqqvWw22DyT/LyUiIqtkZwG7gVHTdSkthZ9/hp9+gjlz3NeffoL208JcPz1CXmVb/CYlMbb4PsKKP8O0DBRyyinVk2GHDlX2Wa2qu6b9iohIdtvQqOm6lJSszs1V8/PnsTDTe0SwTV1uLqmI8X1phD9nhAnmFTJgwOobyDvtBFtvve7c3K0iRP40rVsVEUk12VfAVmmYtLJYXNE7n+8eKmLOL761EuFPP8H8+dUvkZPjitEbcoppuY6pv8/8XxSu05Y1IiIiNbFmwyR/rp/d2+Xz0B5FzP3ZVy0nr/x+wYLq18jJgY4dwRxYjG2y9tTf/7srynUH1Dw3a92qiEhqyr4CNhx2xevK6bqxGMsnRbiyZ5g3K7eZadoUtt0WOnVydWSnTqv/dO7sitemTXEdf09ex9TfunT81bRfERHJUuHZYSLzIqum68bKYrw/K8Jud4Thu+q5uXNnOPzw1Tl5ZX5emZvHfxfk5LFrT/0NtAvUOi6tWxURST3ZV8AWF7uR1yr8xLnthCjXXFpIp07Qvv2Guwmuoo6/IiIi9VY8v5h4WfXcTLM4J1wU5bLdXW5u165muTnUNUR+B039FRHJVNlXwAaDrtCsMmqak+en52kB2KeW19LUXxERkXoLtg/iz60+apqX6+e0vgH27la7a2nqr4hIZsu+ArahR0019VdERKReGnrUVFN/RUQyV/YVsBo1FRERSSkaNRURkZrKvgIWNGoqIiKSYjRqKiIiNVGTVkUiIiIiIiIinlMBKyIiIiIiImlBBayIiIiIiIikBRWwIiIiIiIikhZUwIqIiIiIiEhaUAErIiIiIiIiaUEFrIiIiIiIiKQFFbAiIiIiIiKSFlTAioiIiIiISFpQASsiIiIiIiJpwVhrvY6h1owxC4GfGuBSbYA/GuA6XkjX2BV3cinu5FLcydWQcXey1rZtoGtlJeXmtI0b0jd2xZ1ciju5FPd6cnNaFrANxRgz1Vrb2+s46iJdY1fcyaW4k0txJ1e6xi0blq4/13SNG9I3dsWdXIo7uRT3+mkKsYiIiIiIiKQFFbAiIiIiIiKSFrK9gB3udQD1kK6xK+7kUtzJpbiTK13jlg1L159rusYN6Ru74k4uxZ1cins9snoNrIiIiIiIiKSPbB+BFRERERERkTSRsQWsMaafMWamMWa2MeaqdTzfzBgzpvL5iDGmc5Xnrq48PtMYU5Bicf+fMeZbY8yXxpiJxphOVZ5LGGOilX9eT7G4TzfGLKwS36Aqzw00xsyq/DMwxeK+v0rM3xljFld5zsvP+2ljzO/GmK/X87wxxjxU+ff60hjTq8pzXn7eG4v7lMp4vzLGfGSM2a3Kc3Mqj0eNMVOTF3WN4u5jjFlS5ffhhirPbfB3rDHVIO4rqsT8deXv9OaVz3n5eW9jjHmv8t+6b4wxl6zjnJT8HZcNU25Wbm6guJWbG5Bys3JzDeNOndxsrc24P4AP+B7YDsgFvgB2XuOc84HHKr8/CRhT+f3Olec3A7pUXseXQnEfBLSs/P68lXFXPo6l8Od9OvDIOl67OfBD5dfWld+3TpW41zj/IuBprz/vyvc+AOgFfL2e5w8HwoAB9gIiXn/eNYx7n5XxAKGVcVc+ngO0SdHPuw8wvr6/Y8mOe41z/wH8L0U+7/ZAr8rvWwHfrePflJT8HdefDf5clZtT7/M+HeXmhoxduTm1Pu8+KDc3ZNwpk5szdQR2T2C2tfYHa20Z8CJw1BrnHAWMqvz+ZeAQY4ypPP6itbbUWvsjMLvyeikRt7X2PWvtssqHnwAdkxTbhtTk816fAmCCtfZPa+1fwASgXyPFuabaxn0y8EJSItsIa+0HwJ8bOOUo4BnrfAJsZoxpj7ef90bjttZ+VBkXpM7vd00+7/Wpz38b9VbLuFPp93u+tfbzyu+XAtOBDmuclpK/47JBys3JpdycZMrNyaXcnFyplJsztYDtAPxS5fFc1v6AV51jrV0BLAG2qOFrG0tt3/ss3F2OlZobY6YaYz4xxhzdCPGtT03jPq5yOsHLxphtavnaxlDj966cDtYF+F+Vw1593jWxvr+bl593ba35+22Bd4wx04wx53gU04bsbYz5whgTNsb0qDyWFp+3MaYlLpGMrXI4JT5v46aQBoHIGk9lwu94tlFuVm6uCeXm1P53S7k5SZSb169JXV8o3jLGDAB6AwdWOdzJWjvPGLMd8D9jzFfW2u+9iXAtbwAvWGtLjTH/wt1hP9jjmGrjJOBla22iyrFU/rzTmjHmIFyS3K/K4f0qP+8tgQnGmBmVdzFTwee434eYMeZw4FVgB29DqpV/AFOstVXvCHv+eRtj8nCJ+1Jr7d/JfG+RulBuTjrl5iRSbk465eb1yNQR2HnANlUed6w8ts5zjDFNgE2BRTV8bWOp0XsbYw4FrgWOtNaWrjxurZ1X+fUHYBLuzkgybDRua+2iKrE+Cexe09c2otq890msMYXDw8+7Jtb3d/Py864RY0xP3O/IUdbaRSuPV/m8fwdeIXnTBzfKWvu3tTZW+f1bQFNjTBvS4POutKHfb08+b2NMU1yCfM5aO24dp6Tt73gWU25Wbq4J5eYU/HdLudkTys3rYz1YBNzYf3Ajyz/gppWsXJzdY41zLqB6o4iXKr/vQfVGET+QvEYRNYk7iFt4vsMax1sDzSq/bwPMIkkL0msYd/sq3x8DfFL5/ebAj5Xxt678fvNUibvyvB1xi+ZNKnzeVWLozPobFxxB9UX0n3r9edcw7m1xa9v2WeO4H2hV5fuPgH4pFHe7lb8fuGTyc+VnX6PfMa/irnx+U9xaHH+qfN6Vn90zwAMbOCdlf8f1Z70/M+Vm5eYGibvyPOXm5MWt3JzEuCufV27eUCzJ/GEl+RfjcFx3rO+BayuPDcXdGQVoDvy38j/IT4Htqrz22srXzQRCKRb3u8BvQLTyz+uVx/cBvqr8j/Ar4KwUi/t24JvK+N4Ddqzy2jMrfw6zgTNSKe7KxzcCd6zxOq8/7xeA+UA5bh3BWcC5wLmVzxtgWOXf6yugd4p83huL+0ngryq/31Mrj29X+Vl/Ufl7dG2KxX1hld/vT6iS5Nf1O5YqcVeeczquOU7V13n9ee+HW+fzZZXfhcPT4Xdcfzb6s1VuTq24lZsbNm7l5tSKW7m5YeNOmdy88q6EiIiIiIiISErL1DWwIiIiIiIikmFUwIqIiIiIiEhaUAErIiIiIiIiaUEFrIiIiIiIiKQFFbAiIiIiIiKSFlTAioiIiIiISFpQASsiIiIiIiJpQQWsSIYyxnQ0xvT3Og4RERFxlJtF6k8FrEjmOgTo5XUQIiIisopys0g9GWut1zGISAMzxuwHvAYsBpYCx1prf/A0KBERkSym3CzSMFTAimQoY8zbwGBr7ddexyIiIiLKzSINQVOIRTJXd2CG10GIiIjIKsrNIvWkAlYkAxlj2gBLrLUrvI5FRERElJtFGooKWJHM1Bn41esgREREZJXOKDeL1JsKWJHMNANoY4z52hizj9fBiIiIiHKzSENQEycRERERERFJCxqBFRERERERkbSgAlZERERERETSggpYERERERERSQsqYEVERERERCQtqIAVERERERGRtKACVkRERERERNKCClgRERERERFJCypgRUREREREJC38P0dNtyU2JcYhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# vamos ainda colocar um gráfico do lado do outro para facilitar (note os comandos para fazer isso)\n",
    "fig, ax = plt.subplots(1,2,figsize=(16,6)) # note que ax[0] sera 1 gráfico e ax[1] o outro e você pode controlar o tamanho do gráfico\n",
    "\n",
    "# Euler\n",
    "\n",
    "ax[0].plot(t, yexato, 'b', label='Sol. Exata')\n",
    "ax[0].plot(t, y_euler, 'ro', markersize=5, label='Método de Euler')\n",
    "ax[0].legend()\n",
    "\n",
    "\n",
    "# odeint\n",
    "ax[1].plot(t, yexato, 'b', label='Sol. Exata')\n",
    "ax[1].plot(t, y_solver, 'go', markersize=5, label='Solver odeint')\n",
    "ax[1].legend()\n",
    "\n",
    "for myax in ax: # como temos 2 eixos, faz um laço para escrever o label em cada um deles\n",
    "    myax.set(xlabel=f'$t$', ylabel=f'$y(t)$')\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sL2Djs2Hieqy"
   },
   "source": [
    "Observando o gráfico, é possível ver que utilizando o mesmo passo (h=0.1) na função `odeint`, o resultado é muito mais preciso do que utilizando o método de Euler. Você pode ainda explorar outros solvers no [SciPy](https://docs.scipy.org/doc/scipy/reference/integrate.html). Mas agora, uma coisa importante que você precisa saber:<br>\n",
    "Nas soluções numéricas de equações diferenciais, o que é a **rigidez** ?<br>\n",
    "<br>\n",
    "A rigidez é caracterizada pela existência de duas ou mais escalas distintas para a variável independente e os problemas de valor inicial com equações rígidas são relativamente comuns, particularmente no estudo de ecossistemas e reações biogeoquímicas. Neste caso, outros métodos numéricos podem ser explorados. Você pode pesquisar isso e ver os tutoriais no próprio site do [SciPy.org](https://docs.scipy.org/doc/scipy/reference/integrate.html).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xX5yyvyHaJCp"
   },
   "source": [
    "--**fim** --"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "[EXTRA 02] Exercício_EDO_Solução_Numérica_Euler_e_ODEINT.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}