{
 "metadata": {
  "name": "",
  "signature": "sha256:961d5e162ffa88888a5b06c22419b31141e2616fb9a01c7608a43a1f75bdb12c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Stiff Problems, Implicit Methods and Computational Cost\n",
      "\n",
      "By Thomas P Ogden (<t.p.ogden@durham.ac.uk>)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Stiff Problems\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y_0 = np.array([1.]) # Initial condition\n",
      "\n",
      "N = 10 # Num of steps to take\n",
      "t_max = 5.\n",
      "t = np.linspace(0., t_max, N+1) # Array of time steps\n",
      "\n",
      "solve_args = {}\n",
      "solve_args['a'] = -15."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from exp import exp\n",
      "from ode_int import ode_int_ee\n",
      "\n",
      "y = ode_int_ee(exp, y_0, t, solve_args) # Solve the ODE with Explicit Euler\n",
      "\n",
      "plt.plot(t, y[:,0], 'b-o', label='Explicit Euler')\n",
      "plt.plot(t, np.exp(solve_args['a']*t), 'k--', label='Known')\n",
      "plt.xlabel(r'$t$')\n",
      "plt.ylabel(r'$y(t)$')\n",
      "plt.legend(loc=2)\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAEqCAYAAABuj+WBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYVOX+NvB7DWfEEwcz8YBZKSa6K0lO0hAKmvm6011E\nItrOSivUX3BlVhbqdqtdah46/XZaykvyUpHtFAUNITxhmM0oHjBMRFEUUUpSkMPz/jENgZxmYGAt\nZu7Pda1rYNaznvWd0bx71nrWWpIQAkREROZOJXcBREREHYGBR0REFoGBR0REFoGBR0REFoGBR0RE\nFoGBR0REFoGBR0REFoGBR0REFoGBZ6D09PSgiIiI+H/84x9fX7x4sY/c9RARkXEYeAb65ptvJq9e\nvfq1KVOmJB0+fHik3PUQEZFxrOUuoLMYN25cyvjx43eWlZU5JScnT5C7HiIiMg5HeAASExPDFixY\nsEz/e2VlpU1ERES8j49Plr+///7c3NzBX3755dPZ2dneBw4c8Fu/fn2UnPUSEZHxLDrwhBDS2LFj\nd8+YMWOTJEm1d9GOi4uLdHNzK87KyvJZvnz5G9HR0asmTJiQPH369M2zZs365Omnn/5SzrqJiMh4\nkiU8LWHfvn0BAwYMONevX7/z+t8HDhx41t3dvbC6utoqLi4u8vTp0/cvW7ZsAQA8++yzW2bPnv3x\n6NGj9wohpH79+p2/cOFCX3k/BRERtYVFjPACAgL2JSUlTTl37tyAzMzMwJKSEhd3d/dCALCysqpW\nqVQ1dduXlJS4uLi4lACAJEmi7uiPiIg6J4sIPACYN2/emjVr1sw7ceLE0EmTJv23ubbOzs7XSktL\newC6w54MPCKizs9iAi8lJWVcaGhoqkqlqjl9+vT9zbUNDg5OS0pKmgIAqampoYGBgZkdUyUREbUX\ni7gsYe/evaNVKlVNSEjILgDYsGHDzC5duvyhP6wJ6A5d6n+ePn365sjIyDhvb+9sJyensvj4+Ag5\n6iYiItNR5KSVxMTEMI1G8zf9JJI71dTUqPz9/ffHxsbGhoaGpnZ0fURE1Pko6pBmU5cJ3Gn9+vVR\nubm5g3lujYiIDKWowJMkSaSkpIz76KOPXhZCSI21KSgo6J+SkjJu0qRJ/22qDRER0Z0Udw6vscsE\n6pozZ8661atXv/bee++93twIj6M/IiLLYcgASFEjvJbEx8dHeHl5HfP09DwJtPwBhRBcDFzeffdd\n2WvoTAu/L35f/L6UsxhKcSO85uzbty8gJydnWFBQUPqpU6eGHDly5KGuXbve8PPzOyB3bUREpGyK\nDTz9IcnCwkL3mJiYlQkJCeGffPLJLP3655577vPw8PAEhh0RERlCkZclmIIkScJcP1t7yMjIgFqt\nlruMToPfl3H4fRmH35dhkpMzsW7dLuzatRTCgHN4DDwiIup0kpMzMXduKs6cWQpAMijwOtWkFSIi\nIgBYt27Xn2FnOMWew2svzs7OuH79utxlUDvr2bMnrl27JncZRNROKiqMjy+LC7zr168bNY2VOidJ\n4j0JiMyZnV2V0dvwkCYREXU6c+aEoHfvt4zahoFHRESdzoQJgXjssVAACw3ehoFHRESd0q1bgbj/\n/iUGt2fgKZiHhwesra1hY2NTb7n//mafX9sitVqNDz/8EABw3333ISUlpdn2cXFx8PX1BQDk5+dD\npVLh5s2bTfbdWM22trYt1tVS30REdWm1wN/+Znh7Bl4dycmZCA19G2p1LEJD30ZycusfdG6KviRJ\nwnfffYfKysp6y+nTp1tdl75f/aSOX375BePGjWu2fWRkJA4ePFjvvaYm/kiShHXr1jWo+fbt222q\nuTlVVcafvCaizu3334FffwVGjDB8Gwben/QXMe7a9S/88EMsdu36F+bOTW1VUJmyr6bcvn0bXl5e\nmDlzJgCgvLwcgwcPxtKlS2tHSps2bUK/fv3Qq1cvzJ8/H9XV1Q368fDwQHJyMgDg/PnzeOKJJ9Cj\nRw8MHDgQGzduBABs2rQJ3t7euHnzJu69914Aumn/Z8+eNbru2NhYPPXUU7W/5+TkQKVq/K/h4cOH\n4e/vj65du2LkyJH44YcfatepVCqsXbsW7u7u2Lp1q9F1EFHndvSo7tWYEZ7FXZbQlMYuYjxzZimm\nTVuI4cMDjerr6NFduH69YV/r1y/EhAnG9dXUSMrW1habN2+Gr68vIiMjkZKSgu7du2PBggUoKCgA\nAOzevRunTp1CcXExQkND0bdvX0RFRdXrRz/aE0Jg0qRJmDBhArZu3YqffvoJY8eOxejRo2vbOjo6\n4syZMxg4cCBKS0vh6OhoVM36/Rni8uXLGDduHDZu3IjQ0FAkJSVh4sSJyM/Ph7OzMwAgOTkZWq0W\nrq6uBvVJROZDo9G9coTXCk1dxFhTY2V0XzU1jfdVXm5cX0IITJ48GQ4ODvWWdevWAQAeeughvP76\n65g2bRo++OADbN68ud5oafHixejSpQs8PDzwwgsvICkpqcl9HT58GL/++itiY2NhY2MDHx8fbNmy\nBVZW9Wtu6RpGIQT+53/+p0HN0dHRBm2vbxMXF4fHHnsMkyZNgr29PaZOnYqHHnoI33zzTW27d955\nh2FHZKG0WsDVFejTx/BtOML7U1MXMfr4VKOFOR0NhIZWYdeuhu/b2zc8pNgcSZKwdetWPP744022\nefXVV7F8+XKMHz8enp6e9dYNGDCg3s9FRUVN9pOfn48BAwbUC7iJEycCAPbu3WtUzWvWrMHLL79s\nUPuamsaf9Xv27Fls3boVDg4Ote8JITBmzJja3/UjPSKyPBqNbnRnzD0mOML705w5IRg0qP5FjIMG\nvYmoqLGy9tWSmJgYeHt7Y/fu3UhPT6+37ty5c7U/5+fnw93dvcl+evXqhcuXL9d7b9WqVTh27JhJ\n61WpVPVCrqkQdnd3x7Rp03Dr1q3a5ciRI5g9e7ZJ6yGizqeqCsjJMe78HcARXi39ubX16xeivNwK\n9vbViIoaZ/Q5N1P31dwhwO3bt+Pbb7/FiRMn8Pnnn+P555+vF1BLly7FBx98gF9//RXr16/H22+/\n3WSffn5+6NKlC1avXo3Zs2cjNTUVixYtwrRp0+q1s7bW/ZW5du1aq87hDRgwAJs2bUJ5eTmsra3x\n8ccfN2gjSRLCwsLg5+eH/fv3w9vbGykpKZg2bRo0Gg1cXFya7J+IzN/p00B5uXHn7wBA9kezt+Mj\n30VjmnpfiTw8PISVlZWwtraut3Tt2lWUlpYKd3d3sWrVKiGEELdv3xbDhg0Tr7zyisjPzxeSJIkl\nS5YINzc34ebmJt59993aftVqtfjwww9r95GcnCyEEOLEiRNi9OjRwsnJSXh5eYmdO3cKIYTYtGmT\n8Pb2FkIIUVNTI/z9/YWdnZ04d+5cg5rVanWjNdvY2IgrV66I27dvi8jISHHPPfeIkSNHimXLlgmV\nSiWEEOLs2bNCpVKJP/74QwghRGpqqhg+fLhwcHAQw4cPFzt27Kjdj0qlEsePH2/yu+tMf85EZJwt\nW4QAhDh6VPf7n/+9t5gLFvc8PP2MRHOWn5+Pe+65B2VlZU2OwsydJfw5E1mq+fOBNWuAsjLAxqb2\nv3c+D4+IiMyLVgsMHaoLO2Mw8MwUH49DROZKozF+wgqg0MBLTEwMW7BgwbI736+oqLALCwtLHDVq\n1CFfX9+Du3fvNv20RzPg4eGB6upqiz2cSUTmq6gIuHy5FRNWoLDAE0JIY8eO3T1jxoxNkiQ1OAGT\nkJAQ7urqevXQoUOjtm3bNvGVV175UI46iYhIHlqt7rU1IzxFXZYgSZJISUkZFxcXF3n69OkGjwTw\n8PDIf/jhh38CAHt7+/KysjKnjq+SiIjkog+81ozwFBV4AGBlZVWtUqkavf2GWq3OAICcnJxhL774\n4n9iYmJWNtdXbGxs3W2hVqtNVygREXU4rRbo1SsDa9dmGL2t4gKvJYsXL34nKSlpypo1a+YFBQWl\nN9e2buAREVHnp9EAo0apERurrn1v0aJFBm3bqQIvISEh/PDhwyOzs7O9bW1t2+8Ba0REpDi3bgG5\nucCUKa3bXlGTVurST1opLCx0Dw8PTwCAlJSUcWfPnh0YGhqaGhQUlP7YY4/tkbdKIiLqKMePA9XV\nrTt/B4B3WlEqlUqFnJwcDB06FABQVlaGJ554AiUlJdi9ezd69+4tc4XK1ln+nInIcBs3AjNnAnl5\nwKBBf71v6J1WOtUhTUv122+/Yfz48aipqUFmZiZ69uwpd0lERB1OowGcnICBA1u3vWIPaZLOtWvX\nEBwcDEdHR+zZs6c27GJjY/HMM89g+vTp6Nq1K9zd3fHVV1/Vbrdt2zaMGDEC3bp1g5+fHw4ePAgA\nmDFjBl566SUAQEVFBRwcHPDWW7pHGZWXl8PR0RHnzp3DjBkzEBUVhUmTJsHR0RGDBg0y6rl4RESm\nptXqDmeqWplcHOE1oqnLFzIyMkzS3lCXL19GeHg4AGDHjh2wtbWtt/7rr7/GF198gQ0bNmDFihWY\nO3cunnrqKWg0GkRERCApKQkBAQH4+uuv8cQTT+Do0aMYP348Fi9eDAD48ccfYW1tjczMTADAwYMH\n4eHhUfvg2M8//xzbtm1DYmIi5s6di/nz5+PAgQNt+kxERK0hhC7wIiJa3wdHeAoWFhaGHj164MSJ\nE7UjtLoCAwMRFhYGGxsbTJ48GUVFRaioqMCnn36K8PBwjBkzBvb29oiIiMCDDz6I7777DiEhITh9\n+jRKSkqwd+9evPTSS9BoNKioqEB6ejrGjRtX2//TTz+NoKAg2NvbY9KkSTh79mxHfnwiolr5+cDv\nv7d+wgrAEV6jjB2ZtXUk15SpU6fi/fffx6xZszB16lRoNBq4urrWrq/7IFT9fTOrqqqQn5+PMWPG\n1OurV69euHLlCnr27ImRI0di37592Lt3b+2oLSsrCxkZGbUPiZUkqUH/VVVV7fI5iYhaotHoXltz\nSzE9jvAU7IUXXgAAvP/+++jevTsiIyMN2u6uu+5CQUFBvfdOnTqFgX+e6R03bhz27NmDn376CT4+\nPggKCsLOnTtx7NgxPProo6b9EEREJqDV6s7dDRvW+j4YeJ2Ag4MDtmzZgj179uC9995rtq0kSQgP\nD8fmzZuxb98+3Lp1Cxs3bsT58+cxefJkALrAi4uLw9ChQ2Fvb4+goCB88skn8PX1hZ2dHQBwSj8R\nKYpGA9x/P9CWh8DwkKZC3fk8uxEjRmDFihV4/fXXERgYCEmSGrTR/z527FisWLEC//znP3H+/Hl4\neXnh22+/hZOT7l7b3t7esLKyqh3N+fn5oby8vN75u+b6JyLqaFotMGpU2/rghedklvjnTGQ+SkuB\nnj2BZcuAN95ouN7QC895SJOIiBTt6FHda1tmaAIMPCIiUri2PPS1LgYeEREpmkYDuLkBbb2FMAOP\niIgUTavVje7aOm/O4mZp9uzZk7MNLQBvsE1kHqqqgJwcICqq7X1ZXOBdu3ZN7hKIiMhAublARUXb\nz98BPKRJREQKpr+lWFtnaAIMPCIiUjCtFrC1BQYPbntfDDwiIlIsjUZ3/0wbm7b3pcjAS0xMDFuw\nYMGyO9+vrKy0iYiIiPfx8cny9/ffn5uba4LMJyIiJRJCF3imOJwJKCzwhBDS2LFjd8+YMWOTJEkN\n7gsVFxcX6ebmVpyVleWzfPnyN6Kjo1fJUScREbW/oiKguNg0E1YAhQWeJEkiJSVl3EcfffRyY/dF\nS0tLC548efI3ABAQELBPo9GY6GsgIiKl0d9hxVQjPMVdlmBlZVWtUqlqGltXUlLi4uLiUgLowrGx\nUWBdsbGxtT+r1Wqo1WoTVkpERO2pqRmaGRkZrXrwtuICrznOzs7XSktLewC6w5/GBB4REXUuWi0w\nYADQo0f99+8cwCxatMig/hR1SLMlwcHBaUlJSVMAIDU1NTQwMDBT7pqIiKh96G8pZiqKDTz96K2w\nsNA9PDw8AQCmT5+++eLFi328vb2zV6xYMX/FihXz5a2SiIjaw61burusmOr8HWCBD4AlIiLly84G\nHnkE+OYb4Mknm2/LB8ASEVGnZcpbiukx8IiISHG0WqBbN8DDw3R9MvCIiEhxNBpg+HBAZcKUYuAR\nEZGi1NQAR4+a9nAmwMAjIiKFOXsWuHHDtJckAAw8IiJSGFPfUkyPgUdERIqi0ejO3Q0bZtp+GXhE\nRKQoWq3uga8ODqbtl4FHRESKotGY/vwdwMAjIiIFuX4dKCgw/fk7gIFHREQKop+wwhEeERGZtfaa\noQkw8IiISEG0WuCuu4DevU3fNwOPiIgUQ6Npn9EdwMAjIiKFqKwEjh9vn/N3AAOPiIgU4tQp4PZt\njvCIiMjMteeEFYCBR0RECqHRAHZ2urustAdFBV5lZaVNREREvI+PT5a/v//+3Nzceh9bCCHNnj37\n40cfffSHUaNGHcrIyFDLVCoREZmYVqu7f6a1dfv0r6jAi4uLi3RzcyvOysryWb58+RvR0dGr6q7/\n/vvvx1y/fr3nDz/88OiWLVuenTt37lq5aiUiItMRov1uKaanqMBLS0sLnjx58jcAEBAQsE+j0dT7\n6NbW1lU3btzoKoSQrl275tytW7ff5amUiIhM6dIl4OrV9jt/BwDtNHBsnZKSEhcXF5cSAJAkSUiS\nJOqu9/PzO3Dp0qW7hwwZcurixYt91q5dO7e5/mJjY2t/VqvVUKvV7VA1ERG1lUajezVkhJeRkYGM\njAyj9yEJIVpu1UHCw8MToqKi1vv5+R0QQkgDBgw4V1BQ0F+/ftGiRe9WVVVZL1myZGFxcbHbI488\n8qNWqx3R2EhPkiShpM9GRERNW7YMePNNoLQU6N7duG0lSYIQQmqpnaIOaQYHB6clJSVNAYDU1NTQ\nwMDAzLrrb9++bevm5lYMAN27d//N3t6+/M5RIBERdT4aDTBwoPFhZwxFjfAqKyttIiMj4/Ly8u51\ncnIqi4+PjwCAmJiYlQkJCeHXr1/v+c9//vOz0tLSHhUVFXYvv/zyRxEREfGN9cURHhFR5zFkCODp\nCWzdavy2ho7wFBV4psTAIyLqHP74A+jaFXj3Xd1irE55SJOIiCxPTo7usoT2nKEJMPCIiEhm7fnQ\n17oYeEREJCuNRjdZZcCA9t0PA4+IiGSl1eoOZ0otnoVrGwYeERHJpqbmr8Brbww8IiKSza+/6mZp\nMvCIiMisGXNLsbZi4BERkWy0WsDKCnjggfbfFwOPiIhko9Ho7rJib9/++2LgERGRbDpqwgrAwCMi\nIplcuwacP98x5+8ABh4REclEf4cVjvCIiMis6WdoMvCIiMisabVA797AXXd1zP4YeEREJAuttuPO\n3wEMPCIiksHt28Dx4x13OBNg4BERkQxOnQIqKznCIyIiM9fRE1YABh4REclAq9XdXeW++zpun4oK\nvMrKSpuIiIh4Hx+fLH9///25ubmD72yzYsWK+Q8++ODPI0eOPJycnDxBjjqJiKhtNBrAywuwtu64\nfSoq8OLi4iLd3NyKs7KyfJYvX/5GdHT0qrrrs7Ozvb/66qunsrOzvXfs2PH4/PnzV8hVKxERtY4Q\nHXtLMT1FBV5aWlrw5MmTvwGAgICAfRqNpt7pzB07djw+ffr0zdbW1lW9evW6kpiYGCZPpURE1FqF\nhUBJScdOWAGADhxMtqykpMTFxcWlBAAkSRKSJIm66y9dunR3cXGx2+OPP77j5s2bjq+++uoHDzzw\nwPGm+ouNja39Wa1WQ61Wt1PlRERkqLbeUiwjIwMZGRlGb6eowHN2dr5WWlraAwCEENKdgde1a9cb\nZWVlTjt27Hi8tLS0x4gRI7QhISG7unXr9ntj/dUNPCIiUgb9DM3hw1u3/Z0DmEWLFhm0naIOaQYH\nB6clJSVNAYDU1NTQwMDAzLrrfX19D3bv3v03AHBwcLjl4OBwS6VS1chRKxERtY5WC9xzD9CtW8fu\nVxJCtNyqg1RWVtpERkbG5eXl3evk5FQWHx8fAQAxMTErExISwgHgtddeW63RaP52+/Zt26ioqPVh\nYWGJjfUlSZJQ0mcjIiKd++/XzdBMSjJNf5IkQQghtdjOXEOBgUdEpDxlZbqRXWws8M47punT0MBT\n1CFNIiIyb8eO6S5L6OgZmgADj4iIOlBHP/S1LoNnad66dcshISEh/NixY15VVVXWN2/edFSpVDVd\nu3a9MWrUqENPPfXUV5xAQkREzdFqgR49gP79O37fBp3D+/7778ecOHFi6IQJE5IHDRp0pu46IYR0\n9OjR4WlpacHBwcFpI0aM0LZbtUbgOTwiIuXx9QXs7IBWXEbXJJNNWikvL7e/cOFC33vvvTevpc6O\nHz/+QHMXgnckBh4RkbJUVwPduwPPPw+sXWu6fk02acXe3r68btgVFRX11v988+ZNx7ptlRJ2RESk\nPGfOAH/8Ic+EFcCISSv//ve/39y5c+f4bdu2TdS/d/z48QfS09OD2qc0IiIyJ3JOWAGMmLTy5JNP\nbk1PTw/auHHj8999993/6d27d9EjjzzyY2FhoXtQUFB6exZJRESdn0ajexzQ0KHy7N/gwPP09Dzp\n6el5cuDAgWfHjx+/s6ioqHd2drb3Qw89dKQ9CyQiIvOg1QJDhuge/CqHFietVFRU2N24caOrq6vr\n1ZY6Kygo6N+/f/8Ck1XXBpy0QkSkLH37AkFBwP/9v6bt12STVuzs7CqysrJ8tmzZ8uytW7ccGmtz\n/fr1nv/5z39ePHfu3IDWFEtERObt6lXdc/DkOn8HGHhI84knnth+6dKlu99///3/uXLlSq/y8nL7\nyspKGysrq2pHR8ebffv2vfDCCy98qn+SARERUV36CStyzdAEWnHz6OnTp2/u1avXFT8/vwO+vr4H\ne/fuXdROtbUJD2kSESnH6tVAdDRw5Qrg5mbavtv1aQmnTp0akpWV5XPw4EHfn3766eGnn376y5iY\nmJVKurUYA4+ISDkiI4G0NN1hTVNrt8DLysryEUJIvr6+BwHgq6++emrEiBHazMzMwJkzZ25oZb0m\nx8AjIlKOESMAd3dgxw7T921o4Bl8WYLe999/P8bGxqZyzZo18xwdHW/279+/wNXV9epdd911uXWl\nEhGROauoAE6cACZMkLcOowPv73//+7c3b950nD9//gr9exs2bJjZr1+/86YtjYiIzMHJk0BVlbwz\nNAE+8ZyIiNrZ5s3AjBm64BsyxPT9d8onnldWVtpERETE+/j4ZPn7++/Pzc0d3Fi7mpoala+v78HU\n1NTQjq6RiIiMo9EADg7AfffJW4eiAi8uLi7Szc2tOCsry2f58uVvREdHr2qs3fr166Nyc3MHS5LE\nIRwRkcJptYCXF2BlJW8digq8tLS04MmTJ38DAAEBAfs0Gk2DSxQLCgr6p6SkjJs0adJ/DRnCEhGR\nfITQjfDkvOBcz+hJK+2ppKTExcXFpQTQnYNrbAQ3Z86cdatXr37tvffee72lEV5sbGztz2q1Gmq1\n2sQVExFRcy5cAK5fN+2ElYyMDGS04pHpigo8Z2fna6WlpT0AQAgh3Rlo8fHxEV5eXsc8PT1P6ts0\n11/dwCMioo6n0eheTTnCu3MAs2jRIoO2U1TgBQcHpyUlJU3x8/M7kJqaGhoYGJhZd/2+ffsCcnJy\nhgUFBaWfOnVqyJEjRx7q2rXrDT8/vwNy1UxERE3T30PTy0veOgCFXZZQWVlpExkZGZeXl3evk5NT\nWXx8fAQAxMTErExISAiv2/a55577PDw8PCEkJGRXY33xsgQiIvn94x+60Pvll/bbR7veS7MzYOAR\nEcnvvvt05+++/rr99tEpr8MjIiLzceMGkJenjBmaAAOPiIjaybFjule5bymmx8AjIqJ20R4zNNuC\ngUdERO1CqwV69gT69pW7Eh0GHhERtQuNRnc4U1LIPbEYeEREZHLV1bpzeEo5nAkw8IiIqB3k5QG3\nbilnwgrAwCMionagtAkrAAOPiIjagVYLWFsDnp5yV/IXBh4REZmcRgMMHQrY2cldyV8YeEREZHJa\nrbLO3wEMPCIiMrHiYuDiRWWdvwMYeEREZGL6RwJxhEdERGZNP0OTgUdERGZNqwXc3QFXV7krqY+B\nR0REJqXRKO/8HcDAIyIiEyovB06dUt7hTICBR0REJnTiBFBVxREeERGZOaXO0AQUFniVlZU2ERER\n8T4+Pln+/v77c3NzB9ddX1FRYRcWFpY4atSoQ76+vgd37949Vq5aiYioIa0WcHQEBg2Su5KGFBV4\ncXFxkW5ubsVZWVk+y5cvfyM6OnpV3fUJCQnhrq6uVw8dOjRq27ZtE1955ZUP5aqViIga0miA4cMB\nKyu5K2lIUYGXlpYWPHny5G8AICAgYJ9Go6l3FNjDwyN/1qxZnwCAvb19eVlZmZMcdRIRUUNCKPOW\nYnrWchdQV0lJiYuLi0sJAEiSJCRJEnXXq9XqDADIyckZ9uKLL/4nJiZmZXP9xcbG1t0WarXa1CUT\nEdGfCgqA0tL2n7CSkZGBjIwMo7dTVOA5OztfKy0t7QEAQgjpzsADgMWLF7+TlJQ0Zc2aNfOCgoLS\nm+uvbuAREVH76qgJK3cOYBYtWmTQdooKvODg4LSkpKQpfn5+B1JTU0MDAwMz665PSEgIP3z48Mjs\n7GxvW1vb23LVSUREDWk0gCQBXl5yV9I4SYgGgyjZVFZW2kRGRsbl5eXd6+TkVBYfHx8BADExMSsT\nEhLCp0+fvvnIkSMPubq6XgV0hz337NnzWGN9SZIklPTZiIjM3ZQpwLFjwOnTHbtfSZIghJBabGeu\nocDAIyLqWIMGAQ8/DHz5Zcfu19DAU9QsTSIi6px+/x349VflztAEGHhERGQCR4/qXpV4SzE9Bh4R\nEbWZkm8ppsfAIyKiNtNoABcX3XPwlIqBR0REbaa/w4rU4tQR+TDwiIioTaqqdJcjKPlwJsDAIyKi\nNvrlF92DX5U8YQVg4BERURt1hgkrAAOPiIjaSKMBbGwAT0+5K2keA4+IiNpEqwWGDgVsbeWupHkM\nPCIiahONRvnn7wAGHhERtcHly0BRkfLP3wEMPCIiagP9hBWO8IiIyKx1lhmagMIeAEtEZMmSkzOx\nbt0uVFSzCl0SAAAQwElEQVRYw86uCnPmhGDChEC5y2qWRgP06wc4O8tdScsYeERECpCcnIm5c1Nx\n5szS2vfOnHkLABQdevpbinUGDDwiog526xZw7hxw9iyQn697jYvbhcuXl9Zrd+bMUqxfv1CxgVde\nDpw6BTz5pNyVGIaBR0RkYpWVQEHBX2FWN9jOntXNaqzLzg5QqRr/57i83Krd622t48eB6mqO8Fql\nsrLS5rnnnvs8Ly/vXisrq+rPPvvsn4MHD841dP2dQkPfVswxcKUem2ddrMsS62qr6mqgsPCvELsz\n2C5cAGpq/mpvZQX07w8MHAg8/rjudeBAwMND99q7NzB+fBV27Wq4r8rK6o75UK2g0eheO0vgQQih\nmGXDhg3Pz5s3730hBDIzM0dPmDBhuzHr6y4ABCDEoEFviu3bfxBy2r79BzFo0JsCELUL62JdrEve\n2kJC3hKPPvquCAl5q0FNNTVCXLokxIEDQnzxhRBLlwoxc6YQwcFCDBokhI2NqPe5JEmIvn2FCAgQ\nYto0Id55R4jPPhMiPV2I/HwhKisNq+nO78vKaoGwtf1BJCW1z/fQVlFRQnTpIkR1tbx16KKs5YyR\ndG2V4dlnn90ye/bsj0ePHr1XCCH169fv/IULF/oaur4uSZIE8CgAwNo6H126eAAA/va3jEb3rdGo\nG33fFO2PHn0b16/va/C+tXU+/P3zO7weffs//jiLqqqBd6zJQM+eCzF8+JIOr0evfl1/te/ZcyG8\nvHR1abWN9z9iROP9m6K9rq5zDd7v2XMhamr2NtpPY59Xkpr/fhp7ntjPPzfe/qGHMvDzz2/j+vV/\n3bFGXe/vvd6wYbp66v5zDQDHjzfsXwjg/vsbb5+Xp67XTm/AgIza3wsK3sbNm/q6/mpvZZUPJycP\nSBLwyCMZsLFB7WJrq3tNTVVDknTflUqF2p9nzqzfXr+sXNl4+7VrG7bfvz8Tc+dOQ0XFX3/3ra3P\nols3Z4wa9XPtKK28vP73aWMD2NvrFgcH3ev69RkYOFA3erOzq9Na3fD7BICMjIxG39e3v3btN1y4\ncB01NbrPsXZtHD78MBCHDgH/+hfw5pu6z9Xa/k3dPjBQ92ggW1t565EkCUKIFp/Ep6hDmiUlJS4u\nLi4lgC6wdKFl+PqG8gEA1dWlqKoqhbV1j/You0U1NY1/zXL/v0ZTfz9qauQ9Z9BUXUJYQfXnlaNN\nPWRS1cSVpaZp33Rddds39+d6Z3DcqbqJo1dNtb9927C/X/r67Oz++lkfCpIEWFvXb6d3992Nty8s\nbLyeBx74q01JiTVu3mzYRr++pgYoLdWd77pzKSnRra8/jgIWLmx8v00JCmrs3V0A6v+PXlXVQFy/\nno+iIt09ISdM+Otw48CBwKxZjf9dCQkxrp6WODt3h7Nz99rfp00LxFNPATNnAm+/rTtntnGjaffZ\nWkIAR48Czz4LnDjRsfsuLS1FaWkpYmNjjdvQkGFgRy3PPPNMwv79+/2EEKipqZH69etXYMz6ugv+\nPKQJCBEa+rZpxs2tFBLylmj4ny7rYl2syxg1NULcvi3EH38IUVoqRHGxEBcvCnHunBB5eUKcPCnE\n0aNC/PSTEFlZQuzdK8SePUKkpgqxfbsQW7cK8eWXQnh6vttoXY8++q5pvwATqqkRYtky3aFTb28h\nCgvlrkiIs2d139snn8hdieGHNGUPubrLp59+OvO1115bJYTAzp07x02dOjXemPX1PljtObwFsp8z\naPxcButiXaxLDkr9HwRDfPut7pxZnz5CZGfLW8vWrbrvLStL3jqEMDzwFHUOr7Ky0iYyMjIuLy/v\nXicnp7L4+PgIAIiJiVmZkJAQ3th6d3f3Rg+uSJIkQkPfRlTUWEXMCktOzsT69btRXm4Fe/tq1sW6\nWJeMNd15gfegQW9i7dpxstdmiKNHgYkTgeJiYNMm4Omn5alj0SLdcuMG0KWLPDXoGXoOT1GBZ0qS\nJAlz/WxE1DZKDGJjXLkCTJ4M7N8PvPMO8O67TZ+Pbi9PPgmcPKm78FxuDDwGHhGZsYoK3WSaTZuA\nf/xD99qRI6177gG8vYHExI7bZ1MMDTw+LYGIqBOyswM++wxYuRJISgJGj9Zd8N4RfvtNd5F9Z3gk\nUF0MPCKiTkqSgOhoYPt2IC9PN+I6dKj993v0qO6109xh5U8MPCKiTu7xx4GDBwFHR+DRR4Evvmjf\n/elvKcYRHhERdbgHHtCN7nx8gIgI3V1Z6t7P05S0WsDVVXdjgs6EgUdEZCZcXYFdu4AXXgCWLdPN\n5CwrM/1+NBrd4cym7mCkVAw8IiIzYmsL/O//AmvXAtu2Af7+umfvmUpVFZCT0/kOZwIMPCIisyNJ\nwJw5wM6durDz9tZds2cKubm6SyI624QVgIFHRGS2QkKArCyge3fdjbQ3bWp7n1qt7pUjPCIiUpQh\nQ3STWQIDgeeeA2Jimn4qhyE0Gt1h0yFDTFdjR2HgERGZOWdn3eHNV14BVq0CJk0Cfv+9dX1ptboZ\noTY2pq2xIzDwiIgsgI0N8MEHwEcfASkpgK8v8Ouvxvej1XbO83cAA4+IyKLMnq27dOHSJeCRR4Af\nfjB826Ii4PLlznn+DmDgERFZnMceA378EXBzA8aMAT791LDt9BNWOMIjIqJO4957dbcjCw4GXnwR\nmDtXd41dc/S3FGPgERFRp9Kjh+7G0/PmAevWARMmAKWlTbfXaoH+/YGePTuuRlNi4BERWTBra+D9\n93WHNdPTdffi/OWXxttqNJ33/B3AwCMiIgAzZwLffw9cvQqMGqX7ua5bt3R3WemshzMBBQVeZWWl\nTURERLyPj0+Wv7///tzc3MF3tqmoqLALCwtLHDVq1CFfX9+Du3fvHitHrURE5igwEMjOBvr0AcaN\n013CoJeTo3v6QmcOPGu5C9CLi4uLdHNzK46Pj4/Yu3fv6Ojo6FXbt29/om6bhISEcFdX16uJiYlh\nV69edfXz8ztw+vTp++WqmYjI3AwcCBw4AEydqrtQPScHCA3NxFtv7QJgjbVrq2BvH4IJEwLlLtVo\nigm8tLS04NmzZ38MAAEBAfvCw8MT7mzj4eGR//DDD/8EAPb29uVlZWVOHV0nEZG569YN+PZb3TP1\n3nsvE5s2peLWraUAgL17gYsX3wKAThd6igm8kpISFxcXlxIAkCRJSJIk7myjVqszACAnJ2fYiy++\n+J+YmJiVzfUZGxtbd1uo1WpTlkxEZLasrIAVK4AdO3YhJ2dpvXVnzizF+vULZQu8jIwMZGRkGL2d\nLIG3ZMmShV9++eXTdd8rLi52Ky0t7QEAQgipscADgMWLF7+TlJQ0Zc2aNfOCgoLSm9tP3cAjIiLj\nubg0HhPl5VYdXMlf7hzALFq0yKDtZAm8hQsXLlm4cOGSuu9t2LBhZlJS0hQ/P78DqampoYGBgZl3\nbpeQkBB++PDhkdnZ2d62tra3O65iIiLLZGfX+NXo9vZteOSCTBQzS3P69OmbL1682Mfb2zt7xYoV\n81esWDEfAAoLC9315/NSUlLGnT17dmBoaGhqUFBQ+mOPPbZH3qqJiMzbnDkhGDTorXrvDRr0JqKi\nOt8keUmIRo8cdnqSJAlz/WxERB0pOTkT69fvRnm5FeztqxEVNVZRE1YkSYIQQmqxnbmGAgOPiMgy\nGBp4ijmkSURE1J4YeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeERE\nZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeEREZBEYeAQAyMjI\nkLuEToXfl3H4fRmH31f7UEzgVVZW2kRERMT7+Phk+fv778/NzR3cVNuamhqVr6/vwdTU1NCOrNGc\n8T8w4/D7Mg6/L+Pw+2ofigm8uLi4SDc3t+KsrCyf5cuXvxEdHb2qqbbr16+Pys3NHSxJkujIGomI\nqPNSTOClpaUFT548+RsACAgI2KfRaP7WWLuCgoL+KSkp4yZNmvRfIYTUsVUSEVGnJYRQxBISEpJ6\n/Pjxofrf+/bte76xdpMmTfr2xIkTnjNmzPg8NTU1pKn+AAguXLhw4WIZiyE5Yw0ZLFmyZOGXX375\ndN33iouL3UpLS3tAV7nU2OHK+Pj4CC8vr2Oenp4n9e2a2gdHf0REVJf052hIdhs2bJh58uRJz1Wr\nVkWnpKSMi4+Pj4iPj4+o22bWrFmf5OTkDLOxsak8derUkF69el35+OOPZ/v5+R2Qq24iIuocFBN4\nlZWVNpGRkXF5eXn3Ojk5lcXHx0e4u7sXFhYWusfExKxMSEgIr9v+ueee+zw8PDwhJCRkl1w1ExFR\n56GYwCMiImpPipmlSURE1J7MKvCMuXid6ktMTAxbsGDBMrnrULKKigq7sLCwxFGjRh3y9fU9uHv3\n7rFy16RkN27c6Pr3v//9W7VaneHv77//559/flDumjoD3ljDcIGBgZlBQUHpQUFB6a+++uoHLW4g\n9+UIplw2bNjw/Lx5894XQiAzM3P0hAkTtstdk9KXmpoaacyYMbvt7e1vLViw4N9y16Pk5fPPP5/x\n8ssvfyiEQHFxset99913Wu6alLwsWrTonTVr1swVQmDPnj1BEydO/E7umjrDsmbNmrk9e/a81txl\nV1wEysrKuhj7b7xZjfAMvXid/iJJkkhJSRn30UcfvSx4KUezPDw88mfNmvUJANjb25eXlZU5yV2T\nko0ZM+b7Z5555v8BQElJiUvXrl1vyF2T0vHGGobLzc0dnJ+f7xEcHJwWEhKy6/DhwyNb2sasAq+k\npMTFxcWlBND9Q85bjxnGysqqWqVS1chdh9Kp1eoMLy+vYzk5OcNCQkJ2xcTErJS7JiXz8/M7cNdd\nd10eP378zoiIiPgpU6YkyV2T0s2ZM2fd6tWrXwN0/4bJXY+S2djYVEZFRa1PS0sLXrt27dywsLDE\nmpqaZjNNlgvP24uzs/O1li5eJ2qLxYsXv5OUlDRlzZo184KCgtLlrkfJLly40Pfuu+++tHPnzvHn\nzp0b4Ovre1B/BIYaMubGGgQMGzYsx8vL6xgAeHp6nnR1db1aVFTUu0+fPheb2sasRnjBwcFpSUlJ\nUwAgNTU1NDAwMFPumsh8JCQkhB8+fHhkdna2N8OuZXPmzFmnn3jh4OBwi4c0m7dv376A9PT0oKCg\noPSUlJRxr7/++nsHDhzwk7supVq2bNmC2NjYWAC4dOnS3b///nu3u++++1Jz25jVdXhNXbwud12d\nwebNm6fn5uYO/ve///2m3LUo1fTp0zcfOXLkIVdX16uA7pDTnj17HpO7LqU6efKk50svvfS/KpWq\nprq62mrx4sXv8H8UDMMba7Tsxo0bXadOnfrF9evXe6pUqpqlS5e+FRAQsK+5bcwq8IiIiJpiVoc0\niYiImsLAIyIii8DAIyIii8DAIyIii8DAIyIii8DAIyIii8DAIyIii8DAIyIii8DAIzITJ0+e9OSd\ncoiaxsAjMhPp6elBDz744M9y10GkVAw8IjOwc+fO8Rs3bnz+woULfYuKinrLXQ+REvFemkRmYuLE\nidu2bds2Ue46iJSKIzwiM1BUVNS7d+/eRXLXQaRkDDwiM5Cdne39yCOP/Jidne198+ZNR7nrIVIi\nBh6RGejTp8/FwsJC97KyMidHR8ebctdDpEQ8h0dERBaBIzwiIrIIDDwiIrIIDDwiIrIIDDwiIrII\nDDwiIrIIDDwiIrIIDDwiIrIIDDwiIrII/x/ztwzG3na+EQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x388f110>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Stability\n",
      "\n",
      "[[ Describe stability regions and why the exlicit Euler method is unstable for the above problem ]]"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##\u00a0Implicit Euler"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[[ Write an implicit Euler solver here and solve the above problem]]"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Computational Cost\n",
      "\n",
      "If the Runge-Kutta method is $\\mathcal{O}(h^4)$, why would we ever use a less accurate method like Explicit Euler? The answer of course is the RK method we wrote has more computation steps, so we expect it to be more expensive in computer time. If accuracy were the only concern we could pick a method truncating the Taylor expansion way up the series to some huge order, but it might take an age to run and we likely don't need that level of accuracy. We make a trade-off: a sufficient level of accuracy in our solution and a reasonable amount of time to find it.\n",
      "\n",
      "(**A Note on Notation**: It's confusing, but we use '[Big O Notation][big_o]' to quantify the computational complexity of a method, as well as the order of accuracy like we've already seem. I'll use the blackboard bold $\\mathbb{O}$ to distinguish computational cost from calligraphic $\\mathcal{O}$ for order of accuracy.)\n",
      "\n",
      "As both the Explicit Euler and Runge-Kutta methods loop once through $N$ steps, executing a sequence of instructions, we expect both to have a computational cost of order $\\mathbb{O}(N^1)$ i.e. that the problem will be solved in linear time. When $N$ doubles, so should the running time in both cases. But with its computation of the interval $k$-points, we expect the Runge-Kutta method to be slower. We can check this with Python's `timeit` module.\n",
      "\n",
      "[big_o]: http://en.wikipedia.org/wiki/Big_O_notation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from exp import exp\n",
      "import timeit\n",
      "\n",
      "y_0 = np.array([1.]) # Initial condition\n",
      "\n",
      "solve_args = {}\n",
      "solve_args['a'] = 1.\n",
      "\n",
      "t_max = 5.\n",
      "\n",
      "from ode_int import ode_int_ee, ode_int_rk, ode_int_ab\n",
      "\n",
      "def solve_exp_over_a(method, N, max_a):\n",
      "    num_a_steps = 10\n",
      "    t = np.linspace(0, 5., N+1)\n",
      "    if method == 'ee':\n",
      "        ode_int_exp = lambda a: ode_int_ee(exp, y_0, t, dict(solve_args, a=a))\n",
      "    elif method == 'rk':\n",
      "        ode_int_exp = lambda a: ode_int_rk(exp, y_0, t, dict(solve_args, a=a))\n",
      "    elif method == 'ab':\n",
      "        ode_int_exp = lambda a: ode_int_ab(exp, y_0, t, dict(solve_args, a=a))\n",
      "#     elif method == 'ie':\n",
      "#         ode_int_exp = lambda a: ode_int_ie(exp, y_0, t, dict(solve_args, a=a))\n",
      "    for a in np.linspace(1., max_a, num_a_steps+1):\n",
      "        ode_int_exp(a)               \n",
      "        \n",
      "N_list = [2, 20, 200, 2000]        \n",
      "\n",
      "ee_time = np.zeros(len(N_list))\n",
      "rk_time = np.zeros(len(N_list))\n",
      "ab_time = np.zeros(len(N_list))\n",
      "ie_time = np.zeros(len(N_list))\n",
      "\n",
      "for i, N in enumerate(N_list):\n",
      "\n",
      "    print 'N: %d'%N  \n",
      "  \n",
      "    ee_time[i] = min(timeit.Timer('solve_exp_over_a(\"ee\", %d, 5)'%N,\n",
      "                    setup='from __main__ import solve_exp_over_a').repeat(7, 1))    \n",
      "    rk_time[i] = min(timeit.Timer('solve_exp_over_a(\"rk\", %d, 5)'%N,\n",
      "                    setup='from __main__ import solve_exp_over_a').repeat(7, 1))\n",
      "    ab_time[i] = min(timeit.Timer('solve_exp_over_a(\"ab\", %d, 5)'%N,\n",
      "                    setup='from __main__ import solve_exp_over_a').repeat(7, 1))\n",
      "#     ie_time[i] = min(timeit.Timer('solve_exp_over_a(\"ie\", %d, 5)'%N,\n",
      "#                     setup='from __main__ import solve_exp_over_a').repeat(7, 1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "N: 2\n",
        "N: 20\n",
        "N: 200"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "N: 2000"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(N_list, ee_time, label='Explicit Euler')    \n",
      "plt.plot(N_list, rk_time, label='Runge-Kutta')\n",
      "plt.plot(N_list, ab_time, label='Adams-Bashforth')\n",
      "# plt.plot(N_list, ie_time, label='Implicit Euler')\n",
      "plt.legend(loc=2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "<matplotlib.legend.Legend at 0x3876c10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEQCAYAAADh3jDlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYlFX/BvB7WBRUBEFQxD2T3FITlUVwEFe01zJNKUUt\n7U1xQcXcck8TU99KsczdKLMyy1zAFQVXXNNcUXFfkEVkFZjz++P5DcywowOz3Z/rmivkGYbDw+jd\n+T7n+R6ZEAJERES6ykTbAyAiIioOg4qIiHQag4qIiHQag4qIiHQag4qIiHQag4qIiHQag4qIiHQa\ng4qIiHRaqYNqy5YtA6dNm/Zl/s9nZWWZDx48ONTV1fW4h4fHkatXrzprdohERGTMSgwqIYSsW7du\ne4cNG7ZBJpMVaGOxadMmf3t7+7jjx4+7Llq0aOqkSZOWls9QiYjIGJUYVDKZTISFhfVcuXLlaCGE\nLP/x/fv3+/Tr1+8PAOjUqVPUuXPn2pTHQImIyDiZleZJpqamOSYmJorCjsXHx9vZ2dnFA1KoFTbr\nUiruGBERGZbCJjcv45UXU9ja2iYkJSXZANKgSgojIQQfWnjMnj1b62Mw1gfPPc+9MT406ZWDysfH\nZ//WrVvfA4Dw8PAeXl5eh199WERERJIyBZVytnT//n0nPz+/zQAwdOjQjQ8ePKjTvn376ODg4CnB\nwcFTymOgRERknGSanqIV+81kMlGR34/yREREQC6Xa3sYRonnXnt47rVHJpNBaOgaFYOKiIg0TpNB\nxc4URESk00q1PL282draIjExUdvDIANQo0YNJCQkaHsYRKRBOhFUiYmJGl/OSMZJJtNIpYGIXsGj\nlEcafT2W/oiISCOyFdn45vg3cF6h2ZavDCoiInplkbcj8daqtxAYHgi3um4afW0GFRERvbRHKY/g\nv80fXhu88CzzGf54/w/s/nC3Rr+HTlyjIiIi/ZKtyEbIyRDMipiFjOwMzPCcgeme01HFvIrGvxdn\nVGXUsGFDmJmZwdzcXO3RtGnTV3pduVyOkJAQAMDrr7+OsLCwYp+/adMmuLlJ0+vY2FiYmJggLS2t\nyNcubMyVKlUqcVwlvTYRGZ/8Zb4Loy7giy5flEtIAQyqMpPJZNi+fTuysrLUHteuXXvl11WuWLt+\n/Tp69uxZ7PP9/f1x7Ngxtc8VtXJSJpPh22+/LTDmFy9evNKYi5OdnV1ur01E2lFUma+p3av9j3pJ\n9CKoAgMBuVxzj8DA8hnnixcv0KpVK4wYMQIAkJGRAWdnZyxYsCB3ZrJhwwbUq1cPDg4OmDJlCnJy\ncgq8TsOGDbFz504AwN27d9GnTx/Y2NigUaNGWLt2LQBgw4YNaN++PdLS0tCkSRMA0j1Et27dKvO4\n58yZgwEDBuT++eLFizAxKfytcerUKXh4eMDKygouLi44dOhQ7jETExN88803cHJywrZt28o8DiLS\nTdmKbHx74ls4r3DGln+3YLrndFwafQnvNnu3Qm4J4TWql1DUzKVSpUrYuHEj3Nzc4O/vj7CwMFhb\nW2PatGm4c+cOAGDv3r24cuUK4uLi0KNHD9StWxdjx45Vex3l7EoIgb59+6J3797Ytm0bTp8+jW7d\nusHT0zP3uVWqVMGNGzfQqFEjJCUloUqVwqfexd2nVto32uPHj9GzZ0+sXbsWPXr0wNatW/H2228j\nNjYWtra2AICdO3fi/PnzqFmzZqlek4h0W9SdKATsCsA/j/9B99e6Y3mv5eU+gyqggvcnEYUp6vO6\nqEGDBqJSpUrCwsJC7fHNN9/kPufzzz8X9evXF1ZWVuLSpUtCCCFu3bolZDKZiImJyX3eV199JTp3\n7iyEEEIul4uQkBAhhBANGzYUO3fuFCdPnhTW1tYiOzs792u2b98uYmJixPr164WLi4vaa6emphY6\n5s6dOwtzc/MCY544caIQQojZs2eL/v375z7/woUL4v/7Mua+dkpKili8eLEYMGBAgddevXq1EEII\nmUwmIiMjy35SNUif3ktEuuzR80fCf5u/wByIesvqia2XtgqFQlHqr///v4sayQ7OqMpIJpNh27Zt\n8PX1LfI5Y8aMwaJFi9CrVy80a9ZM7ViDBg3UPn70qOg7uGNjY9GgQQOYmprmfu7tt98GAERGRpZp\nzF9//TVGjx5dqucrFIVu5oxbt25h27ZtsLS0zP2cEAJdu3bN/bNyZkVE+ilbkY2V0Ssx8+BMpGel\nY7rndEzvNB1VK1XV2pj04hqVvgkKCkL79u2xd+9eHDx4UO3Y7du3cz+OjY2Fk5NTka/j4OCAx48f\nq31u6dKluHDhgkbHa2JiohZORYWnk5MThgwZgvT09NzHmTNnMGrUKI2Oh4i0I+pOFNr90A7jw8bD\nta4rLo6+iAVdFmg1pAAG1UsRxVzv2bFjB/78809s2bIF06ZNw8cff4zU1NTc4wsWLEBaWhouXryI\n5cuXY+DAgUW+pru7O6pWrYply5YhPT0df/75J+bOnYtatWqpPc/MTJoYF9eMtbgxN2jQAGfPnkVG\nRgays7Px3XffFXiOTCbDwIEDsWPHDhw5cgQvXrzA9u3b4ebmhuTk5CJfm4h03+OUxxj25zB4rvdE\nYnoifh/wO8I+DKv4a1FFYFC9hL59+xa4J6l69ep49uwZPv30U8ydOxf16tXDtGnTULVqVUyZMiV3\nwULjxo3RsGFDdOnSBR999BE++eQTAIUvaDA3N88NPgcHB8yaNQu//vorHBwc1JazOzk5wd3dHU2a\nNMldtJHf+PHjC72PKi4uDh988AE8PT3RokULuLm5oWPHjmrjUX7cpEkThIaGYvTo0bCxscHMmTPx\nyy+/oFGjRkX+DESku7IV2Vh+YjmcVzjj5ws/Y1qnabgccBnvNX9Pp/4+68TGicoVboYsNjYWjRs3\nRkpKSpEr8+jVGcN7iUgTjtw5gtG7RuOfx/+gW+NuWN5rOZxraq6ZrCY3TuRiCiIiI/I45TGm7JuC\njec3om71uvh9wO/o16yfTs2g8mNQVSBdfiMQkWHLVmTju+jvMPPgTKRlpWFqp6n43PNzrS+UKA2W\n/sig8L1EVNCRO0cQsCsA5x+fL5cyX2FY+iMiohLlL/P9NuA3vNdMtxZKlAaDiojIwGQrsvH9qe/x\n+YHP9a7MVxgGFRGRATl69ygCdgXg3KNz6Nq4K1b0WlHuZb7yxqAiIjIAT1KfYMq+KdhwboNel/kK\nw6AiItJjhZX5ZnjOQLVK1bQ9NI1hZ4oyyL+7b+XKldGwYUPMmzdP20PLpdynStWRI0dgbW2NoKCg\nMn/9sGHDMHnyZI2Pk4he3dG7R9F+dXuM3T0W7Z3a459R/+BLny8NKqQABlWZyGQyrFu3LneH3MzM\nTPz0008IDg7G1q1btT28Qh08eBA9e/bEZ599hiVLlmh7OESkAU9Sn+Cjvz6CxzoPPE17il/7/4o9\ng/fgjZpvaHto5YJB9Yo8PDzQrFkz3LhxQ21nXkDqoj58+HAA0i66gwYNwtChQ2FlZQUnJyf89ttv\nuc/duHEjGjdujBo1amDGjBlo2LBhbt++mJgY9OjRAzY2NmjevHmpQzEsLAxvv/02Fi1ahBkzZuR+\n3sTEBJcuXcr9c//+/TF37lzs3r0bI0aMwOnTp9GoUSNMnToVP/74I5YtW5b7c2zbtg3NmzdHlSpV\n8PrrrxfawJaIykeOIgchJ0PgvMIZP/7zI6Z4TMHlgMsY0GKAQVyLKopeXKMKDAvEuUfnNPZ6bWq3\nwdc9v36pr1W9mTQnJwcHDhzAtWvX4O3tje+++65AM1fVP//+++/46aefsGbNGgQHB2P8+PEYMGAA\njh07hqCgIOzduxdNmzbF9OnTc0MqIyMDXbt2xdSpU/HXX3/h8OHDGDBgAJydndGyZcsix/nXX3/h\n/fffx9y5cxEQEFDsz6QcZ69evbB27VqsWLEC0dHRAKRdfe3t7bF48WJkZmbigw8+wI4dO+Dt7Y3I\nyEh4e3ujf//+sLe3f6nzSUSlc+zuMQTsCsDZR2fRtXFXLO+13GBnUPlxRlUGQgiMHDkSlpaWsLS0\nhIWFBXr06AE/P78C14WUz1cNNi8vLwwcOBDm5ubo168fHj16hMzMTGzcuBH+/v5o06YNqlSpguDg\nYJibmwMA/v77bzg4OODTTz+FhYUFunfvjnfeeQe//PJLkeOMiYmBv78/2rZti59//hkZGRll+hmL\n+pyJiQmOHz8OHx8fJCQkICcnByYmJoiLiyv16xNR2SjLfO7r3PEk9YnBl/kKoxczqped/WiaTCbD\nmjVr4O/vn/u5iIgI+Pj4YMyYMQWer1Ao1GZUdnZ2uR8rO6hnZ2fj9u3b6NmzZ+6xypUro2bNmhBC\n4NatWzhz5kyBXXUHDx6Mw4cPw8fHJ3dsL168ACDN9Hbu3IkWLVrgzTffxLhx4/DDDz8U+XMVtaNv\nYVatWoWwsDA4Ojqibdu2ZfpaIiq9HEUOVp1ehRkHZiDlRQo+8/gMM71mGtxCidLQi6DSZXK5HHXq\n1MHdu3cL3SnXwsKixNeoWrWq2k6+qampePr0KWQyGZycnODl5YUDBw7kHr9x4wYsLCzg5OSErKys\nAq/n7OwMT09PAMCmTZvQvXt3+Pj45G7SKJPJ1Mb58OFDtG7dusjxKcP2xx9/xKFDh3D58mVUrlwZ\niYmJWLlyZYk/HxGVjWqZz6eRD5b3Wo5m9s20PSytYelPA8zMzKBQKFC/fn0cOXIEAHDz5k3s2rWr\nxK+VyWTw9fXFxo0bcenSJaSmpiIoKCg3gHx9fXH58mVs374dmZmZiIqKgqurK2JjY0s1ti5duiAo\nKAiffPIJbty4AQBq4zxx4kTu9Sjlz5KSkoKcnJzcPycmJgKQZn8KhQKpqamIi4tDYGAgTExMkJaW\nVroTRUTFikuNw8fbP84t823pvwV7h+w16pACOKPSCCsrK+zbtw/ffvsthg0bhu3bt8PJyQkDBgzI\n/Qc//8IK5ecA4KOPPsL9+/fh4+OD9PR0fPrpp7CxsYGZmRlq1KiBv//+G+PGjYOfnx/q1q2LJUuW\nwMPDo9CxFPZ9vvjiCxw4cACDBg3CkSNHsHbtWowdOxYhISFo3Lgx+vbtm/tcDw8PpKSkoFmzZrh2\n7Rp8fX3h7+8PmUyGb7/9FmFhYahfvz7q1auH2bNnIysrC++88w7u3bunyVNKZFRY5iset/nQAU+e\nPEFWVhacnJwASKU/W1tbJCcno3LlyloenX4x9vcS6Z/j944jYFcAzjw8gy6NumBFrxUGMYPS5DYf\nLP3pgD179sDLywt37txBYmIipk2bBrlczpAiMmDKMp/bWjc8SnmELf23YN+QfQYRUprG0p8O+PDD\nD3HmzBm4uLggLS0Ncrkc69at0/awiKgc5Chy8MPpHzD9wHSkvEjBZPfJmOk1E1aVrbQ9NJ3F0h8Z\nFL6XSJflL/Mt77Ucze2ba3tY5YI7/BIR6ZG41DhM2z8Na8+uRR2rOvjlvV/wfov3DbrtkSaVeI0q\nKyvLfPDgwaGurq7HPTw8jly9elVtBy4hhGzUqFHfde7c+VDHjh1PREREyMtttEREeiRHkYPvor+D\n8wpnbDy/EZPdJ+NKwBUMbDmQIVUGJc6oNm3a5G9vbx8XGho6ODIy0nPSpElLd+zY0Ud5fN++fV0T\nExNrHDp0qPONGzde69ev3x/nz58v+u5RIiIjcOLeCYzeNRpnHp6Bd0NvrPBdYbBlvvJW4oxq//79\nPv369fsDADp16hR17ty5NqrHzczMsp8/f24lhJAlJCTYVq9ePbm8BktEpOviUuMwYvsIuK51xaOU\nR/jlvV+w338/Q+oVlDijio+Pt7Ozs4sHpMUQMplM7Uq1u7v70YcPHzq+8cYbVx48eFDnm2++GV/c\n682ZMyf3Y7lcDrlc/lIDJyLSJTmKHKw+sxrT90/H8xfPEeQehFles4xmNV9ERAQiIiLK5bVLXPXn\n5+e3eezYscvd3d2PCiFkDRo0uH3nzp36yuNz586dnZ2dbTZ//vyZcXFx9h06dDh5/vz51oXNrLjq\nj8ob30ukDSzzFVShN/z6+Pjs37p163sAEB4e3sPLy+uw6vEXL15Usre3jwMAa2vrZxYWFhn5Z12G\naN68eTAxMcFff/1V5HNSUlJgYqK9e6pjY2NhYmICc3NzmJubw8zMDPXq1cPEiRORnZ39Sq+df5NI\nVVeuXIGzszNq1ar10mNW9g8cNmwYJk+e/EpjJSovT9OeYuTfI+G61hUPnz/E5vc2s8xXDkr8V3To\n0KEbHzx4UKd9+/bRwcHBU4KDg6fcv3/fyc/PbzMABAUFLTl48KC3t7f3QblcHjFjxowFVlZWz8t/\n6NojhMD69evRrl07bNiwQdvDKdGzZ8+QlZWF7OxsbN++HT///DNCQkJe6TUL6ymodOLECchkMrWO\n8KWhGp7KWVFx34dIW3IUOfj+1PdourwpNpzbgCD3IFwdcxWDWg7i+7U8KDf3q4iH9O0KKurzuio8\nPFw0btxYnD17VlSqVEnExcUJIYTIysoSo0ePFtWrVxf29vZi7ty54v/LnUIIIb7//nvRuHFjYWlp\nKVq2bCl+//333GMymUwsXbpUvPbaa8LKykrMmjVLLF68WNjb24s6deqIlStX5j53/vz5ok6dOsLG\nxkb07t1b3Lt3r9Bx3rp1S8hkMpGamqr2+XfffVeMGTNGCCFEcnKy8PPzEzY2NsLa2lrt9VJSUoSf\nn5+oUaOGcHR0FGPGjBGZmZlCCCEaNmwoPv/8c9GqVSthYWEh5HK5ePLkifjrr7+EqampkMlkwt7e\nXgghxOXLl4WPj4+wtrYWzs7OYu3atblj6dy5s5g6dapo06aNGD58uDAzMxMymUyYm5uLGzduiGHD\nhokhQ4aInj17CktLS9G0aVNx7NixIn83+vZeIv1z4t4J0W5VO4E5EPINcnHx8UVtD0kn/f/fRc1k\nh6ZeqFTf7GWDavx4ITp31txj/Pjiv18JBgwYIBYuXCiEEKJNmzbif//7nxBCiCVLlohWrVqJ69ev\ni4SEBNG3b19hYmIihJBCo0qVKuLChQsiJydHbNq0SVSrVk0oFAohhBRU3t7e4unTpyIiIkLIZDLR\nr18/kZycLCIiIoSZmZlITk4Whw4dErVr1xYPHz4UaWlp4uOPPxb9+/cvdJzKoEpJSRFCCKFQKER0\ndLRwcHAQu3fvFkIIMWfOHNGrVy+RnJwsnj17Jv7zn/+IUaNGCSGEmDt3rujZs6dIT08XDx8+FG3b\nthUrVqwQQgjRoEED8cYbb4irV6+KpKQk4ebmJmbMmCGEEGLDhg3CxcVFCCFEamqqcHJyEsuWLRNp\naWnixIkTwtHRUWzfvl0IIQVV/fr1xYULF4QQQsTGxgqZTCbS0tKEEEIMHTpUWFhYiL1794r09HQx\nfvx44eHhUeTvhkFF5SUuNU6M3D5SyObIhOMSR7H5wubcv79UkCaDik1py+jp06fYtWsXhg0bBgAY\nOnRobvlvzZo1mD17Npo0aYIaNWpgyZIluSUse3t7nD59Gi1atMCTJ09gYmKC1NRUpKen57721KlT\nYWdnl7vp4aRJk2BlZYXOnTsjJycHDx48gEwmw7Nnz7Bnzx6kp6dj5cqVWLVqVbFjrlmzJiwtLVG5\ncmV06NABTZo0gYuLCwDpGtDGjRtRuXJlPH78GBYWFnj69CkAwNTUFDdu3MCRI0dgY2ODw4cPY+jQ\noQCkklxgYCCaNm0Ka2tr+Pj44ObNmwDUt7PfuXMnqlWrhgkTJsDS0hIdOnTA6NGjERoamvs6I0eO\nRMuWLdW+VvU13nnnHXTt2hUWFhbo06dP7vchqgg5ihysOrUKziucsf7cekx0m8gyXwXTjxZKX+vG\nVvSAtGNuRkYG3nzzTQDSdZXk5GScPXsWd+7cQePGjXOfW7du3dyPs7OzMW/ePBw7dgz169fHa6+9\nVuC1ldt8KBdg2NjYqB3PycmBp6cnfvzxR6xevRrjxo1DkyZNMHPmTPTt2xdmZma5f3EOHDiAevXq\nAQDi4+NRpUoVANKWIiNGjMCgQYOwb98+3LlzB4GBgUhLS8Prr7+O5ORk2NvbA5CC09zcHDNmzMCl\nS5fg4eGBpUuXonnz5gV+PlNT00J3G7516xaaNm2q9jl7e3u161e2trZFnm/lLsdKZmZmhX4fovJw\n8v5JBOwKwKkHpyBvKMeKXivQwqGFtodldDijKqO1a9ciJCQE58+fx/nz53Hx4kX4+vpi/fr1qFOn\nDm7dupX7XNVdeJcuXYq4uDjcvHkThw4dwtixYwu8dmn+7+zChQto1KgRwsLCkJCQgJEjR2Lw4MEA\npDDMyspCVlZW7qwsPwcHB4wcORInTpwAIHVuDwwMzN1FuHXr1rnjOHjwID744AMcP34cDx48gKOj\nI4KCgso03lq1auHu3btqn7t69apaoJeE/9dKFe1p2lN88vcncF3jivvJ9/Fzv59xwP8AQ0pLGFRl\ncPToUdy+fRtDhgxBnTp1UKdOHTg5OWHgwIHYvHkzhgwZgi+++AKxsbFITEzE1KlTc/+Rzc7ORnZ2\nNtLS0nDnzh3MmDEDAMq8jfuZM2fw4YcfIiYmBjk5OUhLS8udARVFtYx27949hISEoHPnzrnjSklJ\nQXZ2NsLDw7FlyxZkZGQAAH7++WdMmDABiYmJUCgUSEtLQ+3atcs0Xl9fX8TGxuKHH35Aeno6jh07\nho0bN2LEiBGFjs/MTJrkJyQkFDhGVN5Uy3zrzq7DRLeJuDLmCvxa+fF/mLSIQVUGa9euxX/+85/c\nMppSnz59crdvl8vlcHFxQatWrdCnT5/czQ8nTJgAmUwGe3t7vPvuu5g4cSLc3NzUtoFXVdS29UOG\nDEH37t3h6emJGjVq4LfffsOvv/5a7LhtbGxgbm6OSpUqwd3dHXXr1sVPP/0EAFi1ahUWLVoEW1tb\nbNq0CRs2bEBERARWr16NhQsXIjMzE40bN0bDhg2RlZWF4ODgIserHKPqx/b29ti5cyfWrFkDOzs7\nDB48GIsWLYK7u3uhP6uTkxPc3d3RpEkT3L59u9Dl6fwHg8pD9P1ouK51xac7P0Urh1Y49+k5LOm+\nBNUrV9f20Iwe96Mig8L3EpXV07SnmL5/OtacWYPa1WpjSfcl8GvJGdSr4n5URESvKEeRg7Vn12La\n/ml4lvEME9wmYHbn2ZxB6SAGFREZnej70QjYFYDoB9Ho3KAzVviuQEuHltoeFhWBQUVERiM+LR7T\nD0zH6tOrUataLfzU7yeW+fSATgRVjRo1+EYhjahRo4a2h0A6SCEUWHNmDct8ekonFlMQEZUX1TKf\nVwMvhPiGsMxXAbiYgoioBPFp8ZhxYAZ+OP0Dy3x6jkFFRAZFIRRYe2Ytpu6fimcZzxDoGog58jks\n8+kxBhURGYxTD04hYFcATt4/Ca8GXljRawVa1Wql7WHRK2JQEZHey1/mC303FB+0+oBlPgPBoCIi\nvaUQCqw7uw5T901FUkYSxruOx5zOc2BtYa3toZEGMaiISC+plvk863sixDeEZT4DxaAiIr2SkJ6A\nGQdmYNWpVXCo6oAf3/0RH7b6kGU+A8agIiK9wDKf8WJQEZHOO/3gNAJ2BeDE/RMs8xkhBhUR6SyW\n+QhgUBGRDlIIBdafXY8p+6YgKSMJ4zqOw1z5XJb5jBSDioh0imqZr1P9TgjxDcGbtd7U9rBIixhU\nRKQTEtIT8PmBz/H9qe/hUNUBm97ZhMFvDmaZjxhURKRdyjLf1P1TkZCewDIfFcCgIiKtOfPwDEbv\nHM0yHxWLQUVEFY5lPioLBhURVRiFUGDDuQ2Ysm8Ky3xUagwqIqoQqmU+j3oeCPENQevarbU9LNID\nDCoiKleJ6Yn4/ODn+C76O9hXtcfGdzZiyJtDWOajUmNQEVG5yF/mG9txLObK58LGwkbbQyM9w6Ai\nIo078/AMAnYF4Pi94yzz0StjUBGRxiSmJ2LmwZn47tR3qFmlJst8pBEMKiJ6ZQqhwMZzGzFl3xTE\np8cjoH0A5nnPY5mPNIJBRUSv5OzDswjYFYBj947BvZ479vjuQZvabbQ9LDIgDCoieimqZT47Szts\n6LsBQ1oPgYnMRNtDIwPDoCKiMmGZjyoag4qISo1lPtKGEufoWVlZ5oMHDw51dXU97uHhceTq1avO\n+Z8THBw8pW3btmddXFxO7dy5s3f5DJWItCUxPRFjdo2By2oXxCTEYEPfDYgcHsmQogpR4oxq06ZN\n/vb29nGhoaGDIyMjPSdNmrR0x44dfZTHo6Oj2//2228DoqOj2yckJNh26dLlQO/evXeW77CJqCIo\nhAKbzm/CZ3s/Q3x6PEa3H4353vNZ5qMKVeKMav/+/T79+vX7AwA6deoUde7cObX/hdq1a5fv0KFD\nN5qZmWU7ODg82bJly8DyGiwRVZxzj87Bc70nhv81HK/bvY7Tn5zG8l7LGVJU4UqcUcXHx9vZ2dnF\nA4BMJhMymUyoHn/48KFjXFycva+v7660tLQqY8aMWdGiRYt/i3q9OXPm5H4sl8shl8tfevBEpHlJ\nGUmYeXAmVkavhJ2lHdb3XQ//1v5czUfFioiIQERERLm8dolBZWtrm5CUlGQDAEIIWf6gsrKyep6S\nklJt165dvklJSTatW7c+37179z3Vq1dPLuz1VIOKiHRHYWW+efJ5qGFZQ9tDIz2Qf+Ixd+5cjb12\nif+L5OPjs3/r1q3vAUB4eHgPLy+vw6rH3dzcjllbWz8DAEtLy3RLS8t0ExMThcZGSETlTrXM18S2\nCU6NPIXlvZYzpEgnyIQQxT4hKyvL3N/ff1NMTEyTatWqpYSGhg4GgKCgoCWbN2/2A4CJEycuO3fu\nXJsXL15UGjt27PKBAwduKfSbyWSipO9HRBUnKSMJsw7OQkh0COws7bC422KW+UgjZDIZhBAaafJY\nYlBpEoOKSDcohAI/nv8Rn+37DE/TnmKUyyjM957PGRRpjCaDijf8EhmZ84/OI2BXAI7cPQK3um4I\n+zAMbR3bantYREViUBEZCdUyn62lLdb9Zx2GthnKMh/pPAYVkYETQuDHf37E5L2TWeYjvcSgIjJg\nqmU+17qOEC61AAAe2ElEQVSu2P3hbrzl+Ja2h0VUJgwqIgOUlJGE2RGzseLkCtha2mLtf9ZiWJth\nLPORXmJQERkQ1TJfXGocRrWXyny2lrbaHhrRS2NQERmI84/OY8zuMYi6E8UyHxkUBhWRnnuW8Qyz\nImaxzEcGi0FFpKeEEAj9JxST907Gk9Qn+NTlU3zR5QuW+cjgMKiI9NA/j/9BwK4ARN2JQkenjtj5\nwU60q9NO28MiKhcMKiI98izjWe5qvhqWNbDm7TUY3nY4y3xk0BhURHqAZT4yZgwqIh3HMh8ZOwYV\nkY5SLfPZWNiwzEdGi0FFpGOEEPjpwk8I2hPEMh8RGFREOuXC4wsI2BWAyDuR6ODUATs+2AGXOi7a\nHhaRVjGoiHTAs4xnmHNoDpafWA4bCxusfns1Pmr7Ect8RGBQEWmVssw3ee9kPE55jP+6/BcLuixg\nmY9IBYOKSEvyl/n+9vubZT6iQjCoiCpYcmYy5kTMwbcnvmWZj6gUGFREFUQIgZ8v/IygvUF4nPIY\nn7T7BAu6LIBdFTttD41IpzGoiCrAxScXEbArAIdvH0b7Ou1Z5iMqAwYVUTnKX+b7oc8P+Pitj1nm\nIyoDBhVROWCZj0hzGFREGpa/zLd90Ha0d2qv7WER6S0GFZGGqJb5rC2sWeYj0hAGFdErEkJg88XN\nmLRnEh6nPMbIdiOxsMtClvmINIRBRfQK/n3yLwJ2BeDQ7UNwqePCMh9ROWBQEb2E5MxkzD00F98c\n/wbWFtZY1WcVPm77MUxNTLU9NCKDw6AiKgNlmS9oTxAepTximY+oAjCoiEopf5nvz0F/ooNTB20P\ni8jgMaiISvA887lU5jvxDapXrs4yH1EFY1ARFUEIgV8u/oJJeybhUcojjHhrBBb6LETNKjW1PTQi\no8KgIirEv0/+xZjdYxARG8EyH5GWMaiIVKiW+awqWeH73t9jxFsjWOYj0iIGFRGkMt+Wf7dg0p5J\nePj8Ict8RDqEQUVG798n/2Ls7rE4GHsQ7Rzb4Y/3/0DHuh21PSwi+n8MKjJazzOfY97hefj6+Ncs\n8xHpMAYVGR3VMt+D5w8w4q0R+NLnS5b5iHQUg4qMyqW4SxizawzLfER6pMT9B7KysswHDx4c6urq\netzDw+PI1atXnQt7nkKhMHFzczsWHh7eQ/PDJHo1zzOfY/LeyWj9fWuce3QO3/X+DidGnGBIEemB\nEoNq06ZN/vb29nHHjx93XbRo0dRJkyYtLex5y5cvH3v16lVnmUwmND9MopcjhMCWi1vwRsgbWHJ0\nCYa1GYZrY6/hU5dPeS2KSE+UGFT79+/36dev3x8A0KlTp6hz5861yf+cO3fu1A8LC+vZt2/fv4QQ\nsvIYKFFZXYq7hK4/dsWgrYNQu1ptHPv4GFa/vZrXooj0TInXqOLj4+3s7OziAUAmk4nCZkzjxo37\ndtmyZRMXL178WUkzqjlz5uR+LJfLIZfLyzxoouI8z3yO+Yfn43/H/werSlZY6bsSn7T7hDMoonIU\nERGBiIiIcnntEoPK1tY2ISkpyQYAhBCy/EEUGho6uFWrVheaNWt2Wfmc4l5PNaiINEkIgV///RUT\n90zEg+cP8HHbj/Glz5ewr2qv7aERGbz8E4+5c+dq7LVLDCofH5/9W7dufc/d3f1oeHh4Dy8vr8Oq\nx6OiojpdvHixpbe398ErV668cebMmbesrKyeu7u7H9XYKIlKcDnuMsbsHoMDtw7gLce3sPX9rXCt\n66rtYRGRBsiEKH7tQ1ZWlrm/v/+mmJiYJtWqVUsJDQ0dDABBQUFLNm/e7Kf63OHDh6/38/Pb3L17\n9z2FfjOZTJT0/YjKIuVFCuYdmof/Hf8fqlWqhoVdFrLMR6Rt8fGQ1axZYoWttEoMKk1iUJGmCCHw\n26XfMDF8Iu4/v88yH5E23b4NREYCUVHSfy9dggwlXwoqLQYV6R3VMl/b2m0R4hsCt3pu2h4WkXFQ\nKIBLl9SD6e5d6Vj16oC7O+DpCdmMGQwqMj4pL1Iw//B8LDu2jGU+oory4gVw+nReMB05AiQkSMdq\n1wY8PfMerVoBptLfR5lMxqAi45G/zPdR24+wyGcRy3xE5eH5c+DYsbxgOnECSE+XjjVtCnTqlBdM\njRsDssKziEFFRuPK0ysYs2sM9t/azzIfUXl4/DivhBcVBZw9K5X3TEyAtm3zgqlTJ6BWrVK/LIOK\nDF7KixR8cfgLLDu2DFUrVcWCLgvw33b/ZZmP6FUIAdy8KYWSMpiuXZOOWVgArq55weTmBlhZvfS3\nYlCRwRJC4PdLv2Pinom4l3wPw9sMx6Kui+BQ1UHbQyPSPzk5wIUL6sH08KF0rEYNKZSUwdSuHVCp\nksa+tSaDitt8kM648vQKxu4ei30396FN7Tb4tf+vLPMRlUVGBhAdnRdMR48CycnSsXr1AG/vvDJe\n8+ZSeU8PcEZFWpe/zPeF9xfsbk5UGklJUhgpgyk6WlqlB0hBpFz00KkT0KBBhQ6NpT8yCCzzEZXR\n/ft5Cx8iI6WynhCAmZlUulMGk7s7UFO7uwQwqEjv5S/zrfRdyTIfkSohgKtX1W+svXVLOla1qrTY\nQRlMHTpIn9MhDCrSW6plvirmVbCgywKW+YgAIDtbWhquDKaoKCAuTjpmb69+/1KbNtIsSocxqEjv\nCCGw9fJWTAifgHvJ9zCszTAEdw1mmY+MV1oacPx4XjAdOwakpkrHGjdWD6amTYu8sVZXMahIr1x9\nehVjd4/F3pt70aZ2G4T4hsC9nru2h0VUseLj1W+sPX1amkXJZMCbb6rfWOvkpO3RvjIGFemF1Bep\n+CLyCyw9uhRVzKvgiy7Saj4zE90uWRBphLKjuDKYLl2SPl+pknRNSRlM7u6AjY12x1oOGFSk05Rl\nvonhE3E3+S6GtRmGRT6LUKta6duvEOkV1Y7iymBS7Sju4ZEXTO3bS10gDByDinSWapmvda3WCPEN\ngUd9D20Pi0izVDuKR0ZKHcUTE6Vjjo55Jbx8HcWNCYOKdE7qi1QsiFyAJUeXsMxHhke1o3hkpNRR\nPCNDOta0qfqNtcV0FDcmDCrSGUII/HH5D0wIn8AyHxkO1Y7ikZHAuXPqHcWVweThUaaO4saEQUU6\n4Vr8NYzdPRZ7buxhmY/0lxDAjRvqwXT9unRM2VFcGUyurq/UUdyYMKhIq/KX+eZ7z8eo9qNY5iP9\nkJMD/POPeseHR4+kY8qO4spgeustjXYUNybsnk5aIYTAtivbEBgWiLvJdzG09VAEdw1mmY90W0YG\ncPJkXjDl7yjepUteMDVrpjcdxY0JZ1RUKqplvjdrvYkQ3xB0qt9J28MiKigpSVqFpwwm1Y7iLVqo\nz5jq19fuWA0YS39UYVJfpGJh1EJ8deQrWJpb4gvvL1jmI91y/776/UuqHcVdXPKCycMDsLPT9miN\nBoOKyl3+Mp9/a38s7rqYZT7SLtWO4spgUu0o7u6eF0wdOwJVqmh3vEaMQUXl6lr8NYzbPQ7hN8JZ\n5iPtysqSloarBtPTp9Ixe3v1G2v1oKO4MeFiCioXyjLfkqNLYGFmgW96foPR7UezzEcVJzVVuplW\nGUzHj6t3FO/dOy+Y9LCjOL0czqgIQgj8eeVPBIYH4s6zO/Bv7Y/grsGoXa22todGhu7p07yFD5GR\nwJkz6h3FVTs+1Kmj7dFSGbD0RxpzPf46xu4ei/Ab4Wjl0AohviHwbOCp7WGRIRJC6iiuemPt5cvS\nMWVHcWUwubkZZEdxY8KgoleWlpWGhZEL8dXRr2BhZoH53vNZ5iPNUiiAf/9VD6Z796Rjyo7iymBy\ncTGKjuLGhNeo6KXlL/MNeXMIFndbzDIfvboXL4BTp/IWPRTWUVz5aNnSKDuK08thUBmR6/HXMS5s\nHMJiwtDKoRUODzvMMh+9vOTkvI7iUVHqHcWdnYF+/fKCqVEjLnygl8bSnxHIX+abJ5+HgA4BLPNR\n2Tx+rL5MXNlR3NRU6iiuupW6g4O2R0taxmtUVCpCCPx19S8EhgXi9rPbLPNR6Sk7iqsGk7KjuKWl\n1EVcGUzsKE6F4DUqKlH+Mt+hYYfg1cBL28MiXaXaUVwZTMqO4ra2UiiNHMmO4qQVDCoDk5aVhi+j\nvsTiI4tR2bQyvu7xNct8VFB6utSsVRlMR49Ku9gCUqNWH5+8GRM7ipOWsfRnIPKX+Qa/ORiLuy6G\no5WjtodGuiAxUQojZTCdOqXeUVy1FRE7ipMG8BoVqYlJiMG43eOwO2Y3Wjq0RIhvCMt8xu7ePfX7\nly5eVO8orlyN5+7OjuJULhhUBEAq8y2KWoTgI8GobFoZ87znIaB9AMxNzbU9NKpIQgBXrqgHU2ys\ndKxaNanLgzKYOnRgR3GqEAwqIyeEwPar2zE+bDzLfMYoKws4ezZv0UNhHcWVj9at2VGctKJCV/1l\nZWWZDx8+fH1MTEwTU1PTnHXr1n3k7Ox8VXk8MzOzsr+//6bY2NiGJiYminnz5s3q1q3bXk0MjgqK\nSYjB+LDx2HV9F1o6tORqPmOQmip1EVcG07FjQFqadEzZUVwZTK+/zhtryeCUGFSbNm3yt7e3jwsN\nDR0cGRnpOWnSpKU7duzoozy+efNmv5o1az7dsmXLwKdPn9Z0d3c/eu3atablO2zjk7/Mt6z7Mozp\nMIZlPkP09GleGS8qSr2jeOvWwEcfsaM4GZUSg2r//v0+o0aN+g4AOnXqFOXn57dZ9XjDhg1j27Vr\ndxoALCwsMlJSUqqVz1CNU/4y34etPsRX3b5imc9QKDuKq96/pOwoXrmydE1p8uS8hQ/W1todL5EW\nlBhU8fHxdnZ2dvGAdI1JJpOpXWSSy+URAHDx4sWWn3zyyQ9BQUFLinu9OXPmqH4t5HJ52UdtJG4k\n3MC4sHHYdX0XWti3QMTQCHRu2Fnbw6JXoeworhpMyo7i1tZSR/EhQ9hRnPROREQEIiIiyuW1Swwq\nW1vbhKSkJBsAEELI8gcVAMybN2/W1q1b3/v6668Dvb29Dxb3eqpBRYVTlvkWH1mMSqaVWObTZ5mZ\nwOnTecF05AiQlCQdq1NH/f4ldhQnPZZ/4jF37lyNvXaJQeXj47N/69at77m7ux8NDw/v4eXldVj1\n+ObNm/1OnTrlEh0d3b5SpUovNDYyIySEwN/X/sb4sPGITYplmU8fqXYUj4wETp5U7yjev39eMLGj\nOFGplLg8PSsry9zf339TTExMk2rVqqWEhoYOBoCgoKAlmzdv9hs6dOjGM2fOvFWzZs2ngFQePHDg\nQJdCvxmXpxfpRsINjA8bj53Xd6KFfQuE+IawzKcPHj1Sv3/p/Hn1juLK1XgeHuwoTkaF91EZkPSs\ndCw6sgjBUcGoZFoJc+RzMLbDWJb5dJEQQEyMejDFxEjHlB3FlcHk6irdbEtkpBhUBuLvq39jXNg4\nxCbF4oNWH+Crbl+hjhWXG+uMnBxphqS6VDx/R3FlMLVty47iZDQyMoD796V1QMrH3bvqf378mNt8\n6DXVMl9z++Y4OPQg5A3l2h4WpadL15SUoVRYR3FlML3xBjuKk0FKT88LIdXwUf04Lq7g19nYAPXq\nAXXrSjvBrF6tuTFxRlWB0rPSEXwkGIuiFsHc1Bxz5XNZ5tOmxERpFZ4ymKKjpfZEQF5HceWqPHYU\nJwOQlqY+6yksjOLjC36dra0UQHXr5oWR6sdOTgUr3Sz96aG/r0qr+W4l3WKZT1vu3VO/f0nZUdzc\nXLpnSVnK8/CQ/mYS6ZHU1KLDR/lxQkLBr7OzKzx8VB8v08eYQaVHbibexPiw8dhxbQea2zdHiG8I\ny3wVQdlRXDWYVDuKu7vnBRM7ipOOS0kpOnyUf1benqfK3r7kmZClZfmMmUGlB/KX+eZ0noNxHcex\nzFdeVDuKK4NJWcNwcFC/sZYdxUmHJCcXHT7Kj589K/h1Dg7Fz4ScnLTb2IRBpeN2XNuBcbvH4VbS\nLfi19MOS7ktY5tM01Y7ikZHSx8qO4q+9ph5M7ChOWiCEFEIlzYSU63WUZDKgVq3Cw0f55zp1pFaQ\nuoxBpaPyl/lW9FoB70be2h6WYYiLy1v4EBkpdRTPycnrKK4MJnYUpwoghFRqK2kmlJKi/nUyGVC7\ndvEzoTp1DONOBwaVjknPSsfiI4vxZdSXLPNpghDS9STVG2uvXJGOKTuKK1fkubmxozhplBDSooPi\n7hG6ezdvAq9kYgI4OhY/E3J0lNbuGAMGlQ5RLfMNajkIS7otgVN1J20PS78oFNIKPNVgun9fOqbs\nKK4MJhcX3a95kM4SQrp0WdJMKD1d/etMTKSZTnEzIUdHXvpUxaDSATcTbyIwLBB/X/sbzWo2Q4hv\nCMt8pZWZCZw6lRdMhXUUVz5atGBHcSoVIaQKcXH3CN27J739VJmaSgsPipsJ1arFECorBpUW5S/z\nze48G+M7jmeZrzjJyVKXB+VqvPwdxVWDqWFDLnygAhSKvBAqbib0It/+DWZmUggVNxOqVYv/L1Qe\nGFRasvPaTowLG4ebiTdZ5ivOo0fqy8RVO4q/9VbearxOnaQbPcioKRTA48fFz4Tu389rGqJkbl74\n7Ef1YwcHdrrSFgZVBbuVeAvjw8bnlvlW+K5Al0aF7mRifJQdxVWDSbWjuJtbXjCxo7jRycnJC6Gi\nlmnfvw9kZ6t/XaVKRYeP8mFvzxDSZQyqCpKRnZFb5jOVmWKOXFrNV8nUANaOvixlR3HVYHr8WDpm\nZ5e3RNzTU5o9GcsSJyOUkwM8fFj8TOjBA+l5qiwsSp4J1azJCrC+Y1BVANUy38AWA7Gk+xLUrV5X\n28OqeKodxSMjpd1rlXcoNmigfmMtO4objOzsvBAq6prQw4cFQ8jSsuSZkJ0dQ8gYMKjK0a3EWwgM\nD8T2q9uNs8yXkJC38CEyUlqdp7w40LKlejDVq6fdsdJLycqSZjrFzYQePZKuHamqUqXoAFJ+XKMG\nQ4gkDKpykL/MN7vzbIx3HW/4Zb67d9XvX7p4Ufq8sqO4MpjYUVwvvHiRF0JFXRN69Ei6tKiqalUp\nbIqbCdnYMISo9BhUGrbz2k6MDxuPG4k3DLvMJwRw+bJ6MN2+LR1TdhRXLhNv354dxXVMZqb6rqqF\nleSUlwtVVa9e8jWh6tUZQqRZDCoNiU2KRWBYIP66+hfeqPkGVvRaAZ/GPtoeluZkZUk98VS3Us/f\nUVz5ePNN3tGoRapbexc1E3rypODXWVuXfE2oevWK/3mIGFSvKCM7A18d+QoLoxbCVGaKWZ1nIdA1\nUP/LfCkpeR3Fo6IK7yiufDRpwv+FriDp6SX3jXv6tODX1ahRdPjUqyfdyGplVfE/D1FpMKhewa7r\nuzBu9zjcSLyB91u8j6Xdl+pvmS8uTn22VFhHceU1JkdHbY/WICm39i6uW0JxW3sXNRMqbGtvIn3C\noHoJel/mU3YUV71/SbWjeMeOeavx2FFcI1JSSp4JJSYW/LqaNUueCfHyHxk6BlUZ6G2ZT9lRXDWY\nlB3FbWykVXjKYGJH8TJ7/rzkmVBxW3sXNxMqr629ifQJg6qUdl/fjbG7x+pHmU/ZUVwZTEeO5O0/\n7eSkfv9Sy5a8sbYYz56VPBNKTi74dcpdVYtaHVenjna39ibSJwyqEqiW+ZztnLHCdwW6Nu5a7t+3\nTJ49k7o8KIPp5Mm8/QfeeEM9mNhRHIBU/VSGUHEzodJu7a36sT5s7U2kTxhURcjIzsCSo0uwIHKB\n7pX5Hj5Uv3/pn3/UO4qrbqVuhB3FhZCu95S0l1BqqvrXyWR5u6oWdU3I0dEwtvYm0icMqkKolvkG\nNB+Apd2Xop61llr8CAFcv64eTDduSMeqVJG6iCtX5HXsaPDLu1S39i5uJlTarb1VPzamrb2J9AmD\nSkVsUiwmhE/An1f+1F6ZLztb6iiuulQ8f0dxZTC1bWtQ/7IKId0DVNJMSLlPopKpqVRuK+6aUO3a\nvAeZSF8xqKBe5jORmWCW1yxMcJtQMWW+9HTgxIm8YDp6VFrLDOR1FFc+nJ31duGDQpEXQkV1Syhs\na28zs7wQKmomxK29iQyb0QdVhZf5EhKkVXjK2VJhHcWV15j0pKO4QiG15ClpV9X8W3ubm0uLEIub\nCTk4cGtvImNntEGVv8y3vNdydHutmwZH+P/u3lW/f6mwjuKentK9TDVqaP77v6KcnLwQKm5X1eK2\n9i5qJsStvYmoNIwuqJRlvoWRCyGTyTRb5lN2FFcNJmVHcSsrqaO48hpThw5av5szJ0fapqG4e4Qe\nPCi4tXflykXPgpR/rlmTIUREmmFUQRUWE4axu8ciJiEG/Zv3x7Luy16tzKfsKK56Y62yGVutWur3\nL1VwR/HsbCmEilsZV9LW3kXNhLi1NxFVJKMIKtUyX1O7pljRa8XLlfmK6yjepIn6irxy7Ciena2+\nq2phYfTwYcFdVS0ti9/CoV49qcEpQ4iIdIlBB1Vmdmbuaj6ZTIaZXjMxwXUCKpuVsm1AXFzewofI\nyArpKK7c2ru4mVBJW3sXFUbc2puI9JHBBlX+Mt/S7ktR37p+0S8ohHQ9SRlKkZHqHcU7dMgLppfs\nKJ6ZWfJM6PHjglt7V6tW+OxH9WNra4YQERkmgwuq20m3MSF8ArZd2Yamdk2xvNdydH+te8EXUCiA\nf/9VX/hw7550zNpaWoWnnC21b19i8zbl1t7FzYRK2tq7qJkQt/YmImNmMEFVYpnvxYu8juJRUVJJ\nT7kBkKOj+o21LVuq3byTnp63tXdRy7Tj4gqO0dq6+JkQt/YmIiqZQQRVeEw4xuweg5iEGLzX7D0s\n67EM9U1qqHcUP3Eir/dO06a5oZTu4om7Zo1w776syJJccVt7598/qF49w9/aOyIiAnK5XNvDMEo8\n99rDc689mgyqEtdeZ2VlmQ8fPnx9TExME1NT05x169Z95OzsfLW0x/NTLfO5mTfGbzVnos3RZCD4\nXeDcOUChgDAxwfMmbXHX61NcsvPESTMPXE6ohbvRwL1tUqOI/Ozs8sLH1bXwXVWrVn3Z06T/+BdW\ne3jutYfn3jCUGFSbNm3yt7e3jwsNDR0cGRnpOWnSpKU7duzoU9rj+c3/+HW8HavA8rs14RR3E8B8\nZJpY4LylKyLMZ2BfZiccU7gh5ZoVcE36GuWuqg0aSJeh8pfi6tbV+n24RERUTkoMqv379/uMGjXq\nOwDo1KlTlJ+f3+ayHM9vzR9ZSER1RMEV38ATF2t4IqFhO9SuXwl16wJd6gJD85XmuKsqEZERE0IU\n++jevXv4v//+21z557p1694ty3HVBwDBBx988MGHcTxKypfSPkqcUdna2iYkJSXZQPquMplMJspy\nXJWmLqwREZHxKLEFqY+Pz/6tW7e+BwDh4eE9vLy8DpflOBER0asocXl6VlaWub+//6aYmJgm1apV\nSwkNDR0MAEFBQUs2b97sV9hxJyen+xUyeiIiMngVeh8VERFRWXH3ISIi0mnlHlRZWVnmgwcPDnV1\ndT3u4eFx5OrVq87l/T2NlZeX12Fvb++D3t7eB8eMGbMiOzvbrLBzz9+J5mzZsmXgtGnTvgSKPq+F\nfZ6/g1eneu4B9fd/QEBACMBzr2mZmZmVBw4cuKVjx44n3Nzcju3du7dbWf6deelzr6nlg0U91qxZ\n83FgYOD/hBA4fPiwZ+/evXeU9/c0xkdKSkrV/Oe2qHPP38mrPxQKhaxr1657LSws0qdNm7awrOd7\n7dq1H/F3oLlzX9j7n+de84/169cPGz16dIgQAnFxcTWbNGlyvajzqclzX+4/mJ+f38+HDx/2VL7B\nnJyc7mn7ZBvi4/Tp02+1aNHiYpcuXfZ369Ztz8mTJ9sXde75O9HMIzs723TdunXDp06d+mVx57Ww\nz/N3oNlzn//9Hx0d7cJzr/nHwYMH5f/8808rIQSeP39ezdHR8UFFvO/LfZ/1+Ph4Ozs7u3hAakpb\n3H1W9PLMzc2zxo4du/y///3vqsuXLzfr3bv3zkaNGt2ytbVNAPLOvRBCxt+JZpiamuaonrvCzmth\n5xsAEhISbPP/brTxM+grU1PTHBMTk9ytSPO///v06bPj2rVrTXnuNUsul0cAwMWLF1t+8sknPwQF\nBS0JDw/vUZp/Z4CXP/flHlRluSGYXl7Lli0vtmrV6gIANGvW7HLNmjWfmpqa5jx79swayDv3MplM\n8HeiOarnrrDzWtj5NjExUdja2ibk/91o5ycwDIW9/x89elSb517z5s2bN2vr1q3vff3114He3t4H\no6Oj25fm35lXOfflvpiCNwRXjC+//HLanDlz5gDAw4cPHZ8/f241YMCA3wo79/ydlI+izmv+z3t6\nekZ26dLlAH8HmrNw4cLpqu//5OTk6o6Ojg957jVr8+bNfqdOnXKJjo5u7+3tfRCooPd9edc0X7x4\nYT5o0KDNLi4u0XK5/OC9e/ectF1nNcRHcnKy1dtvv729U6dOkV5eXociIyM7FXXu+TvR3GPDhg1D\nlRf0y3K++TvQ7Lkv7P3Pc6/5h7+//8aWLVtekMvlB+Vy+UFvb+8DFfG+5w2/RESk03jDLxER6TQG\nFRER6TQGFRER6TQGFRER6TQGFRER6TQGFRER6TQGFRER6TQGFRER6bT/A938sMQJ3WYxAAAAAElF\nTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x38a2e90>"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Next\n",
      "\n",
      "In Part 4 we'll look at using the **Scipy** ODE integration pack and the **QuTiP** solver for the specific problem of solving the Schr\u00f6dinger equation (for a **closed quantum system**) and Lindblad master equation (for an **open quantum system**)."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}