introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to MAT230\n",
    "The first Jupyter file. \n",
    "\n",
    "You run the individual cells by pushing at the right place or by holding *shift* and hit *enter* on your keyboard.\n",
    "\n",
    "If you're running this Jupyter file through *mybinder.org* you can play around with the code without altering the file on the server. Next time you open it will be back to the one on the net. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Read libraries\n",
    "\n",
    "First, let's read some libraries. Enter the cell below, hold *shift* and hit *enter*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a comment line. \n",
    "from sympy import dsolve, Eq, symbols, Function, lambdify,solve, init_printing, pprint, latex # This is also a comment. \n",
    "from sympy.plotting import plot\n",
    "from scipy.integrate import solve_ivp\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The first problem demonstrating *sympy* and *numpy*\n",
    "Let's solve this problem:\n",
    "$$\\frac{d x(t)}{d t} = 1-x(t)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First solve it symbolically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking...\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} x{\\left(t \\right)} = 1 - x{\\left(t \\right)}$"
      ],
      "text/plain": [
       "Eq(Derivative(x(t), t), 1 - x(t))"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Declare sympy symbols\n",
    "t = symbols('t')\n",
    "x = symbols('x', cls=Function)\n",
    "\n",
    "#define the problem\n",
    "deqn1 = Eq(x(t).diff(t), 1- x(t))\n",
    "\n",
    "print(\"Checking...\")\n",
    "\n",
    "deqn1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the solution is with initial conditions is:\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x{\\left(t \\right)} = C_{1} e^{- t} + 1$"
      ],
      "text/plain": [
       "Eq(x(t), C1*exp(-t) + 1)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solve the problem with sympy using initial conditions:\n",
    "\n",
    "sol1 = dsolve(deqn1, x(t))#,ics=ics)\n",
    "\n",
    "#Print the solution\n",
    "print('the solution is with initial conditions is:')\n",
    "sol1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the general solution. \n",
    "\n",
    "Lets use the initial conditions, $x(0)=x_0$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C =  [X0 - 1]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "C1=symbols('C1')\n",
    "X0=symbols('X0')\n",
    "C=symbols('C')\n",
    "sol1.subs(t,0).subs(x(0),X0)\n",
    "C=solve(sol1.subs(t,0).subs(x(0),X0),C1)\n",
    "print('C = ', C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x{\\left(t \\right)} = C_{1} e^{- t} + 1$"
      ],
      "text/plain": [
       "Eq(x(t), C1*exp(-t) + 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol1.subs(C1,C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the solution is with initial conditions {x(0): 2}  is:\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x{\\left(t \\right)} = 1 + e^{- t}$"
      ],
      "text/plain": [
       "Eq(x(t), 1 + exp(-t))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solve the problem with sympy using initial conditions:\n",
    "ics={x(0):2} \n",
    "\n",
    "sol2 = dsolve(deqn1, x(t),ics=ics)\n",
    "\n",
    "#Print the solution\n",
    "print('the solution is with initial conditions', ics, ' is:')\n",
    "sol2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Lambdify the equation\n",
    "g=lambdify(t,sol2.rhs,'numpy')\n",
    "\n",
    "# Now we can evaluate at numerical $t$s\n",
    "g(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x281231490>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6z0lEQVR4nO3deXxV1b338e8ZkpP5hIQkJCEQ5kEkRJBBQFQUBIq1tWpFC0rtFYtSy9VWtI9tH28vtz7VVsXh2kKpFXFgqrPgQCIKajBhnhMgkISQACfzePbzRyAVSSAnJNnnnHzer9d+Jeyszfmt12rJ173XXstiGIYhAAAAk1jNLgAAAHRuhBEAAGAqwggAADAVYQQAAJiKMAIAAExFGAEAAKYijAAAAFMRRgAAgKnsZhfQEm63W3l5eQoPD5fFYjG7HAAA0AKGYai0tFQJCQmyWpu//+ETYSQvL09JSUlmlwEAAFohNzdX3bt3b/bnPhFGwsPDJTV0JiIiwuRqAABAS5SUlCgpKanx93hzfCKMnHk0ExERQRgBAMDHXGiKBRNYAQCAqQgjAADAVIQRAABgKsIIAAAwFWEEAACYijACAABMRRgBAACmIowAAABTEUYAAICpPAojCxcu1OWXX67w8HDFxsbqxhtv1J49ey54XVpamoYPH66goCD17t1bL774YqsLBgAA/sWjMJKWlqa5c+dq06ZNWrdunerq6jRp0iSVl5c3e01OTo6mTp2q8ePHKzMzU4888ojmzZunlStXXnTxAADA91kMwzBae/Hx48cVGxurtLQ0XXnllU22+fWvf6233npLu3btajw3Z84cbdmyRRs3bmzR55SUlMjpdMrlcrE3DQAAPqKlv78vas6Iy+WSJEVFRTXbZuPGjZo0adJZ5yZPnqyMjAzV1tY2eU11dbVKSkrOOtrDRzuPae6r32h/YVm7/P0AAODCWh1GDMPQ/PnzNW7cOA0ZMqTZdgUFBYqLizvrXFxcnOrq6lRUVNTkNQsXLpTT6Ww8kpKSWlvmeb361WG9uzVfb2/Ja5e/HwAAXFirw8h9992nrVu3avny5Rds+92tg888GWpuS+EFCxbI5XI1Hrm5ua0t87ymp8RLkt7ZmqeLeFoFAAAugr01F91///166623lJ6eru7du5+3bbdu3VRQUHDWucLCQtntdkVHRzd5jcPhkMPhaE1pHrl2UJwcdqsOHC/XrvxSDU5gPgoAAB3NozsjhmHovvvu06pVq/TJJ5+oV69eF7xmzJgxWrdu3Vnn1q5dqxEjRiggIMCzattYeFCArh4QK0l6eyuPagAAMINHYWTu3Ll65ZVX9Oqrryo8PFwFBQUqKChQZWVlY5sFCxZo5syZjX+eM2eODh06pPnz52vXrl1asmSJFi9erAcffLDtenERpqckSJLe3sKjGgAAzOBRGHnhhRfkcrl01VVXKT4+vvF4/fXXG9vk5+fr8OHDjX/u1auX3nvvPa1fv17Dhg3T448/rmeeeUY33XRT2/XiIlwzMFYhgTYdOVmprNxTZpcDAECn49GckZbcOVi6dOk55yZMmKBvvvnGk4/qMMGBNl03OE7/ysrTO1vzldqji9klAQDQqbA3jaTpQxse1byzNU9uN49qAADoSIQRSeP7d1V4kF3HSqr19cETZpcDAECnQhiR5LDbdP0l3STxVg0AAB2NMHLambdq3t9WoLp6t8nVAADQeRBGTruiT7SiQgNVXF6jjdnFZpcDAECnQRg5zW6zauqlpx/VsFcNAAAdhjDyLd87/VbNB9sLVF1Xb3I1AAB0DoSRb7k8OUpxEQ6VVNXps71N7ygMAADaFmHkW2xWi6Zdenp5eN6qAQCgQxBGvmN6Srwk6aOdx1RZw6MaAADaG2HkO4YlRap7l2CV19Tr0z2FZpcDAIDfI4x8h8ViaZzIyls1AAC0P8JIE848qvlkd6FKq2pNrgYAAP9GGGnC4PgI9Y4JVXWdWx/tOmZ2OQAA+DXCSBMsFsu/d/Ldkm9yNQAA+DfCSDPOPKpJ33dcpypqTK4GAAD/RRhpRt/YcA3sFq7aekMf7igwuxwAAPwWYeQ8zuzk+zaPagAAaDeEkfM4M2/kiwNFOl5abXI1AAD4J8LIefSIDlFKUqTchvT+du6OAADQHggjFzB9aMNE1tWZR02uBAAA/0QYuYAbhiXIZrUo8/ApHTheZnY5AAD4HcLIBcSGB+nKfl0lSau+OWJyNQAA+B/CSAv8aHiSJGnVN0dV7zZMrgYAAP9CGGmBiYNiFRFkV76rShsPFJtdDgAAfoUw0gJBAbbGNUdW8qgGAIA2RRhpoR8N7y6p4RVfdvIFAKDtEEZaaFhSpHrHhKqq1q33t7E8PAAAbYUw0kIWi0U3XdZwd2QFj2oAAGgzhBEP/PCyRFks0lc5J3S4uMLscgAA8AuEEQ/EO4M1rm/DmiNMZAUAoG0QRjx05lHNqswjcrPmCAAAF40w4qHJl3RTmMOu3BOV+vrgCbPLAQDA5xFGPBQcaNPUS7tJ4lENAABtgTDSCmeWh393a74qaupMrgYAAN9GGGmFy5O7qEdUiMpr6vXhDtYcAQDgYhBGWsFiseiHlyVKklZuPmpyNQAA+DbCSCudeavm8wNFyjtVaXI1AAD4LsJIKyVFhWhUrygZhrQ6k7sjAAC0FmHkItx0evO8lZuPyDBYcwQAgNYgjFyEqZfGKzjApuyicmXmnjK7HAAAfBJh5CKEOeyaMqRhzZEVm1lzBACA1iCMXKQzj2re2ZKnqtp6k6sBAMD3EEYu0pje0UpwBqmkqk4f7TpmdjkAAPgcwshFslot+sHpNUfezOBRDQAAniKMtIEzy8On7zuuIycrTK4GAADfQhhpA726huqKPtEyDOn1r3PNLgcAAJ9CGGkjM0b1kNQQRmrr3SZXAwCA7yCMtJFJg7upa1igCkur9TETWQEAaDHCSBsJtFt184iGuSPLvjxscjUAAPgOwkgbuu3yhkc1n+0r0qHicpOrAQDANxBG2lCP6BBd2T9GkrT8KyayAgDQEoSRNjZjZMPdkTczclVTx0RWAAAuhDDSxiYOilVchEPF5TX6cEeB2eUAAOD1CCNtLMBm1a2nJ7K+ykRWAAAuiDDSDm4d2UNWi7Qxu1gHjpeZXQ4AAF6NMNIOEiODdfWAWEnScu6OAABwXoSRdnJmRdYV3xxRVW29ydUAAOC9PA4j6enpmj59uhISEmSxWLRmzZoLXrNs2TKlpKQoJCRE8fHxuuuuu1RcXNyaen3GVQNileAM0qmKWr2/Pd/scgAA8Foeh5Hy8nKlpKRo0aJFLWq/YcMGzZw5Uz/96U+1Y8cOvfnmm/r666919913e1ysL7FZLfrx6dd8mcgKAEDz7J5eMGXKFE2ZMqXF7Tdt2qTk5GTNmzdPktSrVy/dc889euKJJzz9aJ9z6+VJevrjffr64EntPVaq/nHhZpcEAIDXafc5I1dccYWOHDmi9957T4Zh6NixY1qxYoWmTZvW7DXV1dUqKSk56/BFcRFBunZQw0RW7o4AANC0Dgkjy5Yt06233qrAwEB169ZNkZGRevbZZ5u9ZuHChXI6nY1HUlJSe5fZbmaM6ilJWvnNEVXWMJEVAIDvavcwsnPnTs2bN0+PPfaYNm/erA8++EA5OTmaM2dOs9csWLBALper8cjN9d19Xsb37aqkqGCVVtXp7a15ZpcDAIDXafcwsnDhQo0dO1YPPfSQhg4dqsmTJ+v555/XkiVLlJ/f9FsmDodDERERZx2+ymq16DYmsgIA0Kx2DyMVFRWyWs/+GJvNJkkyDKO9P94r3Dw8SXarRVm5p7Qjz2V2OQAAeBWPw0hZWZmysrKUlZUlScrJyVFWVpYOH274r/4FCxZo5syZje2nT5+uVatW6YUXXlB2drY+//xzzZs3TyNHjlRCQkLb9MLLxYQ7NHlIN0ncHQEA4Ls8DiMZGRlKTU1VamqqJGn+/PlKTU3VY489JknKz89vDCaSdOedd+qpp57SokWLNGTIEN18880aMGCAVq1a1UZd8A23n35UsybzqEqrak2uBgAA72ExfOBZSUlJiZxOp1wul8/OHzEMQ9c+laYDx8v12PcGa/a4XmaXBABAu2rp72/2pukgFotFd41tCCBLvzioerfXZ0AAADoEYaQD/fCyRDmDA3T4RIU+3nXM7HIAAPAKhJEOFBJob9zNd8nnOSZXAwCAdyCMdLCZY3rKZrVoU/YJXvMFAECEkQ4X7wzW1EvjJUlLNhw0txgAALwAYcQEs8cmS5Le3pKnwtIqc4sBAMBkhBETpPbootQekaqpd2vZJhZBAwB0boQRk8w+/Zrvsi8PqaqW3XwBAJ0XYcQkU4Z0U4IzSEVlNXp7C7v5AgA6L8KISew2q2ZekSxJWvL5wU6zaSAAAN9FGDHRjy9PUnCATbvyS7Qp+4TZ5QAAYArCiIkiQwJ10/BESSyCBgDovAgjJrvzioaJrB/tOqaDReUmVwMAQMcjjJisb2yYrhoQI8No2EAPAIDOhjDiBX46ruHuyJsZuSqpqjW5GgAAOhZhxAuM69tV/WLDVF5Trze+zjW7HAAAOhRhxAtYLBbNPn13ZOkXB1Xv5jVfAEDnQRjxEj9ITVSXkAAdOVmpdTuPmV0OAAAdhjDiJYICbJoxqockXvMFAHQuhBEv8pPRybJbLfoq54S2HjlldjkAAHQIwogX6eYM0g0pCZKkF9YfMLkaAAA6BmHEy8y5qo8k6YMdBdpfWGpyNQAAtD/CiJfpHxeuSYPjZBjSC+uzzS4HAIB2RxjxQj+/uq8k6V9ZR3XkZIXJ1QAA0L4II15oWFKkxvXtqjq3ob+mc3cEAODfCCNe6uen54689nWujpdWm1wNAADthzDipcb0idawpEhV17n1d9YdAQD4McKIl7JYLJp7eu7IPzcekquSDfQAAP6JMOLFJg6MVf+4MJVW1+mVTYfMLgcAgHZBGPFiVqtFP7+q4e7I4g05qqypN7kiAADaHmHEy31vaLySooJ1orxGr3992OxyAABoc4QRL2e3WXXPlQ1v1ryUnq2aOrfJFQEA0LYIIz7gR8O7KybcoTxXldZkHTW7HAAA2hRhxAcEBdj0s/G9JEkvph1QvdswuSIAANoOYcRHzBjVU87gAGUfL9eHOwrMLgcAgDZDGPERYQ67Zl2RLEl67tP9MgzujgAA/ANhxIfcdUWyQgJt2pFXorS9x80uBwCANkEY8SFdQgN128gekqTn1x8wuRoAANoGYcTH/Gx8bwXYLPoq54S+PnjC7HIAALhohBEf080ZpB8N7y5J+vO6vSZXAwDAxSOM+KC5V/dVgM2iLw4U64v9RWaXAwDARSGM+KDuXUIa5448uW4vb9YAAHwaYcRH3Xd1XznsVm0+dFLrebMGAODDCCM+KjYiqHHdkSfX7uHuCADAZxFGfNg9V/ZWaKBN24+W6MMdx8wuBwCAViGM+LDoMIdmj2vYs+apdXvYswYA4JMIIz7u7vG9FRFk195jZXpna57Z5QAA4DHCiI9zBgfoP67sLUn6y0f7VFfvNrkiAAA8QxjxA3eO7aWo0EDlFJVrVeZRs8sBAMAjhBE/EOaw694JfSRJT3+0TzV13B0BAPgOwoifuGN0T8WGO3T0VKVez8g1uxwAAFqMMOInggNtuu+avpKkRZ/sU1VtvckVAQDQMoQRP3Lr5UlKjAzWsZJqvbLpkNnlAADQIoQRP+Kw2zRvYsPdkefXH1B5dZ3JFQEAcGGEET/zw8u6Kzk6RCfKa7T0i4NmlwMAwAURRvxMgM2qB67tL0n637QDclXWmlwRAADnRxjxQ9NTEtQvNkwlVXX6a3q22eUAAHBeHoeR9PR0TZ8+XQkJCbJYLFqzZs0Fr6murtajjz6qnj17yuFwqE+fPlqyZElr6kUL2KwW/eekhrsjf9uQrXxXpckVAQDQPI/DSHl5uVJSUrRo0aIWX3PLLbfo448/1uLFi7Vnzx4tX75cAwcO9PSj4YHJl3TT5cldVFXr1p8+3Gt2OQAANMvu6QVTpkzRlClTWtz+gw8+UFpamrKzsxUVFSVJSk5O9vRj4SGLxaJHpw3Wjc99rpXfHNFdY5M1JNFpdlkAAJyj3eeMvPXWWxoxYoSeeOIJJSYmqn///nrwwQdVWdn8o4Pq6mqVlJScdcBzw5Ii9f1hCZKk/3p3pwzDMLkiAADO1e5hJDs7Wxs2bND27du1evVq/eUvf9GKFSs0d+7cZq9ZuHChnE5n45GUlNTeZfqthyYPUKDdqk3ZJ/TRrkKzywEA4BztHkbcbrcsFouWLVumkSNHaurUqXrqqae0dOnSZu+OLFiwQC6Xq/HIzWWvldbq3iVEd4/rJUla+N4u1daziR4AwLu0exiJj49XYmKinM5/z1cYNGiQDMPQkSNHmrzG4XAoIiLirAOtd+9VfdQ1LFDZReVaxjLxAAAv0+5hZOzYscrLy1NZWVnjub1798pqtap79+7t/fGQFB4U0LgQ2tMf72MhNACAV/E4jJSVlSkrK0tZWVmSpJycHGVlZenw4cOSGh6xzJw5s7H9jBkzFB0drbvuuks7d+5Uenq6HnroIc2ePVvBwcFt0wtc0I8vT1K/2DCdrKjVc5/uN7scAAAaeRxGMjIylJqaqtTUVEnS/PnzlZqaqscee0ySlJ+f3xhMJCksLEzr1q3TqVOnNGLECN1+++2aPn26nnnmmTbqAlrCbrPqkWmDJElLPz+ow8UVJlcEAEADi+ED73uWlJTI6XTK5XIxf+QiGIahmUu+0mf7ijTt0ng9d/tlZpcEAPBjLf39zd40nUjDQmiDZLVI727L1+ZDJ8wuCQAAwkhnM7BbhG4Z0bBuy+Pv7GIhNACA6QgjndD8Sf0VEmhTVu4pvb013+xyAACdHGGkE4oND9K9E/pIkv74/m5V1dabXBEAoDMjjHRSd4/vrW4RQTp6qlJLvzhodjkAgE6MMNJJBQfa9KvrB0iSnv14nwpcVSZXBADorAgjndiNwxJ1WY9IldfU67/e3Wl2OQCAToow0olZrRY9fuMQWS3SO1vztWFfkdklAQA6IcJIJ3dJglMzxyRLkh7713ZV1zGZFQDQsQgj0PxJ/RUT7lB2Ubn+9lmO2eUAADoZwggUERSg35zet+bZT/Yp9wT71gAAOg5hBJKkG1ISNKZ3tKpq3fr920xmBQB0HMIIJDXsW/P4jZfIbrXoo13H9NHOY2aXBADoJAgjaNQ3Nlx3j+8tSfrd2ztUWcNkVgBA+yOM4CzzJvZVgjNIR05W6rlP95tdDgCgEyCM4CwhgXY9Nv0SSdJL6dk6cLzM5IoAAP6OMIJzTL4kTlcNiFFNvVu//dcOGYZhdkkAAD9GGME5LBaLfn/DJQq0W7Vhf5He3ZZvdkkAAD9GGEGTekaH6udX9ZEkPf7OTpVV15lcEQDAXxFG0Kw5E/qoZ3SIjpVU6y/r9ppdDgDATxFG0KygAJt+d0PDZNa/f3FQW4+cMrcgAIBfIozgvK4eEKvvDY1XvdvQr1ZsVU2d2+ySAAB+hjCCC/r9DZcoKjRQuwtKWXsEANDmCCO4oOgwh35/+nHNc5/u1868EpMrAgD4E8IIWuR7Q+M1+ZI41bkNPbRii2rreVwDAGgbhBG0SMNGekPkDA7QjrwSvZSebXZJAAA/QRhBi8WGB+m30wdLkp7+aJ/2His1uSIAgD8gjMAjP0hN1DUDY1VT79ZDK7aqjsc1AICLRBiBRywWi/77B5cqPMiuLbmntOTzHLNLAgD4OMIIPNbNGaT/M63hcc2Ta/cqm519AQAXgTCCVrl5RHeN79dV1XVu/WrFVtW72dkXANA6hBG0isVi0f/cNFShgTZlHDqplzceNLskAICPIoyg1RIjg7Vg6iBJ0hMf7NGh4nKTKwIA+CLCCC7KjJE9NKZ3tCpr6/XrlVvl5nENAMBDhBFcFKvVoj/eNFTBATZtyj6hv39x0OySAAA+hjCCi9YjOkSPTGt4XPPH93ezdw0AwCOEEbSJO0b10LWD4lRT79YvXstUVW292SUBAHwEYQRtwmKx6I83XaqYcIf2FZbpD+/uMrskAICPIIygzUSHOfTULSmSpH9uOqSPdh4zuSIAgC8gjKBNje8Xo7vH9ZIk/WrlVhWWVJlcEQDA2xFG0OYeun6ABsdH6ER5jf7zzS287gsAOC/CCNqcw27TM7cNU1CAVZ/tK2IzPQDAeRFG0C76xobrN6c303vigz3akecyuSIAgLcijKDd3D6qh64bfOZ13yxV1vC6LwDgXIQRtJuG132HKjbcof2FZfrDezvNLgkA4IUII2hXUaGBevL0676vbDqstTsKTK4IAOBtCCNod+P7xehn4xte9/31yq06xuu+AIBvIYygQzw4ueF135MVtbr/1UzV1rvNLgkA4CUII+gQDrtNz91+mcIddn118IT++P5us0sCAHgJwgg6TK+uofp/NzfMH/nbhhy9ty3f5IoAAN6AMIIOdf2QbrpnQm9J0kNvbtH+wjKTKwIAmI0wgg730KQBGt07SuU19br3lc0qr64zuyQAgIkII+hwdptVz952meIiHNpXWKaHV22TYbB/DQB0VoQRmCIm3KHnZlwmu9Wit7fkaekXB80uCQBgEsIITDMiOUqPTB0kSfrDu7u0+dAJkysCAJiBMAJT3TU2Wd8bGq86t6GfL/tGx0urzS4JANDBCCMw1Zn9a/rGhulYSbXmLc9UHQuiAUCnQhiB6UIddr14x3CFBtq0MbtYf1q71+ySAAAdyOMwkp6erunTpyshIUEWi0Vr1qxp8bWff/657Ha7hg0b5unHws/1jQ1rXBDtxbQD+mA7G+oBQGfhcRgpLy9XSkqKFi1a5NF1LpdLM2fO1MSJEz39SHQSUy+N193jGjbUm/9GlnbkuUyuCADQEeyeXjBlyhRNmTLF4w+65557NGPGDNlsNo/upqBzeXjKQO05VqrP9hXp7n9k6F9zxyo2IsjssgAA7ahD5oz8/e9/14EDB/Tb3/62Re2rq6tVUlJy1oHOwW6zatGMy9QnJlT5rir97OUMVdXWm10WAKAdtXsY2bdvnx5++GEtW7ZMdnvLbsQsXLhQTqez8UhKSmrnKuFNnMEBWnLn5eoSEqAtR1z6zze2yO1mhVYA8FftGkbq6+s1Y8YM/f73v1f//v1bfN2CBQvkcrkaj9zc3HasEt6oZ3SoXrxjuAJsFr27LV9/+Yg3bADAX3k8Z8QTpaWlysjIUGZmpu677z5JktvtlmEYstvtWrt2ra655ppzrnM4HHI4HO1ZGnzAqN7R+u8fXKqHVmzVM5/sV++YMN2Ymmh2WQCANtauYSQiIkLbtm0769zzzz+vTz75RCtWrFCvXr3a8+PhB24ekaT9x8v0v2nZ+tXKrUqKCtHwnl3MLgsA0IY8DiNlZWXav39/459zcnKUlZWlqKgo9ejRQwsWLNDRo0f18ssvy2q1asiQIWddHxsbq6CgoHPOA8359eSByjlerrU7j+mef2Zo9c/HKikqxOyyAABtxOM5IxkZGUpNTVVqaqokaf78+UpNTdVjjz0mScrPz9fhw4fbtkp0alarRX++dZgGx0eoqKxGd/8jQ6VVtWaXBQBoIxbDMLz+NYWSkhI5nU65XC5FRESYXQ5Mku+q1PcXfa7C0mpdPSBGf5t1uWxWi9llAQCa0dLf3+xNA58R7wzWX2eOkMNu1ad7juvxd3bKB7I0AOACCCPwKSlJkXrqlmGSpKVfHNQLaQfMLQgAcNEII/A504bG6zfTBkmSnvhgj17/mjlKAODLCCPwSXeP7617r+ojSVqwaps+3MEuvwDgqwgj8Fm/mjxAt4zoLrch3b88U19mF5tdEgCgFQgj8FkWi0X//YNLdd3gONXUuXX3PzK0M49NFQHA1xBG4NPsNquevS1VI5OjVFpdp5lLvtLh4gqzywIAeIAwAp8XFGDTX2eN0MBu4Soqq9ZPlnypwtIqs8sCALQQYQR+wRkcoJdnj1RSVLAOFVfoziVfq4RVWgHAJxBG4DdiI4L0z9mj1DUsUDvzS/Szf2Soqrbe7LIAABdAGIFfSe4aqqV3jVSYw64vc05o3vJM1da7zS4LAHAehBH4nSGJTv115ggF2qxau/OYHngtS3UEEgDwWoQR+KUxfaL1wh2XKcBm0bvb8vXLN7YQSADASxFG4LcmDorTC7cPV4DNore35Gk+gQQAvBJhBH7t2sFxem7GZbJbLXprS54efHOL6t3s9AsA3oQwAr836ZJuWnQ6kKzJytNDBBIA8CqEEXQK1w/ppmdvS5XNatGqzKP61YqtBBIA8BKEEXQaUy6N1zM/bggkK785oodXbpWbQAIApiOMoFOZNjReT/94mGxWi97cfEQLVm0jkACAyQgj6HS+NzRBf751mKwW6fWMXD26hkACAGYijKBTuiHl34Fk+Ve5+s83t7BSKwCYhDCCTuv7wxL151sbHtmszjyqe1/ZzF42AGACwgg6te8PS9RLPxkuh92qj3YVataSr1TKbr8A0KEII+j0Jg6K08uzRyr89OZ6M/76pYrLqs0uCwA6DcIIIGlU72gt/4/Rig4N1LajLt38vxuVd6rS7LIAoFMgjACnDUl06o05Y5TgDFL28XL96IUvdOB4mdllAYDfI4wA39InJkxv3nuFencNVZ6rSre8uFHbj7rMLgsA/BphBPiOxMhgvTFnjC5JiFBxeY1ue2mTvswuNrssAPBbhBGgCV3DHFr+H6M1sleUSqvrNHPJV/pwR4HZZQGAXyKMAM2ICArQy7NHauLAWFXXuTXnlc3622fZMgxWawWAtkQYAc4jKMCm//3JcM0Y1UOGIf3Xu7v0mzXbVcdqrQDQZggjwAXYbVb94cYhenTqIFks0rIvD2v2PzJYHA0A2ghhBGgBi8Win13ZWy/eMVxBAVal7z2um1/cqKOsRQIAF40wAnhg8iXd9MY9YxQT7tDuglLd+Nzn2nrklNllAYBPI4wAHhraPVJr5o7VwG7hOl5arVv+dyNv2gDARSCMAK2QGBmsN+eM0YT+MaqqbXjT5q/pvGkDAK1BGAFaKTwoQItnjdAdoxvetPnDe7v08MptqqqtN7s0APAphBHgIthtVj3+/SH6zbSGN21ez8jVrS9tYpM9APAAYQS4SBaLRXeP762ld41UZEiAtuSe0vRnN2jjAZaQB4CWIIwAbWRC/xi9fd84DY5v2NPmjsVfsmIrALQAYQRoQ0lRIVp57xX6YWqi6t2G/uvdXbp/eaYqaurMLg0AvBZhBGhjwYE2PXlLin5/wyWyWy16Z2u+fvDcFzpYVG52aQDglQgjQDuwWCyadUWylv/HaMWEO7TnWKmmL9qgT3YfM7s0APA6hBGgHV2eHKV37h+n4T27qLSqTrOXZuiptXvYaA8AvoUwArSzuIggLf/ZaP1kdE9J0jOf7NeMv37J678AcBphBOgAgXarHr9xiJ7+8TCFOez66uAJTXn6M32wnWXkAYAwAnSg7w9L1Lvzximlu1OuylrNeWWzfrOGVVsBdG6EEaCD9YwO1ZtzrtA9E3pLkl7ZdFg3LNqgvcdKTa4MAMxBGAFMEGi3asGUQXp59kh1DXNo77EyTX92g5Z9eYhF0gB0OoQRwERX9o/RBw+M14T+Maquc+vR1dv182XfyFVRa3ZpANBhCCOAybqGOfT3Oy/Xb6YNUoDNove3F+j6p9OVvve42aUBQIcgjABewGpt2Gxv1b1j1atrqPJdVZq55CstWLVNZdUsJQ/AvxFGAC9yaXen3p03TndekSxJWv7VYU3+c7o+319kbmEA0I4II4CXCQm063c3XKLlPxutpKhgHT1Vqdv/9qV+s2abyrlLAsAPEUYALzWmT7Q++MWVjSu3vrLpsK5/Ol0bDxSbXBkAtC3CCODFQh12PX7jEL169yglRgYr90SlbvvrJv32X9tVUcNdEgD+gTAC+IAr+nbVh7+8UjNG9ZAk/WPjIV3/l8/02T7euAHg+zwOI+np6Zo+fboSEhJksVi0Zs2a87ZftWqVrrvuOsXExCgiIkJjxozRhx9+2Np6gU4rzGHXf//gUv3zpyOV4AzS4RMV+snirzRveaYKS6vMLg8AWs3jMFJeXq6UlBQtWrSoRe3T09N13XXX6b333tPmzZt19dVXa/r06crMzPS4WADS+H4x+vCXV+quscmyWqS3tuRp4pNp+ufGg6p3s3orAN9jMS5i7WmLxaLVq1frxhtv9Oi6Sy65RLfeeqsee+yxFrUvKSmR0+mUy+VSREREKyoF/NP2oy49snqbth5xSZJSkiL1hxuHaEii0+TKAKDlv787fM6I2+1WaWmpoqKimm1TXV2tkpKSsw4A5xqS6NTqn4/V//3+JQp32LUl95RuWLRB//ftnSyWBsBndHgYefLJJ1VeXq5bbrml2TYLFy6U0+lsPJKSkjqwQsC32KwWzRyTrI//c4KmpyTIbUhLPs/RtU+m6f1t+Wy8B8DrdWgYWb58uX73u9/p9ddfV2xsbLPtFixYIJfL1Xjk5uZ2YJWAb4qNCNKzt6Xq5dkj1TM6RAUlVbp32Tf6yeKvtLuAu4sAvFeHhZHXX39dP/3pT/XGG2/o2muvPW9bh8OhiIiIsw4ALXNl/xh9+MCVmndNXwXardqwv0hTn/5Mj67epuKyarPLA4BzdEgYWb58ue688069+uqrmjZtWkd8JNCpBQXYNH/SAH08f4KmXRovtyEt+/Kwrvp/6/VS+gFV19WbXSIANPI4jJSVlSkrK0tZWVmSpJycHGVlZenw4cOSGh6xzJw5s7H98uXLNXPmTD355JMaPXq0CgoKVFBQIJfL1TY9ANCspKgQPXf7ZXrjnjEakhih0uo6/fd7uzXpz+n6cEcB80kAeAWPX+1dv369rr766nPOz5o1S0uXLtWdd96pgwcPav369ZKkq666Smlpac22bwle7QUuntttaOU3R/TEh3t0vLThcc2Y3tH6P98brMEJ/P8KQNtr6e/vi1pnpKMQRoC2U15dpxfWH9BLn2Wrps4ti0X6YWp3PXBtPyVFhZhdHgA/QhgBcF5HTlbof97frXe25kuSAmwW3T6qp+Ze3Vcx4Q6TqwPgDwgjAFokK/eU/t+Hu/X5/mJJUnCATbPHJes/xveRMyTA5OoA+DLCCACPfL6/SE98uEdbck9JkiKC7LpnQh/dNTZZIYF2c4sD4JMIIwA8ZhiG1u08pj+t3aO9x8okSV3DHLr/mr66bWQPBdo7fNFmAD6MMAKg1erdht7aclRPrdur3BOVkqQEZ5DmXNVHt4xIUlCAzeQKAfgCwgiAi1ZT59brGbl69uN9Kjz9OnBMuEM/G99Lt4/qqVAHj28ANI8wAqDNVNXW682MXL2Ylq2jpxrulESGBOinY3tp5hXJcgYz0RXAuQgjANpcTZ1ba7KO6vlP9+tgcYUkKdxh16wrkjV7XC9FhQaaXCEAb0IYAdBu6urdendbvp77dH/jRNfgAJtmjOqhu8Ymq3sXFk8DQBgB0AHcbkNrdx7Tok/3afvREkmSzWrRlCHddPf43hqWFGlugQBMRRgB0GEMw9D6vcf1t8+yGxdPk6TLk7vop+N667rBcbJZLSZWCMAMhBEAptiZV6LFG3L01pajqq1v+OelR1SIZo9N1s0jkngDB+hECCMATFVYUqWXNx7SK18e0qmKWkkNq7reNqqH7hjVk035gE6AMALAK1TU1GnlN0e1ZEOOcorKJUkWi3RV/xjdMbqnrhoQyyMcwE8RRgB4Fbfb0Me7C/XyxoP6bF9R4/nEyGDNGNVDt16epK5h7BYM+BPCCACvlVNUrle/PKQ3Mo7IVdnwCCfAZtH1Q+J1x6geGtkrShYLd0sAX0cYAeD1qmrr9c7WfL2y6ZCyTu8WLEn948J0y4gk/SA1UdHcLQF8FmEEgE/ZftSlZV8e0prMPFXW1kuS7FaLrhkYq5tHJOmqATEKsLFrMOBLCCMAfFJJVa3+lZWnFRm52nLE1Xi+a5hDP0hN0M0jktQ/LtzECgG0FGEEgM/bU1CqFZtztTrzqIrKahrPp3R36kcjkjR9aLwiQ9gPB/BWhBEAfqO23q31e47rzYxcfbK7UHXuhn+2AmwWTegfo+kpCbpucJxCAllQDfAmhBEAfqmorFprMo9qxeYj2l1Q2ng+OMCm6wbH6YaUBF3ZP0aBduaXAGYjjADwe/uOleqtLXn6V1aeDp+oaDzvDA7Q1Eu7aXpKgkb1imZRNcAkhBEAnYZhGNpyxKW3svL09tY8HS+tbvxZTLhDkwbH6foh3TS6dzRv5AAdiDACoFOqdxv6MrtYb23J03vb8lVSVdf4M2dwgK4d1BBMxvfrqqAAm4mVAv6PMAKg06upc2tjdrE+2J6vtTuOqbj832/khAbadNXAWF1/STddPTBWYewmDLQ5wggAfEu929DXB0/og+0F+nBHgfJdVY0/C7RZNap3lCYOjNU1A+PUI5odhYG2QBgBgGacmWPywfYCfbA9XweLK876eb/YMF0zKFYTB8bpsh6RsjPPBGgVwggAtIBhGDpwvFyf7D6mj3cVKuPQSdW7//3PojM4QBP6x+iagbEa368re+UAHiCMAEAruCpqlb7vuD7ZXahP9xTqVEXtWT8fkhih8f1iNL5vVw1P7iKHnUmwQHMIIwBwkerdhjIPn9RHuwqVvve4duaXnPXz4ACbRvWO0vh+MbqyX1f1jQ2TxcKaJsAZhBEAaGPHS6v1+f4ipe87rs/2FZ21nokkdYsI0pg+0RrTO1qje0crKSqYcIJOjTACAO3IMAztOVaqz/Y2hJOvck6ous59VpvEyGCN6h31rXDCWzroXAgjANCBqmrrlXHwpDZlF2tjdrG25J5q3NDvjO5dgjW6d7RG9orSiJ5d1KtrKHdO4NcIIwBgooqaurPCydYjrrPe0pGk6NBADe/ZRSOSu2hEcpSGJDjZ4A9+hTACAF6kvLpOXx88oU3ZJ7T50AltOeJSzXce6zjsVqUkRWpEzy4a3rOLUpIi1ZVXieHDCCMA4MWq6+q1/ahLXx88qYyDJ7X50Amd/M5rxFLDo51hSZGNx5BEJ3vqwGcQRgDAh5xZfG3zoRP6+uBJZeWe0v7CsnPa2a0WDYwP17CkSA3tHqlLE53qFxvGKrHwSoQRAPBxJVW12prrUlbuSWXlupSVe0pFZdXntHPYrRoYH6FLEyM0JMGpIYlO9Y8LZ/4JTEcYAQA/YxiGjp6qVFbuKWUdPqVtR13akVeisuq6c9oG2qwa0C1cQxIjNCg+QgO7RWhgfLgiggJMqBydFWEEADoBt9vQoRMVDcHkqEvbjrq0/ahLJVXnBhSpYe2TQfERGhQffjqkhKtndKhsVl4xRtsjjABAJ2UYho6crDx958Sl3fml2pVfojxXVZPtgwNs6hsbpn6xYeoXF67+cWHqHxeuxMhgWQkpuAiEEQDAWVwVtdpVUKLd+SXalV+qXQUl2lNQes7KsWeEBJ4JKQ0BpU9MmPrEhimpSzATZtEihBEAwAXVuw0dLC7XvmOl2nusTPsKy7TvWKmyj5erpr7pkBJgs6hndKh6dw1Vn9iGkNI7JlR9uobJGcKcFPwbYQQA0Gp19W4dLK5oDCl7CxsCSvbxsmbvpEhSVGigkqNDlBwdquSuoeoZHaJeXUPVMzpUzmCCSmdDGAEAtDm321Ceq1IHTgeTA8fLlH28XAeOl+lYybmvHX9bVGigep4OKkldgpUUFaKkqBD1iApRXEQQk2j9EGEEANChyqrrdLCoXIeKK3SwuLzx+5zich0vPX9QCbBZ1L1LiLp3CVaP0yElMTJYiV2C1T0yWF3DHEym9UEt/f1t78CaAAB+LMxh15DEhkXXvqu8uq4xpBwqrlDuyQrlnmg4jpysVG29oZyicuUUlTf5dwfarIqPDGoIKKdDSmJksBIig9XNGaQEZ7CCA1km31cRRgAA7S7UYdfghAgNTjj3v47r3YYKSqp0+Dsh5eipSh09WamCkirV1Lt1qLhCh4ormv0MZ3CA4p1BDUdksOIjGr52iwhSXIRDcc4ghTvssli4w+JtCCMAAFPZrJbGOx5jFH3Oz2vr3SpwVenoqUrlnQ4oR081HPmuKuWfqlR5Tb1clbVyVdZqd0Fps58VHGBrCCYRQaePhu9jI4IUE+ZQTLhDsREOQksHI4wAALxagM3aONm1KYZhqLS6TvmnqpTv+ndAyXdVKd9VpWMlDUdJVZ0qa+t1sLhCB89zh0Vq2O8nJvx0ODn9tWvYmSNQ0ae/jw4LJLi0AcIIAMCnWSwWRQQFKKJbgAZ0C2+2XUVNnQpLqhvCSWm1jp0OKgUlVTpeWq3jZdU6Xlqt0qo6Vde5deRkpY6crLzg5wfarIoOC2wMJ1EhgYoKDVSX0EBFhzZ8/+0jIiiAybjfQRgBAHQKIYF2JXe1K7lr6HnbVdbUq6isWoWl1f8OKSVVOl5Wo+KyahWX16iorFrFZTUqq65TTb278S5MS9isFkUGB8gZEqAuIYHqEhKgyLO+NnzvDAlQZHDg6a8BCgm0+e0dGMIIAADfEhxoO+9joW+rqq1vCCel1Sour1ZRWY1OltfoxLePitNfy2pUWl2nereh4vIaFZfXSGr67aGm2K0WRYYEKCK4IZw4Tx8RwQENd4aC7YoIOvdceFCAwoPsCvDiJfwJIwAAtFJQgK1x8m1L1NS5dbKipuEor9WpihqdrKjVyYqaxu/PfHVV1upURa1KKmtVU+9WndtQUVmNispqWlmrtSGYOOwKD/p3SDnz/Y3DEnVp93Nfy+4IhBEAADpIoN3a+CZPSxmGocra+sZw8u2QUlJ15mudXJXfPlfX+LPymnpJUlWtW1W11c0uQJeSFEkYAQAA57JYLAoJtCsk0K54Z8vuwHxbXb1bZdV1Kq06c9Q2fK2u/da5Og2Ia37yb3sjjAAA4MfsNqsiQwIVGRJodinN8ng2S3p6uqZPn66EhARZLBatWbPmgtekpaVp+PDhCgoKUu/evfXiiy+2plYAAOCHPA4j5eXlSklJ0aJFi1rUPicnR1OnTtX48eOVmZmpRx55RPPmzdPKlSs9LhYAAPgfjx/TTJkyRVOmTGlx+xdffFE9evTQX/7yF0nSoEGDlJGRoT/96U+66aabPP14AADgZ9r9peONGzdq0qRJZ52bPHmyMjIyVFtb2+Q11dXVKikpOesAAAD+qd3DSEFBgeLi4s46FxcXp7q6OhUVFTV5zcKFC+V0OhuPpKSk9i4TAACYpEOWY/vu8rWGYTR5/owFCxbI5XI1Hrm5ue1eIwAAMEe7v9rbrVs3FRQUnHWusLBQdrtd0dHnbhUtSQ6HQw6Ho71LAwAAXqDd74yMGTNG69atO+vc2rVrNWLECAUEBLT3xwMAAC/ncRgpKytTVlaWsrKyJDW8upuVlaXDhw9LanjEMnPmzMb2c+bM0aFDhzR//nzt2rVLS5Ys0eLFi/Xggw+2TQ8AAIBP8/gxTUZGhq6++urGP8+fP1+SNGvWLC1dulT5+fmNwUSSevXqpffee0+//OUv9dxzzykhIUHPPPMMr/UCAABJksU4M5vUi5WUlMjpdMrlcikiIsLscgAAQAu09Pd3h7xNAwAA0BzCCAAAMJVP7Np75kkSK7ECAOA7zvzevtCMEJ8II6WlpZLESqwAAPig0tJSOZ3OZn/uExNY3W638vLyFB4e3uyqra1RUlKipKQk5ebm+u3EWH/vI/3zff7eR3/vn+T/faR/rWcYhkpLS5WQkCCrtfmZIT5xZ8Rqtap79+7t9vdHRET45f/Avs3f+0j/fJ+/99Hf+yf5fx/pX+uc747IGUxgBQAApiKMAAAAU3XqMOJwOPTb3/7Wrzfl8/c+0j/f5+999Pf+Sf7fR/rX/nxiAisAAPBfnfrOCAAAMB9hBAAAmIowAgAATEUYAQAApvL7MPL888+rV69eCgoK0vDhw/XZZ5+dt31aWpqGDx+uoKAg9e7dWy+++GIHVdp6nvRx/fr1slgs5xy7d+/uwIpbLj09XdOnT1dCQoIsFovWrFlzwWt8aQw97Z+vjd/ChQt1+eWXKzw8XLGxsbrxxhu1Z8+eC17nK2PYmv752hi+8MILGjp0aOOCWGPGjNH7779/3mt8Zfwkz/vna+P3XQsXLpTFYtEDDzxw3nYdPYZ+HUZef/11PfDAA3r00UeVmZmp8ePHa8qUKTp8+HCT7XNycjR16lSNHz9emZmZeuSRRzRv3jytXLmygytvOU/7eMaePXuUn5/fePTr16+DKvZMeXm5UlJStGjRoha197Ux9LR/Z/jK+KWlpWnu3LnatGmT1q1bp7q6Ok2aNEnl5eXNXuNLY9ia/p3hK2PYvXt3/c///I8yMjKUkZGha665Rt///ve1Y8eOJtv70vhJnvfvDF8Zv2/7+uuv9dJLL2no0KHnbWfKGBp+bOTIkcacOXPOOjdw4EDj4YcfbrL9r371K2PgwIFnnbvnnnuM0aNHt1uNF8vTPn766aeGJOPkyZMdUF3bkmSsXr36vG18cQzPaEn/fHn8DMMwCgsLDUlGWlpas218eQxb0j9fH0PDMIwuXboYf/vb35r8mS+P3xnn65+vjl9paanRr18/Y926dcaECROMX/ziF822NWMM/fbOSE1NjTZv3qxJkyaddX7SpEn64osvmrxm48aN57SfPHmyMjIyVFtb2261tlZr+nhGamqq4uPjNXHiRH366aftWWaH8rUxbC1fHT+XyyVJioqKaraNL49hS/p3hi+OYX19vV577TWVl5drzJgxTbbx5fFrSf/O8LXxmzt3rqZNm6Zrr732gm3NGEO/DSNFRUWqr69XXFzcWefj4uJUUFDQ5DUFBQVNtq+rq1NRUVG71dpareljfHy8XnrpJa1cuVKrVq3SgAEDNHHiRKWnp3dEye3O18bQU748foZhaP78+Ro3bpyGDBnSbDtfHcOW9s8Xx3Dbtm0KCwuTw+HQnDlztHr1ag0ePLjJtr44fp70zxfH77XXXtM333yjhQsXtqi9GWPoE7v2XgyLxXLWnw3DOOfchdo3dd6beNLHAQMGaMCAAY1/HjNmjHJzc/WnP/1JV155ZbvW2VF8cQxbypfH77777tPWrVu1YcOGC7b1xTFsaf98cQwHDBigrKwsnTp1SitXrtSsWbOUlpbW7C9sXxs/T/rna+OXm5urX/ziF1q7dq2CgoJafF1Hj6Hf3hnp2rWrbDbbOXcICgsLz0l8Z3Tr1q3J9na7XdHR0e1Wa2u1po9NGT16tPbt29fW5ZnC18awLfjC+N1///1666239Omnn6p79+7nbeuLY+hJ/5ri7WMYGBiovn37asSIEVq4cKFSUlL09NNPN9nWF8fPk/41xZvHb/PmzSosLNTw4cNlt9tlt9uVlpamZ555Rna7XfX19edcY8YY+m0YCQwM1PDhw7Vu3bqzzq9bt05XXHFFk9eMGTPmnPZr167ViBEjFBAQ0G61tlZr+tiUzMxMxcfHt3V5pvC1MWwL3jx+hmHovvvu06pVq/TJJ5+oV69eF7zGl8awNf1rijePYVMMw1B1dXWTP/Ol8WvO+frXFG8ev4kTJ2rbtm3KyspqPEaMGKHbb79dWVlZstls51xjyhi229RYL/Daa68ZAQEBxuLFi42dO3caDzzwgBEaGmocPHjQMAzDePjhh42f/OQnje2zs7ONkJAQ45e//KWxc+dOY/HixUZAQICxYsUKs7pwQZ728c9//rOxevVqY+/evcb27duNhx9+2JBkrFy50qwunFdpaamRmZlpZGZmGpKMp556ysjMzDQOHTpkGIbvj6Gn/fO18bv33nsNp9NprF+/3sjPz288KioqGtv48hi2pn++NoYLFiww0tPTjZycHGPr1q3GI488YlitVmPt2rWGYfj2+BmG5/3ztfFrynffpvGGMfTrMGIYhvHcc88ZPXv2NAIDA43LLrvsrFfuZs2aZUyYMOGs9uvXrzdSU1ONwMBAIzk52XjhhRc6uGLPedLHP/7xj0afPn2MoKAgo0uXLsa4ceOMd99914SqW+bMa3TfPWbNmmUYhu+Poaf987Xxa6pvkoy///3vjW18eQxb0z9fG8PZs2c3/vsSExNjTJw4sfEXtWH49vgZhuf987Xxa8p3w4g3jKHFME7PSgEAADCB384ZAQAAvoEwAgAATEUYAQAApiKMAAAAUxFGAACAqQgjAADAVIQRAABgKsIIAAAwFWEEAACYijACAABMRRgBAACmIowAAABT/X8nBlK6xOt1KgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make an array of t values and plot\n",
    "tval=np.linspace(0,4)\n",
    "plt.plot(tval,g(tval))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical solution using *numpy*.\n",
    "\n",
    "The *numpy* library already loaded. Using a *lamdbified* rhs of the problem. Solve at the same times as the previously \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs=lambdify(t,deqn1.rhs,'numpy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[10], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mrhs\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m<lambdifygenerated-2>:2\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[0;34m(t)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_lambdifygenerated\u001b[39m(t):\n\u001b[0;32m----> 2\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;241m1\u001b[39m \u001b[38;5;241m-\u001b[39m \u001b[43mx\u001b[49m(t)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "rhs(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.11.3"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "16647ce0c82f44e260e72e9c8032e13b5a31dc8dd577de1db7c73c766b1ea8ae"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}