Solve_ODEs_with_scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve Differential Equations in Python\n",
    "\n",
    "![Solve with ODEINT](https://apmonitor.com/pdc/uploads/Main/solve_odeint.png)\n",
    "\n",
    "See https://apmonitor.com/pdc for [source code](https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations) and additional material."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Differential equations are solved in Python with the Scipy.integrate package using function __odeint__ or __solve_ivp__. Another Python package that solves differential equations is [GEKKO](https://gekko.readthedocs.io/en/latest/). See this [link for the same tutorial in GEKKO](https://apmonitor.com/pdc/Main/PythonDifferentialEquations) versus ODEINT. ODEINT requires three inputs:\n",
    "\n",
    "```python\n",
    "y = odeint(model, y0, t)\n",
    "```\n",
    "\n",
    "__model__: Function name that returns derivative values at requested y and t values as dydt = model(y,t)\n",
    "\n",
    "__y0__: Initial conditions of the differential states\n",
    "\n",
    "__t__: Time points at which the solution should be reported. Additional internal points are often calculated to maintain accuracy of the solution but are not reported."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Problem"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example of using _odeint_ is with the following differential equation with parameter $k=0.3$, the initial condition $y_0=5$ and the following differential equation.\n",
    "\n",
    "$\\frac{dy(t)}{dt} = -k \\; y(t)$\n",
    "\n",
    "The Python code first imports the needed Numpy, Scipy, and Matplotlib packages. The model, initial conditions, and time points are defined as inputs to _odeint_ to numerically calculate $y(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgrklEQVR4nO3deZxU5Z3v8c+vqnrf6KVoaPalQRBFNkUFNApxiRo10ZhJTDIxIZnEieZmuc7NJDeTmckkmVyzO8q4JO5bljHGLGqMIFGkQUQEZBNo1m7W3uit6rl/VIFtpxu6oU+dWr7v16tedbrOqT4/ThXfevqpc57HnHOIiEj6CfhdgIiIeEMBLyKSphTwIiJpSgEvIpKmFPAiImkq5HcBXVVUVLjRo0f7XYaISMpYsWLFPudcuKd1SRXwo0ePpqamxu8yRERShplt622dumhERNKUAl5EJE0p4EVE0pQCXkQkTSngRUTSlKdn0ZjZVqARiACdzrmZXu5PRETekYjTJN/jnNuXgP2IiEgXKd9F09oR4a4XN/PSRn2GiIh05XXAO+BPZrbCzBb2tIGZLTSzGjOrqa+v7/cOsoMBFi3ewpMrak+1VhGRtOJ1wJ/vnJsOXAZ83szmdd/AObfIOTfTOTczHO7xatvjCgSMOdUVvLRpP9GoJi8RETnK04B3zu2K39cBvwbO9mI/c8ZXsK+pjfV7Gr349SIiKcmzgDezAjMrOroMvBdY48W+5lbHWv4vbep/F4+ISLrysgVfCbxkZq8DrwK/c879wYsdDSnJZUJlIUv0RauIyDGenSbpnNsCTPXq93c3Z3yYB5dto7UjQm5WMFG7FRFJWil/muRRcydU0N4ZZfnWA36XIiKSFNIm4M8ZU0Z2MKBuGhGRuLQJ+PzsEDNGlbJ4g75oFRGBNAp4iHXTrN/TSF1jq9+liIj4Lq0Cfl78dMmlm9RNIyKSVgE/eWgxZQXZLNmggBcRSauADwSM88dXsGTTPpzTsAUiktnSKuAB5o6voL6xjbf2atgCEclsaRfwc6orADR8sIhkvLQL+KpBeYwfXMhiBbyIZLi0C3iIjS65bMt+WjsifpciIuKbtAz4eRMqaOuMsmLbQb9LERHxTVoG/DljyskKGos36qpWEclcaRnwBTkhpo8s1fnwIpLR0jLgAeZNCLN2dwP7mtr8LkVExBdpG/Bz46dLatgCEclUaRvwp1eVMCg/i8XqphGRDJW2AR8MGHOrw7y4oY5IVMMWiEjmSduAB5g/aTD7mtpZVavTJUUk86R1wF84cTChgPHs2jq/SxERSbi0DviSvCxmjy3n2bV7/C5FRCTh0jrgARZMrmRzfTNb6pv8LkVEJKHSPuAvnjQYgOfW7fW5EhGRxEr7gB9ems/kocU8u1YBLyKZJe0DHmLdNCu2HWS/rmoVkQySMQEfdfD8ep1NIyKZIyMC/vSqYqpKcnlO3TQikkEyIuDNjPmTK1mycZ8mARGRjJERAQ8wf1IlRzoimqtVRDJGxgT87LHlFOWEdLqkiGSMjAn47FCACyaGeW5dHVENPiYiGcDzgDezoJm9ZmZPe72vE1kwuZJ9TW28VnvI71JERDyXiBb8LcC6BOznhI4OPqZuGhHJBJ4GvJkNB94H3O3lfvqqJC+Lc8aW6apWEckIXrfgfwh8FYj2toGZLTSzGjOrqa+v97gcWDCpkk11Tby9r9nzfYmI+MmzgDezK4A659yK423nnFvknJvpnJsZDoe9KueY+ZMrATSEsIikPS9b8OcDV5nZVuBR4CIze9DD/fXJ8NJ8Jg0t5jlNAiIiac6zgHfO/ZNzbrhzbjRwA/Bn59xHvdpff7x3ciXLtx2grqHV71JERDyTMefBd3Xl1KE4B0+v3u13KSIinklIwDvn/uKcuyIR++qL8YOLmDS0mKde3+V3KSIinsnIFjzAVVOrWFV7iNoDLX6XIiLiiYwN+CunDgVQK15E0lbGBvzw0nxmjCrltwp4EUlTGRvwEOumWb+nkQ17G/0uRURkwGV0wF9+xlACBk+tUiteRNJPRgd8uCiH88dX8NTru3BOQwiLSHrJ6IAHuPLMKrYfaOH1HYf9LkVEZEBlfMBfMmUI2cGAumlEJO1kfMCX5GVxwcQwT6/eRUQzPYlIGsn4gIfY2TR1jW0se3u/36WIiAwYBTwwf1Il+dlBnRMvImlFAQ/kZQdZMLmSZ97YQ3tnr3OTiIikFAV83FVTqzh8pIMlG72fVUpEJBEU8HFzq8OU5GVpbBoRSRsK+LjsUIDLzxjCs2v30tLe6Xc5IiKnTAHfxfvPGkZLe4Q/rNF8rSKS+hTwXZwzpozR5fk8+mqt36WIiJwyBXwXZsaHZo3k1a0H2FTX5Hc5IiKnRAHfzQdmDCMUMB6vUSteRFKbAr6bwUW5XDxpML9csUPnxItISlPA9+CGWSPZ39zOc+v2+l2KiMhJU8D3YN6EMENLcnnk1e1+lyIictIU8D0IBozrZo7gpU37qD3Q4nc5IiInRQHfi+tnDgfgCX3ZKiIpSgHfi+Gl+cytDvN4zQ6NEy8iKUkBfxwfnjWCPQ2tvLihzu9SRET6TQF/HBdPqqS8IFtXtopISlLAH0d2KMAHZwzn+fV11DW0+l2OiEi/KOBP4PpZI4hEHU+u3OF3KSIi/aKAP4Fx4ULOHlPGY8trierLVhFJIZ4FvJnlmtmrZva6mb1pZv/i1b68dsOsEWzb38LLWzQpt4ikDi9b8G3ARc65qcBZwKVmNtvD/Xnm8jOGUpqfxX1Lt/pdiohIn3kW8C7m6Ji7WfFbSvZx5GYF+ejsUTy/fi9v72v2uxwRkT7xtA/ezIJmtgqoA551zi3zcn9eunH2KEIB476lb/tdiohIn3ga8M65iHPuLGA4cLaZTem+jZktNLMaM6upr6/3spxTMrg4l6umDuOJmh0cbunwuxwRkRNKyFk0zrlDwF+AS3tYt8g5N9M5NzMcDieinJN205wxHOmI8LBGmRSRFODlWTRhMxsUX84D5gPrvdpfIkyuKua8ceX84q9b6YhoMhARSW5etuCHAi+Y2WpgObE++Kc93F9CfGruGPY0tPLMG7v9LkVE5LhCXv1i59xqYJpXv98vF04YzNhwAXcveZurplZhZn6XJCLSI13J2k+BgPHJ88fwxs7DLN960O9yRER6pYA/CR+YPpxB+VncvWSL36WIiPRKAX8S8rKDfPScUTy7bi/b9uvCJxFJTgr4k/Sxc49e+LTV71JERHqkgD9Jg4tzuXJqFY/X1HL4iC58EpHko4A/BTfNGUNLe4RHdOGTiCQhBfwpOL2qhLnVFdy9ZAst7Z1+lyMi8i4K+FN0y8XV7Gtq58FXtvldiojIuyjgT9HM0WXMra7grhfViheR5HLcgI/PyvRBM/uRmT1hZveb2VfN7PREFZgKbp1fzf5mteJFJLn0GvBm9k1gKXAusAy4C3gc6AS+Y2bPmtmZiSgy2c0YpVa8iCSf47XglzvnZjjnvuSce9g595xz7mnn3O3OuSuBjwDZCaoz6d06fwL7m9t54GW14kUkOfQa8M653wGY2XXd15nZdc65OudcjZfFpZIZo0qZNyHMXYu30NymVryI+K8vX7L+Ux8fy3i3zq/mQHM7D6gvXkSSQK/DBZvZZcDlwDAz+3GXVcXE+uGlm+kjS7lgQphFi7dw4+xRFOR4NhqziMgJHa8FvwtYAbTG74/engIu8b601HRLvBV/v/riRcRnvTYxnXOvA6+b2UPOOQ220kfvtOI3c+O5oyhUK15EfHK80yR/a2ZX9rJurJl9y8w+6V1pqevW+dUcbOng/pe3+l2KiGSw43XRfBqYC6wzs+Vm9oyZvWBmbxM7J36Fc+7ehFSZYqaNLOU9E8Pc+ZfNHGxu97scEclQx+ui2QN81cxqgZeAXOAIsME515Kg+lLWbZdN4rIfLeZHz2/km1fpwl8RSby+nCZZCTwBfBEYQizk5QQmDiniw2eP5IFXtrGprtHvckQkA50w4J1z/wxUA/cAnwA2mtm3zWycx7WlvP+1YAL5WUH+/Xfr/C5FRDJQn0aTdM45YE/81gmUAk+a2fc8rC3llRfmcPNF43nhrXoWb6j3uxwRyTAnDHgz+4KZrQC+R2zwsTOcc/8AzAA+4HF9Ke8T549mZFk+//a7tXRGon6XIyIZpC8t+ArgWufcJc65J46eE++ciwJXeFpdGsgJBfk/l5/Ghr1NPLq81u9yRCSD9KUP/hvOuR4vy3TOqXO5Dy45fQjnjCnj9mc30NCqa8ZEJDE0o1MCmBlfv2IyB1va+emfN/ldjohkCAV8gkwZVsIHpw/nvqVvs21/s9/liEgGUMAn0FcumUhWMMB/PLPe71JEJAMo4BNocHEun7twHH94cw8vrK/zuxwRSXMK+AT79LyxVA8u5Gu/foMmzfwkIh7yLODNbER8cLJ1Zvammd3i1b5SSU4oyHc+cCa7G1r5/h/f8rscEUljXrbgO4EvOecmAbOBz5vZZA/3lzJmjCrlY7NH8YuXt7Jy+0G/yxGRNOVZwDvndjvnVsaXG4F1wDCv9pdqvnLpaQwpzuW2X66mvVNXuIrIwEtIH7yZjQamAct6WLfQzGrMrKa+PnPGaynMCfFvV09hw94m7nxxs9/liEga8jzgzawQ+CVwq3Ouoft659wi59xM59zMcDjsdTlJ5eJJlVw5tYqf/nmThhQWkQHnacCbWRaxcH/IOfcrL/eVqv7vlZPJzwly2y/fIBp1fpcjImnEy7NojNgY8uucc7d7tZ9UV1GYwz+/bzI12w7y0Kvb/S5HRNKIly3484EbgYvMbFX8drmH+0tZH5g+jLnVFXz39+vZeUgTZonIwPDyLJqXnHPmnDvTOXdW/PaMV/tLZWbGt685A4AvPrpK48aLyIDQlaxJYkRZPv969em8uvUAP31BI06KyKlTwCeRa6YN59ppw/jx8xtZvvWA3+WISIpTwCeZb109hRFl+dzyyGscbtHkICJy8hTwSaYwJ8SPb5hGXWMbt/1qNbH5zkVE+k8Bn4SmjhjEly+ZyO/X7NE8riJy0hTwSWrh3LHMGV/Bv/z2TV3lKiInRQGfpAIB4/brp5KfHeLmh1+jtSPid0kikmIU8ElscHEu37/uTNbvaeSbT72p/ngR6RcFfJK76LRKPv+ecTy6vJb7X97mdzkikkIU8CngSwsmMn/SYL719Fr+ummf3+WISIpQwKeAQMD4wYfOYmxFAZ97eCXb97f4XZKIpAAFfIooys3i7o/PxDn41P3LNWG3iJyQAj6FjCov4I6PTGdzfTNffGyVxo8XkeNSwKeY88dX8PX3TeLZtXv5wXMb/C5HRJJYyO8CpP8+ft5o1u1u5Cd/3kR1ZRFXTa3yuyQRSUIK+BRkZnzr6tPZsq+JLz/+OmX52cyprvC7LBFJMuqiSVE5oSB3f2wWY8MFLHyghte2H/S7JBFJMgr4FFaSn8X9nzybisIc/v7ny9mwV2PWiMg7FPApbnBxLg/edA7ZwQA33rOM2gM6R15EYhTwaWBkeT7333Q2R9oj3HjPMuob2/wuSUSSgAI+TZw2pJj7/v5s9ja08fF7X6WhVbNBiWQ6BXwamTGqlDtvnMHGukY+ed9yGhXyIhlNAZ9mLpgQ5kc3TGNV7SE+evcyDrW0+12SiPhEAZ+GLj9jKHd+dAbrdjdyw6JX2NekPnmRTKSAT1PzJ1dy7ydmsW1/C9ff9TJ7Drf6XZKIJJgCPo3Nqa7g/pvOpq6hjevu+qtOoRTJMAr4NDdrdBkPfeocGo50ct2dL7O5vsnvkkQkQRTwGWDqiEE8unA2ndEoH7rrZQ1rIJIhFPAZYtLQYh77zLnkZQe5YdErPL16l98liYjHFPAZZFy4kN987nzOHF7CzQ+/xo+f34hzmjREJF0p4DNMeWEOD37qHK6ZNozbn93AFx9bRWtHxO+yRMQDngW8md1rZnVmtsarfcjJyQkFuf36qXz5vRP4zapdfPTuZezXufIiacfLFvzPgUs9/P1yCsyMmy+q5md/N503dh7m6juWsnZXg99licgA8izgnXOLgQNe/X4ZGO87cyiPfeZc2jujXH3HUh5atk398iJpwvc+eDNbaGY1ZlZTX1/vdzkZ6awRg3jmC3OZPbacr/16Df/4yGsaqEwkDfge8M65Rc65mc65meFw2O9yMlZ5YQ4//8QsvnLJRJ55YzdX/uQl1uw87HdZInIKfA94SR6BgPH594zn0YXn0toR5dr/+isPvKIuG5FUpYCXv3H2mDJ+94U5nDu2nK//Zg2fvn8FdQ0arEwk1Xh5muQjwMvARDPbYWY3ebUvGXjlhTnc94lZfO3ySSzZWM+CHyzmVyt3qDUvkkIsmf7Dzpw509XU1PhdhnSzpb6Jrz65mpptB7notMF8+5ozGFKS63dZIgKY2Qrn3Mye1qmLRk5obLiQxz5zLl+/YjJ/3byPBT94kcdratWaF0lyCnjpk2DAuGnOGP5wyzwmDSnmq0+u5sZ7XmVTXaPfpYlILxTw0i+jKwp4dOFsvvX+01m94xCX/nAJ//r0Whp03rxI0lHAS78FAsbHzh3NC1++kOtmjuDepW9z0ff/wuM1tUSj6rYRSRYKeDlp5YU5/Me1Z/DU5+cwsiyfrz65mmvuWMpKTSgikhQU8HLKzhhewi//4Tx++KGz2H24lWvv+Cufvr+Gdbs1eJmIn3SapAyo5rZO7n3pbRYt2UJTWydXnFnFrfOrGRcu9Ls0kbR0vNMkFfDiiUMt7SxavIX7lm6lrTPCB6YP5wsXVzOiLN/v0kTSigJefLOvqY07XtjMg8u2EY06rjqrioXzxnLakGK/SxNJCwp48d2uQ0dYtHgLjy2v5UhHhAsnhvnMvHHMHluGmfldnkjKUsBL0jjY3M6Dr2zjFy9vZV9TO2cOL2HhvLFcevoQQkF95y/SXwp4STqtHRF+tXIn/71kC2/va2ZIcS4fmjWCG84ewdCSPL/LE0kZCnhJWpGo48/r63ho2TZe3FCPARdPquQj54xkXnWYQEDdNyLHc7yADyW6GJGuggFjweRKFkyupPZAC4+8up3Ha2p5du1ehpfmcd2MEVw9rYpR5QV+lyqSctSCl6TT3hnlT2v38PCy7by8ZT/OwfSRg7hm2jDed2YVZQXZfpcokjTURSMpa9ehIzz1+i5+vXInb+1tJBQwLpgQ5qqzqnjPaYMpzs3yu0QRXyngJS2s293Ab17byf+s2sWehlaygsZ54yq45PQhLJhcSbgox+8SRRJOAS9pJRp1vFZ7kD+s2cMf39zL9gMtmMHMUaW8d/IQLpwYZvzgQp1fLxlBAS9pyznH+j2N/PHNPfxhzR7W74lNQFJVkssFE8NcMCHMeeMr1JUjaUsBLxlj56EjLN5Qz4tv1bN00z4a2zoJBozpIwdx7thyZo8tZ/qoUnKzgn6XKjIgFPCSkToiUV7bfogXN9SxeMM+3tx1mKiD7GCAqSNKmD22nHPGlDNt5CAKcnTGsKQmBbwI0NDawYqtB3lly35e2bKfNbsaiEQdAYOJQ4qZNnIQ00YMYvqoUsaUF+giK0kJCniRHjS2drBi20FWbj/Ea9sPsqr2EI2tnQCU5GVx5vASTq8qYcqwYqZUlTCyLF+hL0lHV7KK9KAoN4sLJw7mwomDgdjZOZvrm3ht+yFWbj/IGzsPc89LW+iIxBpBRTkhJlcVM7mqmImVRUwYUsSEyiIK1b0jSUrvTJG4QMCoriyiurKI62eNAGJX1W7Y28ibuw6zZmcDa3Yd5rHltbS0R449b9igPCYOKaK6spBx4ULGhQsYW1FIqa64FZ8p4EWOIzsUYMqwEqYMK+FDs2KPRaOOnYeOsH5PIxv2NvJW/H7JxvpjrX2A0vwsxoYLGVtRwKjyfEaWFzCyLJ9RZfkMys/SefriOQW8SD8FAsaIsnxGlOWzYHLlscc7I1F2HDzCln1NbKlvZnN9M1vqm3hxQz11jW3v+h1FOSFGluczbFAew0rzYveD8qiK/1xekK0PADllCniRARIKBhhdUcDoigIuOu3d6460R6g92ML2/S1sO9BC7YEWtu1vZuv+ZpZu2kdzly4fiJ3KObg4hyHFuVSW5DKkOHYbXJxDuCiHwUU5hItyKc4N6YNAeqWAF0mAvOwgEypjX8p255yj4UgnOw8did0OtrC7oZW9h1vZ09DK2l0N/HldHUc6In/z3OxQgHBhDhVFOZQXZFNWkH3svqwgm/LCbErzsxmUn01pfhbFuVk6EyiDKOBFfGZmlORnUZKfxeSqnicjd87R0NpJXUMr9U1t1DfGb12W9za0sm53A/ub22nvjPayr9gpoKX52RTnZVESvxXnhmL3ebEPgaLcEIW5IYpzQxTmxH4uyAlRmBMiqA+IlOFpwJvZpcCPgCBwt3PuO17uTyRdmdmxMK7u4a+ArpxzNLdH2N/Uxv7mdg63dHCwpZ2DLR0cjt8fbGnn8JEODh/poPZACw3x5c7oia+LycsKxsM+dl+QEyI/O0hBdoi87CAF2UHysmOP5WcHyc0KkpcVJC87fsuKPZabFSA31GU5K0hOKKAupwHkWcCbWRD4GbAA2AEsN7OnnHNrvdqniMQ+DArjre3+zITlnKOlPUJjayeNrR00tnUeW25q7aSpLXZrbuukqS1C87HlTg40t7Pj4BFa2jppbo9wpD1Ce6TnvyJOJDsUICcUICcUjN8Hjj2WFYwtZ4cCZHdbProuKxggO2hkBQOEggGy4suxn42soBEKxB4PBuKPxe9DASMUDBAKGMGAvbNN/Oej939zs9h9sn04edmCPxvY5JzbAmBmjwLvBxTwIknIzI61yIeU5J7y7+uIRDnSEaG1PcKRjvgtHv6tnRFaO6K0drxz39YZu2+PvPNzW0eUts7Ycnv81tYZpbG1k45I/LH4fUckSkfEHfvZDwGLTUMZsHeCPxD/EIg9BgHrsj5gmEFFQQ6Pf/bcAa/Hy4AfBtR2+XkHcE73jcxsIbAQYOTIkR6WIyKJdLTV7MdQzc45IlFHR8TREY3S0RmlM+poj993xj8MItHY+s5I7LHO6NHnRYlEXWzbaJRIFCLRd9Z3xp8bie/nXTfniHZfdo5I9J26jj4edRBxjiKProb2MuB7+lvlbzr4nHOLgEUQG4vGw3pEJEOYWazLJQh5ZO7Q0AEPf/cOYESXn4cDuzzcn4iIdOFlwC8Hqs1sjJllAzcAT3m4PxER6cKzLhrnXKeZ3Qz8kdhpkvc65970an8iIvJunp4H75x7BnjGy32IiEjPvOyiERERHyngRUTSlAJeRCRNKeBFRNJUUk26bWb1wLaTfHoFsG8Ayxkoqqt/VFf/qK7+Sce6Rjnnwj2tSKqAPxVmVtPbzOJ+Ul39o7r6R3X1T6bVpS4aEZE0pYAXEUlT6RTwi/wuoBeqq39UV/+orv7JqLrSpg9eRETeLZ1a8CIi0oUCXkQkTaVUwJvZpWb2lpltMrPbelhvZvbj+PrVZjY9QXWNMLMXzGydmb1pZrf0sM2FZnbYzFbFb99IUG1bzeyN+D5relif8GNmZhO7HIdVZtZgZrd22yYhx8vM7jWzOjNb0+WxMjN71sw2xu9Le3nucd+PHtT1n2a2Pv46/drMBvXy3OO+5h7U9U0z29nltbq8l+cm+ng91qWmrWa2qpfnenm8esyGhL3HnHMpcSM25PBmYCyQDbwOTO62zeXA74nNJjUbWJag2oYC0+PLRcCGHmq7EHjah+O2Fag4znpfjlm313UPsYs1En68gHnAdGBNl8e+B9wWX74N+O7JvB89qOu9QCi+/N2e6urLa+5BXd8EvtyH1zmhx6vb+v8HfMOH49VjNiTqPZZKLfhjk3g759qBo5N4d/V+4H4X8wowyMyGel2Yc263c25lfLkRWEdsTtpU4Msx6+JiYLNz7mSvYD4lzrnFwIFuD78f+EV8+RfA1T08tS/vxwGtyzn3J+dcZ/zHV4jNkpZQvRyvvkj48TrKzAy4HnhkoPbXV8fJhoS8x1Ip4HuaxLt7iPZlG0+Z2WhgGrCsh9XnmtnrZvZ7Mzs9QSU54E9mtsJiE5x35/cxu4He/+P5cbwAKp1zuyH2HxQY3MM2fh+3TxL7y6snJ3rNvXBzvOvo3l66G/w8XnOBvc65jb2sT8jx6pYNCXmPpVLA92US7z5N9O0VMysEfgnc6pxr6LZ6JbFuiKnAT4DfJKis851z04HLgM+b2bxu6307ZhabyvEq4IkeVvt1vPrKz+P2NaATeKiXTU70mg+0/wLGAWcBu4l1h3Tn5//ND3P81rvnx+sE2dDr03p4rF/HLJUCvi+TePs20beZZRF7AR9yzv2q+3rnXINzrim+/AyQZWYVXtflnNsVv68Dfk3sz76u/Jwc/TJgpXNub/cVfh2vuL1Hu6ni93U9bOPLcTOzjwNXAB9x8Y7a7vrwmg8o59xe51zEORcF/ruX/fl1vELAtcBjvW3j9fHqJRsS8h5LpYDvyyTeTwEfi58ZMhs4fPTPIC/F+/juAdY5527vZZsh8e0ws7OJHfv9HtdVYGZFR5eJfUm3pttmvhyzuF5bVn4cry6eAj4eX/448D89bJPwSeXN7FLgfwNXOedaetmmL6/5QNfV9Tuba3rZX8KPV9x8YL1zbkdPK70+XsfJhsS8x7z45tirG7EzPjYQ+2b5a/HHPgt8Nr5swM/i698AZiaorjnE/nRaDayK3y7vVtvNwJvEvgl/BTgvAXWNje/v9fi+k+mY5RML7JIujyX8eBH7gNkNdBBrMd0ElAPPAxvj92XxbauAZ473fvS4rk3E+mSPvsfu7F5Xb6+5x3U9EH/vrCYWQEOT4XjFH//50fdUl20Tebx6y4aEvMc0VIGISJpKpS4aERHpBwW8iEiaUsCLiKQpBbyISJpSwIuIpCkFvGQsMxtkZp+LL1eZ2ZN+1yQykHSapGSs+NggTzvnpvhdi4gXQn4XIOKj7wDj4uOEbwQmOeemmNkniI3uFwSmEBtbJRu4EWgDLnfOHTCzccQuEgsDLcCnnXPrE/2PEOmNumgkk91GbKjis4CvdFs3Bfg7YuOS/DvQ4pybBrwMfCy+zSLgH51zM4AvA3ckomiRvlILXqRnL7jY+N2NZnYY+G388TeAM+OjA54HPBEfMgcgJ/FlivROAS/Ss7Yuy9EuP0eJ/b8JAIfirX+RpKQuGslkjcSmUes3FxvT+20zuw6OzW07dSCLEzlVCnjJWM65/cDS+ETN/3kSv+IjwE1mdnQkwgGbgk5kIOg0SRGRNKUWvIhImlLAi4ikKQW8iEiaUsCLiKQpBbyISJpSwIuIpCkFvIhImvr/fPaFYp814RAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# function that returns dy/dt\n",
    "def model(y,t):\n",
    "    k = 0.3\n",
    "    dydt = -k * y\n",
    "    return dydt\n",
    "\n",
    "# initial condition\n",
    "y0 = 5\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,20)\n",
    "\n",
    "# solve ODE\n",
    "y = odeint(model,y0,t)\n",
    "\n",
    "# plot results\n",
    "plt.plot(t,y)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('y(t)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Input Arguments\n",
    "\n",
    "An optional fourth input is _args_ that allows additional information to be passed into the _model_ function. The _args_ input is a tuple sequence of values. The argument $k$ is now an input to the _model_ function by including an addition argument.\n",
    "\n",
    "```python\n",
    "y = odeint(model, y0, t, args)\n",
    "```\n",
    "\n",
    "__args__: Additional inputs to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+z0lEQVR4nO3dd3hUZfbA8e+bhBASQg099C69g2AFQYoUFQsKuoqoqCurqyL+VhEL6upakaLCAioqIAisBQRRiiAdKULohJrQEgghJHl/f5x0EkiZmZvMnM/zzJPJnTtzTybJuXfecl5jrUUppZT38XM6AKWUUu6hCV4ppbyUJnillPJSmuCVUspLaYJXSikvFeB0ABmFhYXZWrVqOR2GUkoVGevWrYu21lbI7rFCleBr1arF2rVrnQ5DKaWKDGPM/pwe0yYapZTyUprglVLKS2mCV0opL1Wo2uCVUio7Fy9eJDIykvj4eKdDcUxQUBDh4eEUK1Ys18/RBK+UKvQiIyMJDQ2lVq1aGGOcDsfjrLWcOHGCyMhIateunevnuTXBG2P2AbFAEpBorW3rzuMppbxTfHy8zyZ3AGMM5cuXJyoqKk/P88QV/A3W2mi3HiEhAQICwE+7FJTyVr6a3FPl5+cv+hnx0CG47jouvPKW05EopVSh4u4Eb4GFxph1xphh2e1gjBlmjFlrjFmb148fAAdWLaNa3ThC9y4g7rtFBY1XKaWytW/fPpo2bZrn561bt45mzZpRr149/v73v5PdGhwnTpzghhtuoGTJkjz++OOuCBdwf4LvbK1tDfQEHjPGXJt1B2vtJGttW2tt2woVsp1te1mV+g3gcN2tXKz5O18MmQ379hU8aqWUcpFHH32USZMmERERQUREBD/++OMl+wQFBfHKK6/w9ttvu/TYbk3w1trDKV+PA3OA9q4+RvGA4owImw/jtvNZ7GC49VY4f97Vh1FKqTR79uyhVatWrFmz5rL7HTlyhJiYGDp16oQxhiFDhjB37txL9gsJCaFLly4EBQW5NE63dbIaY0IAP2ttbMr97sAYdxzr1Qd6MvkFy2rbgM0bEmn+6KMwZQr4eKeMUl7JXf/XuVy+dMeOHdx1111MmTKFEiVK0LJly2z3W7p0KYcOHSI8PDxtW3h4OIcOHXJFtLnizlE0lYA5KT2/AcCX1tpLP5u4QEgI3DvY8PHHMMl/OB9NfRQ6dIBHH3XH4ZRSPioqKop+/foxe/ZsmjRpAsDGjRtz3D+79nZPjgZyW4K31u4BWrjr9TOKjIkk8ZpPYHMQ09c9zVvnnyL4ySehZUvo1MkTISilPCWXV9ruULp0aapXr86KFSto0qQJO3bs4M4778x236VLlxIeHk5kZGTatsjISKpWreqpcL1jJmvshVgm7RhD8atrEbP8eZbe8h96zX8Ubr8d1q+HSpWcDlEp5QUCAwOZO3cuPXr0oGTJkgwaNOiyV/BlypQhNDSUVatW0aFDB6ZNm8YTTzzhsXi9IsHXK1ePkZ1HEnyuKXdHWOrVeABu/BxWrIA77oCff4Y81G9QSqmchISEsGDBAm666SZCQkLo16/fZfcfP348999/P+fPn6dnz5707NkTgHnz5rF27VrGjJGuyVq1ahETE0NCQgJz585l4cKFXHXVVQWK1WTXRuSUtm3bWpct+HHkCLRuDUePwvDhMG6ca15XKeVx27dvp3Hjxk6H4bjs3gdjzLqcysAU/Zms2bAW9idUgTlzIDAQPv4YJk50OiyllPIor0nwJ+JOsGDnAuZvWULz5tC8OZxr1hEmTZIdHn8cfvvN2SCVUsqDvCbBLz+wnFtm3MJHG98kNBRiYuCbb4D77oOnnoLERLjtNp3pqpTyGV6T4FtUbkG3Ot3oXL0zw1Kq3qS1yrz5JvToAdHR0K8fnD3rWJxKKeUpXpPga5WpxaLBi3jxuhe54w4oXRpWr4ZNm5BSwjNmQP36sHkz3H8/JCc7HbJSSrmV1yT4jIKDYfBguZ/aBE/ZsjBvHpQqBbNnwyuvOBafUkp5glcleGstR88e5XDsYR5+WLZNnw6xsSk7NGokV/LGwOjRKY30Sil1Ze4sF7xo0SLatGlDs2bNaNOmDUuWLHFFyN6V4D9Y/QFV3qnCG8vfoGlTuOYaKF4ctm7NsFOvXvBWyuIg990Hq1Y5EqtSyjfkplxwWFgY8+fP588//2Tq1KkMTm2CKCCvSvANyjegVPFSaWfIadPg4EHo2DHLjk8/DQ89BPHx0umqI2uUUnng6nLBrVq1SqtR06RJE+Lj47lw4UKB4/SKUgWputftzunnTqdVa6tVK4cdjZGZrXv3ShmDPn2krEHp0h6LVSmVf5cryDhxImkj6SZNIq25Njv5mcjv7nLBs2fPplWrVhQvXjzvwWXhVQne388/2+0xMbBwodQeS1OsGMycCVdfLW04d9wB//ufjLhRSqlsuLtc8NatW3nuuedYuHBhgWMFL0vwGVlrMcZw8SI0bCglabZsgZTfiShTRpJ6hw5yBnjiCSlroAuFKFWo5fbKe9iw9Kt5V3BnueDIyEgGDBjAtGnTqFu3rkvi9ao2eID3V71PjXdrMG6NFBcrVgwGDEh57P1snlC7NsydK72xEybAe+95KlSlVBGTWi542rRpfPnllzRs2JCNGzdmeytTpgxVqlRJKxdsrWXatGnZVp88ffo0vXv3ZuzYsXTu3Nll8XpdgrdYDsYcZFvUtrRtf/+7fJ0+HU6cyOZJV18tS/yBdMBm0wmilFKQXi743Xff5bvvvrvi/uPHj2fo0KHUq1ePunXrZioX/OKLLwLw0UcfsWvXLl555RVatmxJy5YtOX78eIFj9bpywcfPHedM/BnqlK2TqU2+Z0/48Ud4/XV4/vkcnvzKK/DiixAUBEuW6GpQShUSWi5Y+Hy54IohFalfvv4lHa4jRsjXcePg4sUcnvx//wdDh8rwyT59YMcOt8aqlFLu5HUJPifdu8tE1kOH4Ntvc9jJGBg/Hnr3hpMn4eabpXdWKaWKIK9M8FM2TOHWr29l+YHladuMkbb4zp2hQoXLPDkgAL7+Gtq1kwlQvXtnqHWglFJFh1cm+DWH1zDnrzn8ceiPTNsffhiWL4cbb7zCC4SEwIIFULeuLNo9cOBl2nWUUqpw8soEP6TFEL649Qtua3xbpu1+eflpK1aUXtkKFeCnn2QwbSHqkFZKqSvxygTfMbwjg5oNomaZmtk+vmkTPPigtMdfVr16MhEqOBj++1/phFVKqSLCKxP8lbz2GkyeDO++m4ud27WTssL+/jLGMldPUkp5G3eWC963b19aXZuWLVvyyCOPuCJk703wP+/5mbHLxnL07KWjYEaOlK8TJshgmSvq3Tt9ItRTT8HUqa4LVCnl1XJTLhigbt26abNgJ0yY4JJje22Cf3PFm4xaMoq1hy+dONW6tQybPHcOPvooly84eHB6GYMHH4RczGBTSnknV5cLdhevTfC3N76dER1GUL1U9WwfT53N+sEHkuhz5ckn4V//gqQkuPNOWLrUJbEqpfLGvGwwL2cuCnjLjFswLxvm75iftm3SukmYlw3D5qdXHDscexjzsqHqO9kX/bqSHTt2cNtttzFlyhRKlSqV1qyS9Xb69Ok8lQveu3cvrVq14rrrrmPZsmX5ii0rr60m+XDbyxSBBq67TopIrl4Nn34quTtXXn5Z2nXGjYO+feGXX6BNm4IHrJQq9NxVLrhKlSocOHCA8uXLs27dOvr378/WrVspVapUgeL12gR/JcbIVXz//tJv+vjj0o+aqyd+8IEk+RkzZLbr8uVSk1gp5RH2pUsT5/y751+ybVibYQxrk7lecNXQqtk+PzfcVS64ePHiaQt8tGnThrp167Jz507ats22xEyueXWCj46LZs2hNdxU9yYC/C79UW+5RUY+DhmSy+Seys9POlpPn4YffoBu3eC336T0sFLKa6WWC+7RowclS5Zk0KBBl72CL1OmTFq54A4dOjBt2jSeeOKJS/aLioqiXLly+Pv7s2fPHiIiIqhTp06B43V7G7wxxt8Ys8EYs8Ddx8qq46cd6fVlL7Ye35rt435+UkCyfv18vHixYjBrFnTpApGR0LWrfFVKeTV3lAv+7bffaN68OS1atOD2229nwoQJlCtXrsCxur1csDHmKaAtUMpa2+dy+7qiXHBGQ+cNZceJHfz7pn/TMTzrytuZWSu14sPC8niQmBi46Sb44w9o0AB+/RUqV85/0EqpS2i5YJHXcsFubaIxxoQDvYHXgKfceazsfNr301ztt2sX3HWXNK//8UceV+wrVUpKGtx4I2zcKM01S5fm40yhlFKu5e4mmveAZ4HknHYwxgwzxqw1xqyNiopyczjZq1YNDhyAtWth8eJ8vEDZsrKm61VXyQLe3btL+7xSSjnIbQneGNMHOG6tXXe5/ay1k6y1ba21bStcto5v/kWdiyIhKSHHx0uUSF8QZOzYfB6kQgX4+WepX7NhgywhpWWGlXKZwrT6nBPy8/O78wq+M9DXGLMP+Aq40RjzuRuPl61eX/Si4tsVs53RmtHw4RAaKiv1rV6dz4NVqSIfAWrWhFWrZFWoXM+iUkrlJCgoiBMnTvhskrfWcuLECYKCgvL0PLe1wVtrnweeBzDGXA/801p7r7uOl5OqoVUJLhbMwTMHIftJrQCUKQOPPQZvvCFzmb7/Pp8HrFFDzhLXXitDJ3v1koqUJUvm8wWVUqnjyZ1qxi0MgoKCMs2KzQ2PLLqdIcF7dBQNwOn405QMLJntOPisoqNlKPvZs7BiBVx9dQEOHBEBN9wgNYmvuUbOGJrklVIu5vii29bapVdK7u5SJqhMrpI7yMCXESOgXDkXDGmvX19G01SrBsuWyZX82bMFfFGllMo9ry02lp3cfFp59lnYuxfuuMMFB6xXT5J8eLgkee14VUp5kE8k+Bd/eZFa79Vi0Z5FV9w3NFSGtrtMvXpSkCw8XGrWaJJXSnmITyT4M/Fn2H9m/xVH0mQUFwdvvy25ucAyXsmvWCFJPibGBS+slFI580gna265o5MVYPfJ3cQnxtMorBH+frmrKvbuu7J4U/v2MuIxT7Nbcwxkt3S8HjwoSwH++KM0+CulVD453snqtLrl6tKkYpNcJ3eAYcOgYkUpXfC//7kqkLpSq6Z2bVizRpL98eMuenGllMrMJxJ8foSEpK/d+uKLkJxjsYU8ql1bxsc3bAibN8t4+RxWeFFKqYLwmQQ/a9ssBs4cyKLdV+5oTfXII1C1qlQemDPHhcGEh8uVfLNmsGOHjJPfu9eFB1BKKR9K8BuObGDWtln8si/3vaYlSsALL8j9l16SpVhdplIl6Xht106S+7XXSrJXSikX8ZkEP7DJQKb0m8LQ1kPz9LwHH5TqA1u3urAtPlW5clKgLHXRkGuvlWYbpZRyAZ8YRVNQc+dCYiLcequsAuVy587J4rA//wylS8OCBZL0lVLqCnx+FE1B9e8Pt9/upuQO0qM7f76cQc6ckRWiFnh8hUOllJfxqQS/+dhm3lz+Jr/t/y3frxERAadOuTCoVEFB8M03MHQoxMfLWWXaNDccSCnlK3wqwS/YuYCRi0cye9vsfD1/4kRZtOnVV10cWCp/f5g0CZ5/Xnp077tPZlwppVQ++FSC71anG8PbDqdHvR75en779pJ3P/wQ9uxxcXCpjIHXX4d33pHvn3oKRo2SVcGVUioPtJM1j+6/H6ZOhYEDpUXFraZNgwcekLPKAw/AhAlQrJibD6qUKkq0k9WFXn1VxsfPnAm//+7mgw0ZIkN4goJg8mTo21dryiulcs3nEnzMhRiW7F3CxqMb8/X88HB4+mm5//TTHmg56dNHSlqGhUlxsuuug6NH3XxQpZQ38LkEP2XDFLpO68pHf3yU79d49lkpRPb77zA7f/21edOxI6xcKcXK1q+HTp3gr788cGClVFHmcwm+U/VOtKvajvrl6uf7NUJDYcwYGDAAWrRwYXCXU7++JPn27WHfPlkwdvlyDx1cKVUUaSdrPlnrohrxeXXuHNx9t0yMKl4cpk+XHl+llE/STlY3yJjck5Ml73pESAh8+y08+ihcuCCLx44dq8MolVKX8NkEfzr+NDuiC169cft26NwZHnvMBUHlVkAAjBsHb74pZ5pRo+Bvf4OEBA8GoZQq7HwywS/dt5Tyb5Vn6Py8VZbMTkCA9HtOnerhJnFjpLd39mwIDpYAbroJoqM9GIRSqjDzyQTfsnJLAv0DCfQPJNkWbKmm+vXhuefk/mOPSdVJjxowAJYtk5VJfvtNRtxoXXmlFD6a4MsEleHUc6dYPGQxfqbgb8HIkVCrlpRyHzeu4PHlWevWsnhsq1aysHfHjrBkiQOBKKUKE59M8ABBAUEue63gYPjgA7n/4otw5IjLXjr3qlWTK/n+/eH0aejeXc422vmqlM/y2QSf6uT5ky55nVtukVtMjDSNOyIkRNrkn31W6tc8/jgMGyajbZRSPsdnE3xSchItJrSg4r8rcib+jEte8/33oUIFaNnSwQtnPz8ZXfPFF1LD5tNP4cYbtbyBUj7IZxO8v58/JQNLEhQQxPbo7S55zdq1Yf9+qVHjyCSojAYNkmE94eEyA7ZtWygik8iUUq7hswke4Ovbv+bkcyfpGN7RZa9ZokT6/YsXXfay+dOmjST1zp3h0CFZ5/Xzzx0OSinlKW5L8MaYIGPMH8aYTcaYrcaYl911rPwKLxVOoH+gy1/XWunfrFNHrugdVamSjKh56CFpix88GEaMKARnH6WUu7nzCv4CcKO1tgXQErjZGOO6S2UXK+h4+IyMkQq/kZHw8MOFYCBLYKCsN/jxx7JgyPvvww03ODTcRynlKW5L8Fakrk5RLOXmdKq7xNsr36bWe7WYs32OS1/3o4+gbFn46adCsna2MVK/5tdfZVLUihUybv63/C9ArpQq3NzaBm+M8TfGbASOA4ustauz2WeYMWatMWZtVFSUO8PJ1vmL59l/Zj8rDq5w6etWrgzvvSf3R4woRBfLnTpJbYXrr4djx2SEzX/+Uwg+ZiilXM0j5YKNMWWAOcAT1totOe3nRLngg2cOEh0XTYvKLVwyqzUja6F3b/jhB6koMHt2IRhdkyoxEV54Ad56S74fOBA++0yK3SuligzHywVba08DS4GbPXG8vKheujqtqrRyeXIHSeYTJ0rOnDMHZs1y+SHyLyBAxsvPmiUBzpwpQyk3b3Y6MqWUi7hzFE2FlCt3jDElgG6Az60zV706/PvfMvmpfv4XkXKf226DNWugaVPYuRM6dIBJk7TJRikv4M4r+CrAL8aYzcAapA1+gRuPl28bj27kzll3MvLnkW55/YceklpgLVu65eULrmFDWL0ahg6F+HgZ+nPPPRAb63RkSqkCcOcoms3W2lbW2ubW2qbW2jHuOlZBJSQl8M3Wb5i93T0raPv5yejEVA70JV9ZcDB88oksARgSAjNmyESpTZucjkwplU8+PZM1VesqrZnQewIL7nbvB4zYWLjzTrmSP33arYfKv3vvldmvzZpBRIQ02Ywfr002ShVBmuCBAL8AHm77MA3DGrr1OMHBcOAAHD4sQ9ILbc5s1EiabFJnvw4fLsOAdLUopYoUTfAe5O8vk56Cg+Grr6TgY6FVooR0tn71FZQuDd99By1a6EIiShUhmuBTxCfG887Kd7hv7n24c25A/frpi4M89hjs3eu2Q7nGnXfCxo1SsOzwYejWTZaw0gW+lSr0NMGnCPQPZOzysUzbNI3dp3a79VgPPAC33iqLgwwe7MA6rnlVqxYsXQqjR8vg/jfflIQfEeFwYEqpy9EEn8LP+DHmhjFM6z+NCsEV3HosY6T1I7UkzPjxbj2cawQEwEsvSe2amjWlI7ZlS5nJVWg7E5TybR4pVZBbTpQqcNLixfDttzIRKjjY6Wjy4PRp6XidMUO+79lTyhxUqeJoWEr5IsdLFajsde0qdeOLVHIHKFMGvvxSEnzZslJsp2lTKXeglCo0NMFn8Vf0X4xeOppFuxd59LgxMTBlikcPWXB33QV//gk9esDJk3DHHTKO/tQppyNTSnGFBJ+yKtPtxpj3jTEzjTHTjDHPGmOaeCpAT/s+4nte/vVlpmz0XLZNTISrr5bO16+/9thhXaNaNbmC//hj+SjyxRcySeqHH5yOTCmfl2OCN8aMBlYAnYDVwETgGyAReMMYs8gY09wTQXrSgEYDGNZ6GENbD/XYMQMCZOITSDmYv4paSbbUxUQ2boSOHWX911694G9/06t5pRyUYyerMaa3tfZ/OT7RmIpADWuty3pFfa2TNSNr4e675Qq+SROZSBoS4nRU+ZCYCO++C//6l8yCrVpVRtr06eN0ZEp5pXx1sqYmd2PMwGxecKC19rgrk7uvMwY+/VSqBGzdWkjWcs2PgAB45hm5mu/USSZH3XILDBki7fRKKY/JTSfr87nc5jWSbTKL9yzmuUXPuXQx7ispWVJWfQoJkabsCRM8dmjXa9QIli2Dd96BoCCpUtmkiSwwUiTPXEoVPZdrg+9pjPkQqGaM+SDD7b9IO7zXMhgemPcAb618iz8O/eHRY191lVTtBdixw6OHdj1/f3jqKVklqksXOHpUlgbs3x8iI52OTimvd7kr+MPAOiA+5WvqbR7Qw/2hOccYw+PtHuf5Ls9TuWRljx//7rtlgZDURbuLvPr14ddfZdB/aCjMmydnsnHjINlzn5CU8jVXnMlqjClmrb3oiWB8uZP1cmJipPnGzxtmLURGwuOPS3VKkHb6Tz6R5hulVJ7lq5PVGDPfGHNLDo/VMcaMMcY84KogVfY2bYLWreHll52OxEXCw9NXIK9cGX7/HVq1ghdegLg4p6NTyqtc7prwIeAaYLsxZo0x5ntjzC/GmL3ImPh11trJHonSIbEXYpm5dSbzdsxzLIZjx2DPHhgzxosqARgji31v3y7DhS5ehNdfl6v4BYVy2V6liqTcNNE8ASwHgoDzwE5rrVsutQpbE813f31H/6/7065qO/54yLOdrRn95z/w9NOyBsdvv0HbbD+MFWErV8pEqc2b5fv+/aUDomZNJ6NSqkgoaLGxSsBM4B9AZSTJ+4Sb6t7EjbVv5K6md7l1EZAr+cc/4MEH4fx56NdPJop6lauvhnXr5ExWsiTMnSudsG++qQuLKFUAuSoXbIwxQHfgb0BbpGTBZ9Zal66MUdiu4AuThAS46Sa5gm/TRr4WuSqUuXHokAyt/OYb+b5RI3j/feje3dm4lCqkClwu2MpZ4GjKLREoC8wyxrzlsijVZQUGyiSoOnXkYneec90C7lWtmtRr+OknGV75119SrbJ/f+mMUErl2hUTvDHm78aYdcBbSPGxZtbaR4E2wG1ujq9QOBx7mPFrxnPqvLOFs8LCYP58mRR6112OhuJ+3btLKeI335Rmm+++k2ab//s/OHfO6eiUKhJycwUfBtxqre1hrZ2ZOibeWpsM+EQFqcFzBjP8++H8LyLH2msec9VVUnI9VaFfz7UgiheHZ5+VKb2DB0vxstdek2abr77SkgdKXcEVE7y19kVr7f4cHtvu+pAKn9sb307v+r2pXaa206FksnEjNG4Mv/zidCRuVrUqTJsmC9i2bi2Tpe6+WzpnV61yOjqlCi1vmBvpdo+2e5QFgxbQuUZnp0PJ5IsvYNcuGVmzYYPT0XjA1VdLDYdPPoFKlSS5d+okyX7fPqejU6rQ0QRfhL3xhqySFxsr6177RB+kv7+sihIRIbNfg4KkuaZRI3j+eanroJQCNMHnyY7oHYxdNtbRMfEZ+ftLy0XXrjLjtXt3+eoTQkPh1VelfX7QIGmff+MNqFcPPvpIx88rhSb4XEu2ydw47UZGLRnF8gPLnQ4nTfHi8O230jS9e7eslBcb63RUHlSjhrRVrVolTThRUfDEE9Ib/fXXWq1S+TRN8LnkZ/wY3nY4D7R8gAohFZwOJ5NSpeD776FuXVi/HhYudDoiB3ToAMuXy9muYUM52911l2xfssTp6JRyRK5msubrhY2pDkxDyhskA5Oste9f7jk6k7Vg9uyBtWulXd6nJSbClCnw0ktw5Ihs69FDCpq1bu1sbEq5WIFnsuZTIvC0tbYx0BF4zBhzlRuP5/Pq1Mmc3Pfv99EWioAAeOghGWL02mvyEeenn6TGw+23y6K3SvkAtyV4a+0Ra+36lPuxwHagmruO5ynWWn7b/xtjfh3jdCiXtWmT5LOHH/bRJA9SrGfUKGmuefppGXEzezY0ayYTp3a7tJSSUoWOR9rgjTG1gFbA6mweG2aMWWuMWRsVFeWJcAok7mIcfb7sw0tLX2J7VOGd53X6tKyf8emnMHy4j0/6DAuDt9+WhP7oo3KF//nnMrTy4Yfh4EGnI1TKLdye4I0xJYHZwAhr7SWDlK21k6y1ba21bStUKFydl9kJCQzhyQ5PMqrLKEoHlXY6nBxdd53UrQkKgokTZWCJTyd5kBmxH38sQyvvv18+2kyaJL3Tjz4KBw44HaFSLuW2TlaQ9VyBBcBP1tr/XGl/7WR1vYULoW9fGSb+5JPw7ruyoJJCKlW+/LIMp7QWihWDBx6QCVO62IgqIhzpZE2pIf8ZsD03yV25R/fuMnIwMFDKqj/7rNMRFSKNGsGMGdLpevfdMvpm4kQpUzxsmJY/UEWeO5toOgODgRuNMRtTbr3ceDyPSkhKYNa2Wby8tPCvht2rl6xxXby4zP9RWTRuDF9+Cdu2wT33QFKS1LupVw/uu0/WjlWqCHJrE01eFaUmmqhzUVT7TzWSbBIH/3GQqqFVnQ7pig4dkvU01BXs2CHDK7/8UpK9MXDrrTIiR8fRq0LGqXHwXq1CSAWe6vQU7/Z4l5KBJZ0OJ1cyJvcNG+TiND7euXgKrYYNpcjPzp3wyCPSNj97tow77dlT1kssRBdGSuVEr+B9UFISNGkiF6pdu8oa1yWLxjnKGYcPy4LgEyakrybVsSM884zUavb3dzY+5dP0Cl5l4u8PM2dC5cqweDF06wYnTzodVSFWtaqMo9+/H158EcqVk+Jmt90m7fcTJ8L5805HqdQlNMEXUEJSAuPXjKffV/0KTRnh3GjWDJYtk9GAq1fLuPnUsi0qB+XLy7DKAwfggw+gVi2pS//II/JGvvIKREc7HaVSaTTBF5C1lld+e4V5O+axaM8ip8PJk3r1pABjo0awZYu0OmzZ4nRURUBIiMwci4iQxUZat5YyxS++CNWryxDLbducjlIpTfAFVTygOO90f4dvbv+GrrW7Oh1OnoWHy5V8p04yY1/Ls+RBQADceaeU8Fy8WMajxsfLEMsmTeDmm6XIWRH6ZKe8i3ayKkCakJcsgd69nY6kiNuxQ2aU/fe/6e3yV10lV/z33qu92crltJPVgy4kXnA6hHwpUSJzcl++XFocfLYSZX41bCj1biIjYexYGZu6bZvUuqlWDUaMkOGXSnmAJngXenvl24S/G87GoxudDqVA4uKkbPorr8jETh0rnw/lysHIkbB3r5RD6NxZFgR//305CfTsCf/7n55BlVtpgnehyJhIouOimbVtltOhFEhwMEydKutaf/UV3HCDDAVX+VCsmCwduHy5rKf44INS4vPHH6FPH6lkOXasD62WrjxJ2+Bd6FDMISJORnBdzeswXlCycfNmuOUWGRVYqZLUs+nSxemovMCJEzB5sjTlpBY0K1YMBgyQIZfXX68lP1WuXa4NXhO8uqyoKBko8ssvMmhk4kSpqKtcIDlZ6jlPmCDF+1Obaxo0kIVIBg+GIrBGgnKWdrI64Pi545yIO+F0GAVWoYLkoKeekvxTo4bTEXkRPz8ZSjl3rsySHT1aOmJ37pQlBqtVg4EDZahlUpLT0aoiSBO8G0zeMJma79XkjeVvOB2KSwQEwDvvyGCQbt3St8fFOReT1wkPh5dekiabuXOlfT4pSdrFbr4ZateWE8D+/Q4HqooSTfBu0LxSc+IT4zlz4YzTobhUw4bp9xctkpyzYIFz8XilgAApYDZ/vnR+vPYa1Kkjs9Beflne9G7dZE3Z1MJnSuVA2+DdZOvxrTSp2MTpMNzmnnukXDrI0O433pAFRZQbJCfD0qWygvq338r6iyDDnAYOlPVlu3TRjlkfpZ2syuWSk6XZZtQoWemudWsZUlm/vtORebnTp2UN2alT4fff07fXrSszZe+5R38JPkY7WR206+QuRi0eVaQqTeaGn5+UQ1++XIoqrl8vSf7zz52OzMuVKSMjbFaulEXDn39eOmN375YmnAYNoEMH+PBDOH7c6WiVw/QK3o0SkxOp834dDsYcZMZtM7ir6V1Oh+QWZ85IAcVvvpEa8zt2QKlSTkflQ5KSZBzrF19Ip+zZs7Ld3x969IBBg6BvX2nSUV5Hm2gcNHnDZJYfWM7rXV+ncsnKTofjNtbCZ59JH2DXrunbtFnYw+LiYN48SfY//ijtZyDFhvr0kVm1PXvK98oraIJXjhg5UsqvvPWWFlF0RFSUfKz66itpS0sVGgr9+8sMtm7dtHe8iNMEX0hYa4lPjKdEMe+/ejp6VK7m4+NllN+UKXDttU5H5cMOHpRkP2MGrFuXvr10aWm+GTgQbrpJ6uSoIkU7WQuBPaf20OPzHtz/3f1Oh+IRlSvLsqUtWsCePVJe5R//0KHbjqleXWbHrl0rM2VfeQWaN5cOlOnTJclXrCgjcebO1VlsXkKv4D3k4JmDNPyoIUEBQWx/bDuVSlZyOiSPSEiQXDJ2rPQF1qoF48bJ4keqENi5UzpmZ86EjRvTtwcHywzaW2+VhQLKlHEqQnUF2kRTSPwQ8QPtqrUjLDjM6VA8bu1aGDoUNm2SnDF7ttMRqUvs2iXJ/ttvYc2a9O0BAXDjjVLtsm9fqFrVuRjVJTTBq0IhMVGGZw8cKKVXQMqgh4XJiD5ViBw8KE01c+bAb79lLnbWtq0k+r59pZlHh0o5ShN8IWOt5bMNn9GqcivaVG3jdDiOSUqCjh0lP3z8seQNVQhFR0ttnO++k9KiqWvNgpQX7dtXmnGuv147aR2gCb6Q+XT9pzw0/yEalG/A5kc2UzzAN4ep7dolOeHQIUny990Hr78OVao4HZnK0fnzsHixjLWfP1+GS6UKDpZhl717SydL6sc05Vaa4AuZ8xfP0216N4a3Hc49ze9xOhxHxcbCmDGyVOnFizJeftQoGXGjF4OFXHKydK7Mmyfry2bspAUZQtWzp3TWXn21rFqlXE4TfCFkrfWKZf1cJSIC/vlPyRUAbdpIP5++RUXIoUPw/feS7H/+OfOY2NBQmeJ8881SPqFWLcfC9DaOJHhjzGSgD3DcWts0N8/xpQSf0eHYw/gZP68uZZBbP/8s5YeHDYO//122acmDIujCBemc/fFHuW3blvnxBg1kYtVNN8mq7lq8KN+cSvDXAmeBaZrgc7Zs/zIGfD2ANlXb8MM9P+BndO5ZavmUgAD5OnasFE987TUZtKGKoAMHZOnBH3+Us3hMTPpj/v5SATM14bdvr805eeDITFZr7W/ASXe9vreoV65e2v2zCWcdjKTwCAhIT+6JiTB+vKwc1bKlTLTcs8fR8FR+1KgBDz0kEyBOnIAVK2QJws6d5fGVK6XccZcuULasdNK+/TZs2JC+GLnKM7e2wRtjagELLncFb4wZBgwDqFGjRpv9Prjm5K6Tu6hTto5evefg+HEZXTN+vMyMDQiQSVMjR0LNmk5HpwosJkZWrFq0SEbobN+e+fFy5WS4VeqtSRNZkEABDnay5ibBZ+SLTTRZJdtkTsefplyJck6HUujs3y8XfdOmyUVdQIBc4DXN1V+XKjIOH5b69kuWSMLPetFXvrxUrrv+erjuOmjWzKcTvib4IiI+MZ4hc4aw88ROlv1tGaHFdYGG7GzfLu3xu3bJqnWpHbBHj0qRM+VFrIW9eyXZ//qrXOlHRmbep0wZaeq55hq5tW0LgYFOROsITfBFxIm4E3T6rBPHzh3j58E/065aO6dDKtQSEtL/jyMi4KqrpFzKM89AO33rvFNqwl+6ND3hHziQeZ+gIOm07dJFEn+nTl5dLM2pUTQzgOuBMOAY8JK19rPLPcfXEzzA7pO7ibsYR7NKzZwOpUiZPl3a5RMS5PtrrpHquLfc4tOf3n3DgQOwbFn6LeuQTGPk7N+5c3rCr1fPa8be6kSnIiw6Ltonq0/mx6FDMiN24sT0UXj168Ozz0ryVz4iOlpG6axYIaNz1qxJP/OnKl9eCiF16iRf27UrsmPxNcEXUQt3L2TgzIGM7z2eQc0GOR1OkREbK+vDvvee9M8NGQJTpzodlXLMhQuyitXKlZL0f/9dyphmZIyMzmnfPv3WtGmRGI+vCb6Iem/Ve/zjp39wX4v7+G///zodTpGTmCilzZs1g8aNZdvs2bJ84PDhMmNeyxT7IGvlzP/777Ls2KpVMhzr4sXM+5UoAa1by9V9u3bSeVuvXqFr89MEX0RZa1mwcwG9G/TWMfIu0r27DLcGWTP2kUfg/vtltTrlw86fl2Jpf/yRftu169L9SpWSpN+2rdxat4a6dR1N+prgvUR8YjwrD67kxto3Oh1KkRUVJVfw48fDvn2yLSBAOmNHjNCFwVUGJ05Itcw//pAmnrVrpaMnq9BQaNVKkn3r1nK/UaP06dhupgneCyQlJ3H7zNuZt2MeU/tP5d7m9zodUpGWlCSlUcaPhx9+kO/ffz+9wFlSkjbfqGwcOZKe7Neulaadw4cv3S8oSNrwW7ZMvzVvLicDF7tcgvfMKUYVmJ/xo02VNvy671eaV9KKWwXl7y/lTnr1kv/Z6dNhUIZ+7P/7P5lbM2gQ3HmnTqBSKapUgT595Jbq6FFJ9OvXy9d16+TjYepJIKO6dSXRN2smX5s3hzp13HY1oVfwRcyxs8eoVLKS02F4NWulU3bHDvnez0/WnB40SBYML13a2fhUEXD6NGzeLO36qbctWy7tyAVZCatpU/j0U0n8eaRNNF5q4e6FzN8xn3d6vEOgv+9MzfaEuDhZke7LL6UJJ/X/snhx+OADqVevVJ4kJMDOnZL4M95S2/UPHIDq1fP8stpE44ViLsRw16y7OBV/inbV2jGkxRCnQ/IqwcHSNHPnnXDypAyv/PJLmR1fL73CM4sWwcGD0K+fzJ1RKkeBgXKl3rRp5vbAkyfhzz/dsoatXsEXYasjVzN983Q+7PmhLv/nIZGR0h6fOkCid29Zpc7fX4ob9u0rzbN16jgapvIhjiz4odyvQ3gHPur1UVpyPx1/mrdXvk1ScpLDkXmv8PDMo99uv10WIQKpbPvkk9KP1qQJTJ7sTIxKpdIE70UenPcgzyx6hn/89A+nQ/EZf/sbLFwoi5JMnw533CFzYbZty7zm9J9/yvj77IZRK+Uu2gbvRR5r9xhbj29lRMcRTofic8qVk+UE771XOmSXLZMChqk+/xzeekvuN2kiZRK6d5eql8HBzsSsvJ+2wXuZpOQk/P3Sx9S+s/IdBjQeQJ2y2ijspJkz5Qp/yZLMV/bFiklH7vTpzsWmijZtg/chGZP7wt0L+eeif9L+k/bEXYxzMCo1cCDMmycDJpYuheefhzZtpCBaxiv4Q4fghhtkacLFi+GsrsOuCkCbaLxYi0otuKfZPTSr2IzgYtoOUBgEBsoyotddJwuJnzolY+5T/fKLnACWLpXv/f1llnuXLnLr00dmwSuVG9pE4wOstWkjbeZsn8P0zdMZc8MYmlbU1aoLm1OnJMkvXy639eulLg7I6J0zZ9Kv+GfNgkqV5JOAtuP7Lp3o5OMyjpF/Y8Ub/HHoD7rV6aYJvhAqW1bKIdx6q3x/7hysXi3JPjo6PZFbK7NpT52Sq/xmzSTRp96aN9crfaVX8D7nSOwRxq0Zx7+u/RfFA4oDsGTvEqqGVqVRWCOHo1O5de4cPPWUJP8//4Tk5MyPjxsni5qAzI4/cgRatPDqtad9ltaiUTm6kHiBBh81IDImkjUPraF1ldZOh6Ty6Nw5acpZv14KGa5bJ5OsOnSQx0eNgrFj5X6NGpLoW7RIL2jYSM/rRZo20agcnU88T696vfh1/6+0rNwybfsve3+hReUWlCtRzrngVK6EhMh4+muuyf7xKlVk8aEtW6Se1YEDUkgNJNFv3Cj3rZWO34YNJenXq6fNPEWdXsErAC4mXaSYvywwfC7hHFX/U5WEpAR2PbGLaqWqORydcoXERIiIgE2b5LZ1KzRoAG+/LY8fPChX+KmMkWUNGzWS28MPy/6qcNEreHVFqckdICouik7hnTibcDZTcv94zce0qtyKDuEddI3YIiggQOrcN24Md9116eP+/vDMM7B9u9TC370b9uyR2/ffS92dVP/6F8ydC/Xry5V+3bpyMqhTR04SgVq9ulDQK3iVowuJF9I6YqPORVHlnSoYYzj2z2NpTTcZh2Aq73LhgiT5v/6S2/Dh6Z20/frJxK3sdOwIv/8u9xMT4dVXoWbN9Fv16noCcCW9glf5kprcARKTE3myw5OcuXAmU7v81ZOvpmpoVcb3Hk/FkIpOhKncpHhxqaeTsaZOqmnTYNcuafLZtQv27k2/2q9bN32/yEh4+eXMzzVGSi5Xrw7vvQedOsn2HTtknetq1aTfQE8CBacJXuVKldAqvNPjnUzbDsUcYlXkKkIDQykbVDZt+/RN0/EzfvRu0JsyQWU8HKnyhNKl08fcZ5WxUSAoSEbx7N+ffjt0SIZtHjkiyyGmGjcOPvww/fsKFaBqVUn4rVrJJ4FUq1dDxYpyoihRwvU/n7fQBK/yrVqpauz5+x52nNiR1oZvrWXk4pEcjj3M1uFb0xL8mkNriE+Mp3WV1oQEhjgYtXK3jC12lSvDa69lfvziRUnuBw/K4kapqlWDdu3g8GFZxzoqSm6bNskM3lQXLkgzUKrQUEn2FSvKzN4RI6QUBMgnir17ISxMThhhYb71yUATvCqQ2mVrU7ts7bTvE5ISeLzd42w4uiHTxKmxy8cy5685fD7gc+5pfg8AO0/sZOvxrbSr1o7wUq5frkwVTsWKSUdsxhE7AM89JzeQ8gzHj0uyP3xYmotSnT0L7dvLSeDoUYiNldvu3fL44MHp+37zjRR2yyg0VJZXrFFDlmBM9dFHMmGsXDmZUVyuXPr9MmWK5olBE7xyqeIBxXn+mucv2d44rDF7K++lVZVWadtmbZvFC0teYESHEbx787sAHD17lPdXvU/rKq0Z2GSgx+JWhYu/v7TDV6lyaTNQ+fLSRAPSHHTmjJwMjh2Trxmv7qtWlaUUo6Pl00B0dPoJIevs3zFjZJ/svPSSVPgEOSmMGiVJv3TpS2/3359eUmLvXjlZlSolt+LFM3/CcTdN8MojXuv6Gq91zfxZvWbpmnSv253ONTqnbfvz2J+8seINrq15baYE3/6T9pQOKs2cO+dQMrAkAFuPbyXZJlO3XF2tlumjjJFEW6ZM9mP0hwyRW6rkZIiJkc7cuCwVtB97TBL8qVNS1vnkyfT7YWHp+x08CCtX5hzTPfek3x86VNYASBUQACVLyqeI226Dd9/Ny0+bd5rglWPuaX5PWnNNqpplajL6utGZxt/HXohlzeE1BAUEEVIsvf3+6YVP89Pun5h/93z6NOgDwNy/5vLZhs+4rfFt3N/yfgDiLsbxy95fqFSyEm2rZjuaTPkIP7/0E0JWL72Uu9e4+WZZsevUKfn0kPUWGpq+b6VKMj8gNlYeu3gRTp9Ov7mbWxO8MeZm4H3AH/jUWvuGO4+nir4G5Rvw0vWZ/9OCiwWz5dEtRMVFZRpzX7N0TRqHNaZ6qepp27Yc38KCnQtoVrFZ2ra9p/bSZ0YfGoc1Zttj29K2NxvfjOi4aFY9uIqaZWoC8N+N/+WXfb9wX4v7uLH2jQAcOHOAn3b9RJ2ydehap2va8zce3Uhx/+I0CmuUFldiciL+xl/nBnixsDCpzZ8bX36Z+fsLF6QPITY2c7+Cu7htOqIxxh8YB/QErgLuNsZkM6JWqcvz9/OnScUmXF/r+kzbJ94ykW2PbaNF5RZp2wY1G8TcO+cyqNmgTM/vWa8nXWpk/q88EnuEo2ePUqJY+ji75QeWM23TNHaf3J22bePRjQxbMIx3V6V/nrbW0mpiK676+Cos6eMC+3/VH78xfizYuSBt25d/fkmDDxsweunotG3RcdF0m9aNe7+9N1NM/1ryLx747gEOnDmQtu2nXT/x3KLnWLxncdq2w7GHeWP5G0zflHmtv6kbpzJp3aRMK3itilzFjD9nZPqZjsQe4bu/vmPNoTWZfqZFuxdlOg7ISXPZ/mWcOn8qbduxs8dYHbmafaf3pW1LSEpgw5ENbD2+NdPzI05EsOX4FuIT49O2RZ2L4q/ov4iOi07bdv7ieSJORHDwzMFMz993eh97Tu0hKTkp0/u37/Q+ziakL3kVdzGOA2cOEHUuvSHdWktkTCSHYjKvdh4dF82hmEOZYjqXcI4jsUc4E58+ZCcpOYmjZ49y/NzxS55/9OxRLiZdTNt2NuEsx84eyxTTxaSLHD93nJPnTwKS1MuXh5AKUfiXOk6yzdIR4GrWWrfcgE7ATxm+fx54/nLPadOmjVXKU85eOGsjz0TaxKTEtG2rDq6yUzZMsTuid6RtW3d4nX3wuwft+6veT9t2IfGCbTmhpW32cbNMr9ljeg/LaOxPu35K2/bBqg8so7GP/e+xtG17T+21jMbWeLdGpuc3/bipZTR209FNadtG/TzKMhr76q+vpm37/eDvltHYDp90yPT88m+Wt4zGHj97PG3b3+b+zTIa+9n6z9K2zd422zIa2/+r/mnbEhITLKOx/i/7Z/sz/RDxQ9q2D1d/aBmNHb5geNq2faf2WUZjq/+neqbnN/u4Wa5+ppUHVl72Z4o6F5W27f6591tGYyevn3zJzzTgqwGX/EwBYwKu+DN9tPqjHH+mrL+n3P5MOf2ewt4Ku+T3lF/AWptDTnVnE001IOOpOBLokHUnY8wwYBhAjazjppRyo5DAkEvG5HcI70CH8Mx/pq2rtObTvp9m2hboH8iGhzdc8po/3vsjSclJmZpohrQYQo96PQgNTG+crRhSkYX3LsxUAwjg1RteJTouOtOw0R71elCqeCmuqZleLrJKySo8e/WzVC9dPdPzh7QYwtmEs5k+lbSv1p7zieczLbxeuWRlbmlwC+2qtkvbZoyha+2umdb1BWhSoQlnE85mmsxWKaQS7aq2S2vaAqln1KJSCyqXrJzp+fXK1SPZJlPcP71NIiw4jIblG2aaFR0UEES9cvUu+ZlqlqlJqeKlMtU/Kl+iPDVL18z0+ysRUILqpaoTFhyW6fnVQqtd8jOVDy5P1dCqmWIKLhZM5ZKVKVW8VNo2P+NHpZBKVAiucMnzK4ZUJMAvPYWGBIZQMaRippgC/AKoEFyBsiXKZnp+WHAYBuP2pjy31aIxxgwEelhrh6Z8Pxhob619IqfnaC0apZTKm8vVonFnScBIIOOpOBw47MbjKaWUysCdCX4NUN8YU9sYEwjcBeRQf04ppZSrua0N3lqbaIx5HPgJGSY52Vq79QpPU0op5SJuHQdvrf0e+N6dx1BKKZU9XZZHKaW8lCZ4pZTyUprglVLKS2mCV0opL1WoFt02xkQB+/P59DAg+op7eZ7GlTcaV95oXHnjjXHVtNZWyO6BQpXgC8IYszan2VxO0rjyRuPKG40rb3wtLm2iUUopL6UJXimlvJQ3JfhJTgeQA40rbzSuvNG48san4vKaNnillFKZedMVvFJKqQw0wSullJcqUgneGHOzMWaHMWaXMWZkNo8bY8wHKY9vNsa09lBc1Y0xvxhjthtjthpjnsxmn+uNMWeMMRtTbi96KLZ9xpg/U455yWoqTrxnxpiGGd6HjcaYGGPMiCz7eOT9MsZMNsYcN8ZsybCtnDFmkTEmIuVr2Ryee9m/RzfE9W9jzF8pv6c5xpgyOTz3sr9zN8Q12hhzKMPvqlcOz/X0+/V1hpj2GWM25vBcd75f2eYGj/2N5bSWX2G7ISWHdwN1gEBgE3BVln16AT8ABugIrPZQbFWA1in3Q4Gd2cR2PbDAgfdtHxB2mccdec+y/F6PIpM1PP5+AdcCrYEtGba9BYxMuT8SeDM/f49uiKs7EJBy/83s4srN79wNcY0G/pmL37NH368sj78DvOjA+5VtbvDU31hRuoJvD+yy1u6x1iYAXwH9suzTD5hmxSqgjDGmirsDs9YesdauT7kfC2xH1qQtChx5zzLoCuy21uZ3BnOBWGt/A05m2dwPmJpyfyrQP5un5ubv0aVxWWsXWmsTU75dhayS5lE5vF+54fH3K5UxxgB3ADNcdbzcukxu8MjfWFFK8Nkt4p01ieZmH7cyxtQCWgGrs3m4kzFmkzHmB2NMEw+FZIGFxph1RhY4z8rp9+wucv7Hc+L9AqhkrT0C8g8KVMxmH6fftweQT17ZudLv3B0eT2k6mpxDc4OT79c1wDFrbUQOj3vk/cqSGzzyN1aUEnx2y49nHeOZm33cxhhTEpgNjLDWxmR5eD3SDNEC+BCY66GwOltrWwM9gceMMddmedyx98zIUo59gZnZPOzU+5VbTr5vLwCJwBc57HKl37mrjQfqAi2BI0hzSFZO/m/ezeWv3t3+fl0hN+T4tGy25ek9K0oJPjeLeDu20LcxphjyC/zCWvtt1settTHW2rMp978Hihljwtwdl7X2cMrX48Ac5GNfRk4ujt4TWG+tPZb1AaferxTHUpupUr4ez2YfR943Y8x9QB/gHpvSUJtVLn7nLmWtPWatTbLWJgOf5HA8p96vAOBW4Ouc9nH3+5VDbvDI31hRSvC5WcR7HjAkZWRIR+BM6scgd0pp4/sM2G6t/U8O+1RO2Q9jTHvkvT/h5rhCjDGhqfeRTrotWXZz5D1LkeOVlRPvVwbzgPtS7t8HfJfNPh5fVN4YczPwHNDXWhuXwz65+Z27Oq6MfTYDcjiex9+vFN2Av6y1kdk96O736zK5wTN/Y+7oOXbXDRnxsRPpWX4hZdsjwCMp9w0wLuXxP4G2HoqrC/LRaTOwMeXWK0tsjwNbkZ7wVcDVHoirTsrxNqUcuzC9Z8FIwi6dYZvH3y/kBHMEuIhcMT0IlAcWAxEpX8ul7FsV+P5yf49ujmsX0iab+jc2IWtcOf3O3RzX9JS/nc1IAqpSGN6vlO3/Tf2byrCvJ9+vnHKDR/7GtFSBUkp5qaLURKOUUioPNMErpZSX0gSvlFJeShO8Ukp5KU3wSinlpTTBK59ljCljjBmecr+qMWaW0zEp5Uo6TFL5rJTaIAustU2djkUpdwhwOgClHPQGUDelTngE0Nha29QYcz9S3c8faIrUVgkEBgMXgF7W2pPGmLrIJLEKQBzwkLX2L0//EErlRJtolC8biZQqbgk8k+WxpsAgpC7Ja0CctbYV8DswJGWfScAT1to2wD+Bjz0RtFK5pVfwSmXvFyv1u2ONMWeA+Snb/wSap1QHvBqYmVIyB6C458NUKmea4JXK3oUM95MzfJ+M/N/4AadTrv6VKpS0iUb5slhkGbU8s1LTe68xZiCkrW3bwpXBKVVQmuCVz7LWngBWpCzU/O98vMQ9wIPGmNRKhC5bgk4pV9Bhkkop5aX0Cl4ppbyUJnillPJSmuCVUspLaYJXSikvpQleKaW8lCZ4pZTyUprglVLKS/0/4LrrbyChkoAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# function that returns dy/dt\n",
    "def model(y,t,k):\n",
    "    dydt = -k * y\n",
    "    return dydt\n",
    "\n",
    "# initial condition\n",
    "y0 = 5\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,20)\n",
    "\n",
    "# solve ODEs\n",
    "k = 0.1\n",
    "y1 = odeint(model,y0,t,args=(k,))\n",
    "k = 0.2\n",
    "y2 = odeint(model,y0,t,args=(k,))\n",
    "k = 0.5\n",
    "y3 = odeint(model,y0,t,args=(k,))\n",
    "\n",
    "# plot results\n",
    "plt.plot(t,y1,'r-',linewidth=2,label='k=0.1')\n",
    "plt.plot(t,y2,'b--',linewidth=2,label='k=0.2')\n",
    "plt.plot(t,y3,'g:',linewidth=2,label='k=0.5')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('y(t)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "Find a numerical solution to the following differential equations with the associated initial conditions. Expand the requested time horizon until the solution reaches a steady state. Show a plot of the states ($x(t)$ and/or $y(t)$). Report the final value of each state as $t \\to \\infty$.\n",
    "\n",
    "## Problem 1\n",
    "\n",
    "$\\frac{dy(t)}{dt} = -y(t) + 1$\n",
    "\n",
    "$y(0) = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns dy/dt\n",
    "def model(y,t):\n",
    "    dydt = -y + 1.0\n",
    "    return dydt\n",
    "\n",
    "# initial condition\n",
    "y0 = 0\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,5)\n",
    "\n",
    "# solve ODE\n",
    "y = odeint(model,y0,t)\n",
    "\n",
    "# plot results\n",
    "plt.plot(t,y)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('y(t)')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "\n",
    "$5 \\; \\frac{dy(t)}{dt} = -y(t) + u(t)$\n",
    "\n",
    "$y(0) = 1$\n",
    "\n",
    "$u$ steps from $0$ to $2$ at $t=10$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns dy/dt\n",
    "def model(y,t):\n",
    "    # u steps from 0 to 2 at t=10\n",
    "    if t<10.0:\n",
    "        u = 0\n",
    "    else:\n",
    "        u = 2\n",
    "    dydt = (-y + u)/5.0\n",
    "    return dydt\n",
    "\n",
    "# initial condition\n",
    "y0 = 1\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,40,1000)\n",
    "\n",
    "# solve ODE\n",
    "y = odeint(model,y0,t)\n",
    "\n",
    "# plot results\n",
    "plt.plot(t,y,'r-',label='Output (y(t))')\n",
    "plt.plot([0,10,10,40],[0,0,2,2],'b-',label='Input (u(t))')\n",
    "plt.ylabel('values')\n",
    "plt.xlabel('time')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "\n",
    "Solve for $x(t)$ and $y(t)$ and show that the solutions are equivalent.\n",
    "\n",
    "$\\frac{dx(t)}{dt} = 3 \\; exp(-t)$\n",
    "\n",
    "$\\frac{dy(t)}{dt} = 3 - y(t)$\n",
    "\n",
    "$x(0) = 0$\n",
    "\n",
    "$y(0) = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns dz/dt\n",
    "def model(z,t):\n",
    "    dxdt = 3.0 * np.exp(-t)\n",
    "    dydt = -z[1] + 3\n",
    "    dzdt = [dxdt,dydt]\n",
    "    return dzdt\n",
    "\n",
    "# initial condition\n",
    "z0 = [0,0]\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,5)\n",
    "\n",
    "# solve ODE\n",
    "z = odeint(model,z0,t)\n",
    "\n",
    "# plot results\n",
    "plt.plot(t,z[:,0],'b-',label=r'$\\frac{dx}{dt}=3 \\; \\exp(-t)$')\n",
    "plt.plot(t,z[:,1],'r--',label=r'$\\frac{dy}{dt}=-y+3$')\n",
    "plt.ylabel('response')\n",
    "plt.xlabel('time')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4\n",
    "\n",
    "$2 \\; \\frac{dx(t)}{dt} = -x(t) + u(t)$\n",
    "\n",
    "$5 \\; \\frac{dy(t)}{dt} = -y(t) + x(t)$\n",
    "\n",
    "$u = 2 \\, S(t-5)$\n",
    "\n",
    "$x(0) = 0$\n",
    "\n",
    "$y(0) = 0$\n",
    "\n",
    "where $S(t-5)$ is a step function that changes from zero to one at $t=5$. When it is multiplied by two, it changes from zero to two at that same time, $t=5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns dz/dt\n",
    "def model(z,t,u):\n",
    "    x = z[0]\n",
    "    y = z[1]\n",
    "    dxdt = (-x + u)/2.0\n",
    "    dydt = (-y + x)/5.0\n",
    "    dzdt = [dxdt,dydt]\n",
    "    return dzdt\n",
    "\n",
    "# initial condition\n",
    "z0 = [0,0]\n",
    "\n",
    "# number of time points\n",
    "n = 401\n",
    "\n",
    "# time points\n",
    "t = np.linspace(0,40,n)\n",
    "\n",
    "# step input\n",
    "u = np.zeros(n)\n",
    "# change to 2.0 at time = 5.0\n",
    "u[51:] = 2.0\n",
    "\n",
    "# store solution\n",
    "x = np.empty_like(t)\n",
    "y = np.empty_like(t)\n",
    "# record initial conditions\n",
    "x[0] = z0[0]\n",
    "y[0] = z0[1]\n",
    "\n",
    "# solve ODE\n",
    "for i in range(1,n):\n",
    "    # span for next time step\n",
    "    tspan = [t[i-1],t[i]]\n",
    "    # solve for next step\n",
    "    z = odeint(model,z0,tspan,args=(u[i],))\n",
    "    # store solution for plotting\n",
    "    x[i] = z[1][0]\n",
    "    y[i] = z[1][1]\n",
    "    # next initial condition\n",
    "    z0 = z[1] \n",
    "\n",
    "# plot results\n",
    "plt.plot(t,u,'g:',label='u(t)')\n",
    "plt.plot(t,x,'b-',label='x(t)')\n",
    "plt.plot(t,y,'r--',label='y(t)')\n",
    "plt.ylabel('values')\n",
    "plt.xlabel('time')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}