MAT301_DifferentialEquationSolving2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nlr1Ru5KqTRX"
   },
   "source": [
    "# Differential Equation Solving 2\n",
    "\n",
    "Notebook to implement the 4th order Runge-Kutta scheme to solve differential equations for Abertay University undergrad course MAT301 (JT, 2022).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Nopr9bZQUGCZ"
   },
   "source": [
    "## 4th Order Runge Kutta Algorithm\n",
    "As seen in lectures, a highly common approach to solve differential equations is to use the RK4 algorithm. This is one of a family of Runge-Kutta (RK) algorithms, which are based on the same mathematical ideas. RK4 is an industry standard, representing a good trade-off between accuracy and complexity.\n",
    "\n",
    "We apply Runge-Kutta algorithms to differential equations which take the same form, namely:\n",
    "\n",
    "$$ \n",
    "y^{'} = f(x,y). \n",
    "$$\n",
    "\n",
    "The algorithm uses the values of the independent variable $x$ and dependent variable $y$ at specific values ($x_n$ and $y_n$) in order to approximate the next value of $y$ using:\n",
    "\n",
    "$$ \n",
    "y_{n+1}=y_n+hk,\n",
    "$$\n",
    "\n",
    "for $n=0,...,N-1$. The form of the coeffient $k$ varies depending on the particular RK scheme used; for RK4 we note: \n",
    "\n",
    "$$\n",
    "k=\\frac{1}{6}\\left[ k_1+2k_2+2k_3+k_4\\right],\n",
    "$$\n",
    "and individual components $k_1\\rightarrow k_4$ are determined by evaluating the function at different locations according to:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "k_1&=f\\left(x_n,y_n\\right),~~~~~~~~ k_2=f\\left(x_n+\\frac{1}{2}h,y_n+\\frac{1}{2}hk_1\\right), \\\\\n",
    "k_3&=f\\left(x_n+\\frac{1}{2}h, y_n+\\frac{1}{2}hk_2\\right),~k_4=f\\left(x_n+h, y_n+hk_3\\right).\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "In all of these examples $h$ is the stepsize with which we step along in $x$.\n",
    "\n",
    "We have already seen in the lectures how the method can be applied to different differential equations. However, the real power of this method lies in its iterative behaviour; we can use computers to crank up the number of steps of calculation (or decrease the step size) to calculate an (approximate) solution over a large range of $x$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RaEolGImUbMw"
   },
   "source": [
    "To illustrate how this method can be implemented numerically, we will apply it to three initial value problems. These happen to be models of population growth over time, but could in reality be for any type of ordinary differential equation (ODE, or indeed for a system of several coupled ODEs). I have selected these equations because they are broadly related, but gradually increase in complexity. The equations which we will investigate are the following:\n",
    "### 1. Linear Population Equation\n",
    "Our first ODE will be linear:\n",
    "\n",
    "$$\n",
    "y^{'}=0.1y, \\ \\ (2000 \\leq t \\leq 2020), \n",
    "$$\n",
    "\n",
    "with the initial condition,\n",
    "\n",
    "$$\n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "### 2. Non-Linear Population Equation \n",
    "Our second ODE will contain more terms and a higher power in $y$:\n",
    "\n",
    "$$\n",
    "y^{'}=0.2y-0.01y^2, \\ \\ (2000 \\leq t \\leq 2020), \n",
    "$$\n",
    "\n",
    "with the initial condition,\n",
    "\n",
    "$$\n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "### 3. Non-Linear Population Equation with an oscillation \n",
    "Our final ODE will include an extra complication of an oscillatory part added to the previous ODE: \n",
    "\n",
    "$$\n",
    "y^{'}=0.2y-0.01y^2+\\sin(2\\pi t), \\ \\ (2000 \\leq t \\leq 2020), \n",
    "$$\n",
    "\n",
    "with the initial condition,\n",
    "\n",
    "$$\n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "We will tackle each of these equations in turn, applying the RK4 algorithm to each problem and comparing results with any known analytical solutions of the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5BKUkz0jqTRZ"
   },
   "source": [
    "#### Setting up Libraries\n",
    "\n",
    "As always, we'll make use of some clever python tools for plotting and maths, so our first step is to load in the libraries that store these tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "F6L7k675qTRa"
   },
   "outputs": [],
   "source": [
    "## Libraries\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import math \n",
    "import pandas as pd #good library for tables\n",
    "import matplotlib.pyplot as plt \n",
    "import matplotlib.gridspec as gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pzykh1q-qTRf"
   },
   "source": [
    "## Discrete Interval\n",
    "The independent variable, which for these equations is time, $t$, has a particular range $a\\leq t \\leq b $. We divide this interval into $N$ points separated by a constant stepsize using the equation stated in our lectures:\n",
    "\n",
    "$$\n",
    "h=\\frac{b-a}{N}.\n",
    "$$\n",
    "\n",
    "Here, the interval is $2000\\leq t \\leq 2020,$ \n",
    "\n",
    "$$\n",
    "h=\\frac{2020-2000}{200}=0.1.\n",
    "$$\n",
    "\n",
    "This gives the 201 discrete points:\n",
    "\n",
    "$$ \n",
    "t_0=2000, \\ t_1=2000.1, \\ ... t_{200}=2020, \n",
    "$$\n",
    "\n",
    "which we can generalise to\n",
    "\n",
    "$$ \n",
    "t_n=2000+0.1n, \\ \\ \\ n=0,1,...,200.\n",
    "$$\n",
    "\n",
    "Let's go ahead and set this range of $t$ values as a list in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "OYX2NEUfqTRg"
   },
   "outputs": [],
   "source": [
    "N = 200\n",
    "t_end = 2020.0\n",
    "t_start = 2000.0\n",
    "h = ((t_end - t_start) / N)\n",
    "t = np.arange(t_start, t_end + h / 2, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q6U6K7iZa-nV"
   },
   "source": [
    "Normally, I'd print these, but for this massive range of times I think it would be better to visualise these with a plot: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "executionInfo": {
     "elapsed": 245,
     "status": "ok",
     "timestamp": 1646995514491,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "RKt46sT4a_XP",
    "outputId": "0681760b-0b5e-4ba8-9002-970f2d49e678"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAEICAYAAABoNzG1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAarklEQVR4nO3debRlZX3m8e8jpdA4MRUyU2gwCRqcrlM7xQGjOKCJ3Q1NEG0V6Wha13K1YtAMzro6ii60ScUJFcVZiWKrcUpURAshRERliEpBCSUygwP46z/2vnK4nFO3bp3z1j2H+/2sddfdZ+93v8PeZ3Oeevc+l1QVkiRJauN2y90BSZKk2zLDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2NKKlORZSb4+8LqS/N5y9mmUJCckeeUytPs/k1ya5NokO29G+R8neVy//FdJ3tm+l+0kOTzJF5a7H8stycOSnNe/D542gfrem+Q1E+iaNDMMW1oxBsNAo/rH/hBZGAIBquroqnr1eL1bcj9uD7wZeHxV3amqLl/K/lX1uqp6bpve3VqSNX1gXjWp/avqpKp6/OR6ObNeBRzfvw8+tdydmdefs68kuT7JDzZ1bSd5dF/2qiQ/3ordlADDlrTVbGkQWCZ3A7YDzlnujsDMHbuxTOFY92UL3weNx/Ih4ExgZ+BY4GNJVo8oex3wbuB/N+yPNJJhS1ogyVeTPHfg9e9mm9J5S5LL+n8ln53k3kmOAg4HXtrfbvmnvvyPk7wsydnAdUlWJTkmyQVJrkny/SRP78v+IXAC8NC+jiv79beYMUvyvCTnJ/lFklOS7DGwrZIc3d/2uSLJ25NkxDi3TXJckkv6n+P6dfcEftgXuzLJl0fsf0SSnyS5PMmxC7b9bZIP9MvbJflAX+7KJN9Jcrd+205J3tO3f0WST/Xr/zjJ+v7Y/Qx4T5LbDRy7y5N8JMlOfZP/MtDfa5M8tK/nfyQ5t6/780n2HXHab7X/wlnG/tj+RX9sr0ny6iT3SHJakqv7/txhoPyTk5zVj/mbSQ4c0fZ83S9Ich5wXr/urUku6us+I8kjFhzfjyR5X9+Xc5LMDWy/f5Iz+20fTfLhBe+hzepbkguAuwP/1B+XbZPs0b/vftG/D5+3oF8f68/31cCzRgx5xySf7ft3epJ7jDo2I/p1T+D+wN9U1Q1V9XHg34E/G1a+qr5dVe8HLlxKO9KkGLakpXk88EjgnsAOwH8DLq+qtcBJwJv62y1PGdjnMOBJwA5VdSNwAfAI4K7A3wEfSLJ7VZ0LHA2c1texw8LGkzwGeD3wX4HdgZ8AJy8o9mTggcB9+nJ/MmIsxwIPAe7bl30Q8Iqq+hFwr77MDlX1mCH9OAD4v8ARwB50swt7jWjnyH6se/fljgZu6Le9H9i+b29X4C0D++0G7EQ3s3IU8L+ApwGP6tu8Anh7X/aRA/29U1Wdlu75or8C/hRYDfwr3WzIMLfaf0S5JwAPoDtuLwXW0oXsvYF7051rktyfbibl+f2Y/wE4Jcm2I+qlH9uDgQP619+hOzc7AR8EPppku4HyT6U79zsApwDH923fAfgk8N5+3w8BT5/faSl9q6p7AD8FntIfl1/19a2nOwfPAF6X5LEDux0CfKzv10kjxnoY3Xt/R+B84LUD/Tu7D4HDft7RF7sXcGFVXTNQ579x8/tWmiqGLWlpfgPcGfgDIFV1blVtWGSft1XVRVV1A0BVfbSqLqmq31bVh+lmMh60me0fDry7qr7bf/C9nG4mbM1AmTdU1ZVV9VPgK3Qf2KPqelVVXVZVG+k+/I7YzH48A/hMVf1L349XAr8dUfY3dB/qv1dVN1XVGVV1dZLdgScCR1fVFVX1m6r62sB+v6WbufhVf+yeDxxbVev7Nv8WeEZG36p6PvD6/hzdCLwOuO8mZrc2xxur6uqqOgf4HvCFqrqwqq4CPgfcry/3POAfqur0fswnAr+iC2mjvL6qfjHwPvlAVV1eVTdW1d8D2wK/P1D+61V1alXdRBda79Ovfwiwiu5995uq+gTw7YH9tqRvACTZG3g48LKq+mVVnQW8k1u+b06rqk/17+8bhtUDfKKfbbqRLpDdd35DVR1YVTuM+PmLvtidgKsW1HkV3bUpTR3DlrQEVfVluhmEtwOXJlmb5C6L7HbR4Iskzxy4hXMl3YzILpvZhT3oZrPm+3MtcDmw50CZnw0sX0/3wbRoXf3yHiPKDtv3d+Oqquv6fgzzfuDzwMn97cI3pXsAf2/gF1V1xYj9NlbVLwde7wt8cuC4nQvcRPd82TD7Am8dKP8LINzyWC3VpQPLNwx5PX+s9wVeMjgrQzfeTR3fhe+Tl/S3QK/q978rt3yfLDzP2/XBcw/g4qqqEXVvSd/m7UF3zgZnlH7CLY/pRSxuc9+jo1wLLLzu7gJcM6SstOwMW9KtXUd3a2veboMbq+ptVfUAulsW9+Tmh24HP9xuscv8Qj+r8o/AC4Gd+1uF36MLAZuqY94ldB+W8/XdkW7W6OJF9lu0LmCfft3m2ED3AT3fj+37ftxKP7vyd1V1APCf6W5zPpPuQ3mnJDuMaGPhsbgIeOKCmY7tquriIWXnyz9/Qfn/VFXf3Iy2xnUR8NoFbW9fVaNuY96iD/3zWS+juw28Y/8+uYqb3yebsgHYM7nFs3p7DyxvSd/mXUJ3zgZnkPbhlu+/sY5l//zZtSN+TuiLnQPcfUE/7sOUfKFDWsiwJd3aWcCfJtk+3d/ees78hiQPTPLgfmbmOuCXdLMr0M1y3H2Ruu9I92G0sa/v2XQzW/MuBfYafNB6gQ8Cz05y3/4Zm9cBp1fVj5cwvnkfAl6RZHWSXYC/Bj6wmft+DHhykof3fX0VI/57ku5r93+UZBvgarrbijf1t18/B7wjyY5Jbp/kkcPq6J0AvHb+NmDf70P6bRvpbjvefUH5lye5V1/+rkn+y4i6h+0/jn8Eju7fK0lyxyRPWhAONuXOwI19v1Yl+WtuPZMzyml078kXpvtCxiHc8jb1Fvetqi4Cvgm8Pt0XHw6kuz5GPZu1ZFV1r/75sGE/R/dlfkR3nf5N34+nAwcCHx9WZ7ovV2wH3L57me02cY1JE2fYkm7tLcCv6YLPidzyg+QudB9WV9DdPrkc+D/9tncBB/S3Zj41rOKq+j7w93QfiJcCfwR8Y6DIl+n+df6zJD8fsv+X6J6P+jjdDMY9gEO3ZJDAa4B1wNl03+T6br9uUf0zSy+gC38b6I7H+hHFd6MLZ1fT3fr7GjeHuiPowtcPgMuAF2+i2bfSPQj+hSTXAN+ie6Ccqrqe7iHrb/TH/yFV9UngjXS3L6+mm0F84ojx3Gr/xY7BplTVOrpno46nOzbnM/qbecN8ni6I/ojuffZLNu/2HFX1a7ovBTwHuBL4c+AzdM9lTaJvhwFr6Ga5Pkn3XN0Xl7D/pBwKzNGN4Q3AM/pnD0nyiCTXDpR9JN1t3lPpZuJuAFb8H6zV1pNb3taXJN3WJDkdOKGq3rPcfZFWIme2JOk2JsmjkuzW30Y8ku4W2/9b7n5JK9W0/aViSdL4fh/4CN23/C6gu8W22J8okdSItxElSZIa8jaiJElSQ1N9G3GXXXapNWvWLHc3JEmSFnXGGWf8vKpu9T9En+qwtWbNGtatW7fc3ZAkSVpUkp8MW+9tREmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhiYStpI8IckPk5yf5JhNlHtgkpuSPGMS7UqSJE27scNWkm2AtwNPBA4ADktywIhybwQ+P26bkiRJs2ISM1sPAs6vqgur6tfAycAhQ8r9JfBx4LIJtClJkjQTJhG29gQuGni9vl/3O0n2BJ4OnLBYZUmOSrIuybqNGzdOoHuSJEnLZxJhK0PW1YLXxwEvq6qbFqusqtZW1VxVza1evXoC3ZMkSVo+qyZQx3pg74HXewGXLCgzB5ycBGAX4OAkN1bVpybQviRJ0tSaRNj6DrB/kv2Ai4FDgf8+WKCq9ptfTvJe4DMGLUmStBKMHbaq6sYkL6T7luE2wLur6pwkR/fbF31OS5Ik6bZqEjNbVNWpwKkL1g0NWVX1rEm0KUmSNAv8C/KSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGppI2EryhCQ/THJ+kmOGbD88ydn9zzeT3GcS7UqSJE27scNWkm2AtwNPBA4ADktywIJi/wE8qqoOBF4NrB23XUmSpFkwiZmtBwHnV9WFVfVr4GTgkMECVfXNqrqif/ktYK8JtCtJkjT1JhG29gQuGni9vl83ynOAz02gXUmSpKm3agJ1ZMi6GloweTRd2Hr4yMqSo4CjAPbZZ58JdE+SJGn5TGJmaz2w98DrvYBLFhZKciDwTuCQqrp8VGVVtbaq5qpqbvXq1RPoniRJ0vKZRNj6DrB/kv2S3AE4FDhlsECSfYBPAEdU1Y8m0KYkSdJMGPs2YlXdmOSFwOeBbYB3V9U5SY7ut58A/DWwM/COJAA3VtXcuG1LkiRNu1QNfbxqKszNzdW6deuWuxuSJEmLSnLGsMkk/4K8JElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ1NJGwleUKSHyY5P8kxQ7Ynydv67Wcnuf8k2h3LSSfBmjWQwKpV3e9ddul+bne7m5eHbd8a66ahD9PaL/sw3f2yD9Pdr2now7T2yz5Md7/G7cOaNd1n/zJIVY1XQbIN8CPgIGA98B3gsKr6/kCZg4G/BA4GHgy8taoevFjdc3NztW7durH6N9RJJ8FRR8H110++bkmSNJ223x7WroXDD29SfZIzqmpu4fpJzGw9CDi/qi6sql8DJwOHLChzCPC+6nwL2CHJ7hNoe8sce6xBS5Kkleb667sMsJVNImztCVw08Hp9v26pZQBIclSSdUnWbdy4cQLdG+KnP21TryRJmm7LkAEmEbYyZN3Ce5ObU6ZbWbW2quaqam716tVjd26offZpU68kSZpuy5ABJhG21gN7D7zeC7hkC8psPa99bXffVpIkrRzbb99lgK1sEmHrO8D+SfZLcgfgUOCUBWVOAZ7ZfyvxIcBVVbVhAm1vmcMP7x6Q23ff7vU223S/d965+0luXh62fWusm4Y+TGu/7MN098s+THe/pqEP09ov+zDd/Rq3D/vu2/Th+E1ZNW4FVXVjkhcCnwe2Ad5dVeckObrffgJwKt03Ec8HrgeePW67Yzv88GU54JIkaWUZO2wBVNWpdIFqcN0JA8sFvGASbUmSJM0S/4K8JElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIbGCltJdkryxSTn9b93HFJm7yRfSXJuknOSvGicNiVJkmbJuDNbxwBfqqr9gS/1rxe6EXhJVf0h8BDgBUkOGLNdSZKkmTBu2DoEOLFfPhF42sICVbWhqr7bL18DnAvsOWa7kiRJM2HcsHW3qtoAXagCdt1U4SRrgPsBp4/ZriRJ0kxYtViBJP8M7DZk07FLaSjJnYCPAy+uqqs3Ue4o4CiAffbZZylNSJIkTZ1Fw1ZVPW7UtiSXJtm9qjYk2R24bES529MFrZOq6hOLtLcWWAswNzdXi/VPkiRpmo17G/EU4Mh++Ujg0wsLJAnwLuDcqnrzmO1JkiTNlHHD1huAg5KcBxzUvybJHklO7cs8DDgCeEySs/qfg8dsV5IkaSYsehtxU6rqcuCxQ9ZfAhzcL38dyDjtSJIkzSr/grwkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDY0VtpLslOSLSc7rf++4ibLbJDkzyWfGaVOSJGmWjDuzdQzwparaH/hS/3qUFwHnjtmeJEnSTBk3bB0CnNgvnwg8bVihJHsBTwLeOWZ7kiRJM2XcsHW3qtoA0P/edUS544CXAr9drMIkRyVZl2Tdxo0bx+yeJEnS8lq1WIEk/wzsNmTTsZvTQJInA5dV1RlJ/nix8lW1FlgLMDc3V5vThiRJ0rRaNGxV1eNGbUtyaZLdq2pDkt2By4YUexjw1CQHA9sBd0nygar68y3utSRJ0owY9zbiKcCR/fKRwKcXFqiql1fVXlW1BjgU+LJBS5IkrRTjhq03AAclOQ84qH9Nkj2SnDpu5yRJkmbdorcRN6WqLgceO2T9JcDBQ9Z/FfjqOG1KkiTNEv+CvCRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpoVTVcvdhpCQbgZ80bmYX4OeN25hWK3nssLLHv5LHDit7/I595VrJ499aY9+3qlYvXDnVYWtrSLKuquaWux/LYSWPHVb2+Ffy2GFlj9+xr8yxw8oe/3KP3duIkiRJDRm2JEmSGjJswdrl7sAyWsljh5U9/pU8dljZ43fsK9dKHv+yjn3FP7MlSZLUkjNbkiRJDRm2JEmSGprpsJVk7yRfSXJuknOSvKhfv1OSLyY5r/+948A+L09yfpIfJvmTgfUPSPLv/ba3JcmINofuvxwmNf4k2yf5bJIf9PW8YUR7a5LckOSs/ueErTPSoX2Z5Ln/ar9ufly7jmhzKs79BM/7nQfGfFaSnyc5bkh7U3Pe+/4safxJdu7LX5vk+AV1zdR1P6mxz+I13/dnkuf+Nn3db+Lcz9x1vwVjPyjJGf21fUaSxwzUtTzXfFXN7A+wO3D/fvnOwI+AA4A3Acf0648B3tgvHwD8G7AtsB9wAbBNv+3bwEOBAJ8DnjikvZH7z/L4ge2BR/dl7gD864jxrwG+t9znvcG5/yowt0h7U3PuJzn2BfWeATxyms/7Fo7/jsDDgaOB4xfUNVPX/aTGPovXfINzf1u/7keOfUG9U3/db8HY7wfs0S/fG7h4oK5lueZnemarqjZU1Xf75WuAc4E9gUOAE/tiJwJP65cPAU6uql9V1X8A5wMPSrI7cJeqOq26I/2+gX0GDd2/yeA2w6TGX1XXV9VX+np+DXwX2GurDWQLTGrsS2hyas59i7En2R/Yle5Dd6otdfxVdV1VfR345WA9s3jdT2rss3jNw+TGvwS3uXM/aFau+y0Y+5lVdUm//hxguyTbLuc1P9Nha1CSNXRp9nTgblW1AbqTRPdmgu7kXDSw2/p+3Z798sL1C43af9mNOf7BenYAngJ8aURT+yU5M8nXkjxiYgMYw4TG/p5+qvyVI6aVp/LcT+q8A4cBH+7/AzTM1J132OzxjzLT1/2YYx+sZwdm7JqHiY3/tnzdb46Zu+63YOx/BpxZVb9iGa/5VePsPC2S3An4OPDiqrp6xC1Y6KYNF6pNrN/c/ZfVBMY/X88q4EPA26rqwiFlNwD7VNXlSR4AfCrJvarq6vFGsOUmNPbDq+riJHfu6zqC7l88m7v/spjUee8dSjfuYabuvMOSxj+yiiHrZuK6n8DY5+uZuWseJjb+2/p1vzlm6rpf6tiT3At4I/D4+VVDim2Va37mZ7aS3J7u4J9UVZ/oV1/aTxfO3yq4rF+/Hth7YPe9gEv69XsNWb/QqP2XzYTGP28tcF5VHTesrX5K9fJ++Qy6+9j3nNBQlmxSY6+qi/vf1wAfZPh08VSd+0me9yT3AVb15/RWpu28w5LHP8pMXvcTGvu8mbrmYXLjXwHX/WJ1zdR1v9SxJ9kL+CTwzKq6oF+9bNf8TIetftr3XcC5VfXmgU2nAEf2y0cCnx5Yf2h/73Y/YH/g2/304zVJHtLX+cyBfQYN3X/iA9tMkxp/X9drgLsCL95Ee6uTbNMv373ff9i/hpub1NiTrEqyS1/n7YEnA98b0uTUnPtJnvfeYXSzG6Pam5rz3vdhqeMfahav+0mNva9rpq75vg8TGf8Kue4XMzPX/VLHnu7W+GeBl1fVN+YLL+s1X1PwTYMt/aH7pkUBZwNn9T8HAzvTPX9wXv97p4F9jqVL6D9k4FsIwBzdxXYBcDz87q/rPxV41WL7z/L46VJ70T10OF/PcxeOn+7e9zl039L4LvCU28DY70j3bZyz+7G9lZu/pTiV536S7/t+24XAHyxYN5XnfYzx/xj4BXAt3b9aD+jXz9R1P6mxM4PX/ITHv1Ku+6Hv+37bzFz3Sx078ArguoGyZwG79tuW5Zr3f9cjSZLU0EzfRpQkSZp2hi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLU0P8Ha0wgQvgqND4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10,4))\n",
    "plt.plot(t,0*t,'o:',color='red')\n",
    "plt.title('Illustration of discrete time range for h=%s'%(h))\n",
    "plt.ylim(-0.5,0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gjjlaXtObPAJ"
   },
   "source": [
    "So it looks like our time range works, starting and ending at the correct values and with the correct step size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GRQBKvheqTRj"
   },
   "source": [
    "## 1. Linear Population Equation\n",
    "\n",
    "### Tailoring our 4th Order Runge Kutta scheme\n",
    "To write the specific 4th Order Runge Kutta method for the linear population equation we need \n",
    "\n",
    "$$\n",
    "f(x,y)=0.1y.\n",
    "$$\n",
    "\n",
    "(where I've switched back to calling the dependent variable 'x' rather than time 't' for consistency with lecture notes).\n",
    "\n",
    "So lets go ahead and create a function which returns the right hand side of this equation for a given value of $x$ and $y$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "g36MCt9kqTRj"
   },
   "outputs": [],
   "source": [
    "def l_fun(x,y):\n",
    "    ftw = 0.1 * y\n",
    "    return ftw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2tri4OS9qTRm"
   },
   "source": [
    "We now need to evaluate each of the coefficients, beginning with $k_1$. In the lectures we manually calculated formulae to describe each coefficient and evaluate these in a table before combining.\n",
    "\n",
    "The first coefficient, $k_1$, would look like:\n",
    "\n",
    "$$\n",
    "k_1=f\\left(x_n,y_n\\right)=0.1 y_n.\n",
    "$$\n",
    "\n",
    "Subsequent coefficients are more complicated to evaluate mathematically, e.g. $k_2$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "k_2&=f\\left(x_n+\\frac{h}{2},y_n+\\frac{h}{2}k_1\\right)=0.1\\left(y_n+0.05k_1\\right) \\\\\n",
    "&=0.1\\left(y_n+0.05(0.1y_n)\\right)=0.1y_n(1+0.005) \\\\\n",
    "&=0.1y_n(1.005)=0.1005y_n.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This is a lot of mathematical effort to then plug into our calculators!\n",
    "\n",
    "The beauty of this method relies on making the computer to do the hard work for us! Now that we've defined a function to spit out values of $f(x,y)$, we can also evaluate *shifted values* of $x$ and $y$ as input into that function to evaluate $k_1-k_4$, and hence the next iteration of the solution, $y_{n+1}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "THwmJrL2qTRm"
   },
   "outputs": [],
   "source": [
    "yn = np.zeros(N + 1)\n",
    "yn[0] = 6.0\n",
    "## 4th Order Runge Kutta\n",
    "for n in range (0,N):\n",
    "    k1 = l_fun(t[n], yn[n])\n",
    "    k2 = l_fun(t[n] + 0.5 * h, yn[n] + 0.5 * h * k1)\n",
    "    k3 = l_fun(t[n] + 0.5 * h, yn[n] + 0.5 * h * k2)\n",
    "    k4 = l_fun(t[n] + h, yn[n] + h * k3)\n",
    "    yn[n + 1] =yn[n] + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5EQeGB9yFXyS"
   },
   "source": [
    "### Exact Solution \n",
    "From MAT201 mathematics, we should also know that this equation can be separated and integrated, using the initial conditions to evaluate any constants of integration.\n",
    "\n",
    "This simple separable first order ODE \n",
    "\n",
    "$$ \n",
    "y^{'}=0.1y, \\ \\ (2000 \\leq t \\leq 2020),\n",
    "$$\n",
    "\n",
    "with the initial condition,\n",
    "\n",
    "$$ \n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "has a simple exact (analytic) solution\n",
    "\n",
    "$$ \n",
    "y(t)=6e^{0.1(t-2000)},\n",
    "$$\n",
    "\n",
    "which you should be able to derive analytically yourself - see MAT201 notes for a refresher on this!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tLegBPPeqTRp"
   },
   "source": [
    "### Plotting Results\n",
    "\n",
    "We'll now plot the approximate Runge-Kutta solution against the exact solution over time, to see how they behave:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "executionInfo": {
     "elapsed": 6,
     "status": "ok",
     "timestamp": 1646995572546,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "4LdkLbOmqTRp",
    "outputId": "d87017e1-8112-4549-f5e8-5a6534039d65"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD4CAYAAAA0JjXXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwqUlEQVR4nO3deXxU1d3H8c+PEHZkCYsRhCCK7AQIuKCCIAiKCooLAqK00MdKxVqtqH1q20cste5FbXEBlKioqLhhRQUV3AgYkEVUZBFECGE3ICQ5zx9zJ5lJJgtkJplJvu/XK6/M3Dn33nMzXn/8zjn3HHPOISIiIuWrWkVXQEREpCpSABYREakACsAiIiIVQAFYRESkAigAi4iIVIDq5XmyJk2auKSkpPI8pYiISIVZtmzZTudc01CflWsATkpKIi0trTxPKSIiUmHMbFNRn6kJWkREpAIoAIuIiFSAUgdgM4szsy/N7E3v/V/MbKuZpXs/F0SumiIiIpXL0fQBTwLWAscFbHvQOXdfWSpw5MgRtmzZwqFDh8pyGDlGtWrVomXLlsTHx1d0VUREqpRSBWAzawlcCEwBbg5nBbZs2UL9+vVJSkrCzMJ5aCmBc47MzEy2bNlCmzZtKro6IiJVSmmboB8C/gjkFtg+0cxWmtnTZtYo1I5mNsHM0swsLSMjo9Dnhw4dIiEhQcG3ApgZCQkJan0QEUlNhaQkqFbN9zs1NeKnLDEAm9lQYIdzblmBjx4H2gLJwDbg/lD7O+emO+dSnHMpTZuGfBRKwbcC6W8vIlVeaipMmACbNoFzvt8TJkQ8CJcmA+4DXGxmG4EXgP5mNts5t905l+OcywWeAHpHsJ4iIiKRceedkJXFLgIacrOyfNsjqMQA7Jy73TnX0jmXBFwFfOCcG21miQHFhgOrIlTHiIuLiyM5OZnOnTtz0UUXsWfPngqpx6JFixg6dGje+z/96U+cf/75/PLLL0Xuc8899+S93rNnD4899lhE6ygiUqmkpvoyXuA1hvE2Q/I/27w5oqcuy3PA95rZV2a2EjgX+H2Y6lS8CLTT165dm/T0dFatWkXjxo159NFHy3zMspoyZQpLlizhtddeo2bNmkWWUwAWETlGXtNzjhcKxzGDQbyb/3mrVhE9/VEFYOfcIufcUO/1GOdcF+dcV+fcxc65bZGpYoByaKc/44wz2Lp1KwD9+vXLmzpz586d+OexnjlzJpdeeimDBw/mlFNO4Y9//GPe/k899RTt2rWjX79+jB8/nokTJwKQkZHBZZddRq9evejVqxdLliwpsg73338/b7/9Nm+88Qa1a9dm5syZeccBGDp0KIsWLWLy5MkcPHiQ5ORkRo0axeTJk1m/fj3JycnceuutHDhwgAEDBtCjRw+6dOnCvHnzwvZ3EhGJaampMHYsW7Ma0oWveJeBAFQnx/d5nTowZUpk6+CcK7efnj17uoLWrFkTvKFvX+dmzPC9PnzY9/7ZZ33vTzzROV/oDf5p0sT3eUaGr/zrr/veb9tW6Hyh1K1b1znnXHZ2thsxYoSbP3++V5W+bunSpd6hM1zr1q2dc87NmDHDtWnTxu3Zs8cdPHjQtWrVym3evNlt3brVtW7d2mVmZrrDhw+7s846y91www3OOedGjhzpPv74Y+ecc5s2bXLt27cvVI+FCxe6hg0bupNPPtnt3bs3b/uMGTPyjuOccxdeeKFbuHBhUN2dc27Dhg2uU6dOee+PHDmSd5yMjAzXtm1bl5ubW+i8hb4DEZHKbPZs5+rUcbto6Dqz0tVnr1tOcnBcmT07LKcC0lwRMbFcF2Mosy1bQm/PzCzTYf1Z5MaNG+nZsycDBw4scZ8BAwbQoEEDADp27MimTZvYuXMnffv2pXHjxgBcfvnlfPPNNwC89957rFmzJm//ffv2sX//furXrx903JNPPpndu3fz7rvvMmLEiDJdl3OOO+64g48++ohq1aqxdetWtm/fzvHHH1+m44qIxCwv8z2UU51LmMc6TuUdBtOd9PwyrVvDqFERr0r0BeBFi/Jfx8cHv2/VKq+zPIi/nb5Jk+DypQw0/j7gvXv3MnToUB599FFuvPFGqlevTm6u79Hngs/KBvbLxsXFkZ2dje8fO6Hl5uby6aefUrt27WLr0rx5c1JTUxkwYAAJCQmce+65QfUIVZeipKamkpGRwbJly4iPjycpKUnP/IpI1eXv881xjCKVxZzF84ykPwvzy5RH07MnthZjmDLF98cJFMY/VoMGDXjkkUe47777OHLkCElJSSxb5nv8+eWXXy5x/969e/Phhx+ye/dusrOzmTt3bt5ngwYNYtq0aXnv09PTizxOu3bteOWVVxg9ejTp6ekkJSWRnp5Obm4uP/zwA1988UVe2fj4eI4cOQJA/fr12b9/f95ne/fupVmzZsTHx7Nw4UI2hfrHi4hIVeBlvmRl4TAasZuHuIkreTG/TFwcTJ9eLtkvxFoAHjXK98dp3RrMfL/D/Mfq3r073bp144UXXuCWW27h8ccf58wzz2Tnzp0l7tuiRQvuuOMOTjvtNM477zw6duyY10z9yCOPkJaWRteuXenYsSP//ve/iz1Wr169mDFjBhdffDGJiYm0adOGLl26cMstt9CjR4+8chMmTKBr166MGjWKhIQE+vTpQ+fOnbn11lsZNWoUaWlppKSkkJqaSvv27cv2xxERiUX+Abw5OeyjPtXJ4QnGcyP/yi9Tpw7MmlVuwRfAims2DbeUlBTnH1Xst3btWjp06FBudYi0AwcOUK9ePbKzsxk+fDjjxo1j+PDhFV2tYlW270BEJI8/883J4V5u5TF+y6ecQSI/5ZeJi4tY8DWzZc65lFCfxVYGHAP+8pe/5E3q0aZNG4YNG1bRVRIRqZoCMt8n+DW3cS+n8xnN2JFfpgIyX7/oG4QV4+67r0yrM4qISLh4U0y+yOX8hv8whLd5hmuI868rVM59vgUpAxYRkcrHm2JyEX0ZzWz6sISXGUENfINWKzLz9VMAFhGRysXf9Ax0YwXXMpM3uIg6HPR9XsGZr58CsIiIVB7eoKtvslpwiJo0Yg/T+Q0N2ev7PAoyXz8FYBERqRy8zPebnJM4m4+ZwPTCZaIg8/VTACZ/OUL/z9SpU8N27PT0dN5+++2wHU9ERELwMt9NWU0YyAIcxp0UmKSpnKaYLC2NgiZ/KspISE9PJy0tjQsuuCAixxcRqfK8zHdLzvH05wP2cRwf0J9T+Sa/TDlOMVlayoCLsHfvXk499VTWrVsHwMiRI3niiScAuP7660lJSaFTp07cddddefssXbqUM888k27dutG7d2/27t3Ln//8Z+bMmUNycjJz5sypkGsREam0vMzXZWVxFS+QQVPeZVDw4gpRMuiqoKjLgPv1K7ztiivgt7+FrCwIlUhee63vZ+dOKLiAUODaDEXxr4bkd/vtt3PllVcybdo0rr32WiZNmsTu3bsZP348AFOmTKFx48bk5OQwYMAAVq5cSfv27bnyyiuZM2cOvXr1Yt++fdSpU4e//e1vpKWlBc0DLSIiYRAw0YYB/+Z/2Mdx9CJgxsU6daIy+MJRBGAziwPSgK3OuaFm1hiYAyQBG4ErnHO7I1HJSCuqCXrgwIG89NJL3HDDDaxYsSJv+4svvsj06dPJzs5m27ZtrFmzBjMjMTGRXr16AXDccceVV/VFRKoeL/PdlXMcLzCW63mczqwOLhOlma/f0WTAk4C1gD+yTAbed85NNbPJ3vvbylqh4jLWOnWK/7zgaoRllZuby9q1a6lduza7du2iZcuWbNiwgfvuu4+lS5fSqFEjrr32Wg4dOoRzDjML38lFRCQ0L/Pdk1OPQbzLKjpzHu/Rjm/zy0Rx5utXqj5gM2sJXAg8GbD5EmCW93oWMCysNYsCDz74IB06dOD5559n3LhxHDlyhH379lG3bl0aNGjA9u3bmT9/PgDt27fnxx9/ZOnSpQDs37+f7OzsQksEiohIGXiZ7/6sagxhPivpylwuCw6+UZ75+pU2A34I+CNQP2Bbc+fcNgDn3DYzaxbmupWbgn3AgwcPZty4cTz55JN88cUX1K9fn3POOYe7776bv/71r3Tv3p1OnTpx0kkn0adPHwBq1KjBnDlz+N3vfsfBgwepXbs27733Hueeey5Tp04lOTk5r29ZRESOgZf5/pxTkwt5i6X04iUu50ICHvWMgczXr8TlCM1sKHCBc+63ZtYPuMXrA97jnGsYUG63c65RiP0nABMAWrVq1bPgovBaCq/i6TsQkagXsKzguwzkIt7gWcZwBS/ll4ngsoLHqrjlCEuTAfcBLjazC4BawHFmNhvYbmaJXvabCIHrO+Vzzk0H33QkKSkp5bf4sIiIVA5e5uu80c6DWMB3nMyJbMkvE0OZr1+JfcDOududcy2dc0nAVcAHzrnRwOvAWK/YWGBexGopIiJVk5f5/pwFg3mHNxgKEBx8Y6TPt6CyTMQxFRhoZt8CA733x6SkZnCJHP3tRSRqeZnvgZxaDGE+73EeB6gXXCaKFlc4Wkc1EYdzbhGwyHudCQwoawVq1apFZmYmCQkJeoynnDnnyMzMpFatWhVdFRGRYP7Rzjm1uYC3+ZQzeI6ruZIX88vEaObrV+EzYbVs2ZItW7aQkZFR0VWpkmrVqkXLli0ruhoiIvm8zDcrpwaDeYfPOY3nGcnlvJxfJgb7fAuq8AAcHx9PmzZtKroaIiISDQJGO9cGerGUm3mAy3glv0yMZ75+FR6ARUREgKAZrnbRmJPYwEP8PrhMJch8/RSARUSk4gXM7TyYd8gkgbV0oAZH8stUkszXTwFYREQqlpf5bs9JYBDv8jXtmctlwcG3EmW+fgrAIiJScbzM94ecRM7jPbbQkjcZykDeyy9TyTJfPwVgERGpGAHr+d7O3/mJ4/kv53MWS/LLVMLM168sE3GIiIgcG/9o56wsAB7jt3zEOcHBt5Jmvn4KwCIiUr68zPfLnC5cxstkUZvj2E83VuaXieEZrkpLAVhERMqPl/l+ktWNc1lIGilk0DS4TCXPfP0UgEVEpHx4me/7OX0ZyAKasYPFnEVrNueXqQKZr58CsIiIRFZqKiQlwejRvJ3Vlwt5i5P4no84p1KsanSsFIBFRCRy/COdN20C4GS+YyALWEQ/jmd7frkqlPn6KQCLiEhkeP29LiuL+QzGAe34lje4mAR25ZerYpmvnwKwiIiEV2oqNGkCo0eTm5PLzTzABcznZUYULlsFM18/TcQhIiLh429yzsriCNW5jhmkMppJPMRlzA0u27o1TJlSJYMvKACLiEi4BCwl+DN1uJyXmM8F3MPtTGYq5i9XiWe3OholNkGbWS0z+8LMVpjZajP7q7f9L2a21czSvZ8LIl9dERGJSgHTSgIsoycf0J8n+DW3BwbfKtrfG0ppMuBfgP7OuQNmFg8sNrP53mcPOufui1z1REQk6gVkvgepRW0OcQ4fs562tODH/HLKfIOUmAE7nwPe23jvx0W0ViIiEhsCMt8VdOVU1vEqwwCCg29CgoJvAaUaBW1mcWaWDuwAFjjnPvc+mmhmK83saTNrVMS+E8wszczSMjIywlNrERGpeAELKrzLQM7mYxzGyXyXXyYuDmbPhp07FXwLKFUAds7lOOeSgZZAbzPrDDwOtAWSgW3A/UXsO905l+KcS2natGmoIiIiEmsCMt+ZjOVC3qING/iUM+jCKl+ZKvyIUWkc1XPAzrk9wCJgsHNuuxeYc4EngN7hr56IiESdgMz3U07nOmbSj0V8zNm0ZKuvjAZblag0o6CbmllD73Vt4DzgazNLDCg2HPz/5BERkUopYIIN/2jn0/mMVK7mLS7kOPb7yinzLZXSZMCJwEIzWwksxdcH/CZwr5l95W0/F/h9BOspIiIVyd/knJnJPupzOS+yki4YcDXPU4MjvnLKfEutxMeQnHMrge4hto+JSI1ERCS6BDxmtJUTuJC3WEVnLuUVuvJVfjk9ZnRUNBOWiIgULWCw1XK6czGvs5cGvMWFnM+7+eWU+R41BWAREQktIPNNoyd9+ZAm7GQJfZT5hoFWQxIRkWAhBlt1ZSUTmM4X9A4Ovppg45gpAIuISL6AwVaHied27mEnCdTgCA9yM83Z4SunCTbKTAFYRER8Ap7v3UkCA1nAVG7nTYYGl9NjRmGhPmAREQkabLWW9gzlTbbSgucYyUheyC+nwVZhowAsIlLVBQy2+ozTGMw71OQXFtGP0/k8v5wGW4WVmqBFRKqyAuv4nsx3nMNHfEHv4OCrwVZhpwAsIlJVeZnvwaxcpnIbR6hOEzJ5nUtozWZfGQ22ihgFYBGRqibgMaMtOcdzDh9xO1N5nwHB5TTYKqLUBywiUpX4m5yzsljCmVzGXLKowzwuZjD/zS+nwVYRpwxYRKSqCHjM6Hmu4lwWchz7+JzTuJg38ssp8y0XCsAiIpVdiJmtOrGai3iDzzmNDnydX1aZb7lRABYRqcwCZrbaxvE8yE0AdOUr5jKCRuzJL6vMt1wpAIuIVFYBTc4fcxY9WM6fuJsNJBUuq8eMyp0CsIhIZeRlvi4nhwf4fVB/bxs25pfTY0YVRgFYRKSyCch8x/MEf+ABLuZ1ltKLzqzOL6cm5wpV4mNIZlYL+Aio6ZV/2Tl3l5k1BuYAScBG4Arn3O7IVVVERIqVmgqTJkFmZt6mC3ibU1nHLdyHBZZNSICHH1bwrUCleQ74F6C/c+6AmcUDi81sPnAp8L5zbqqZTQYmA7dFsK4iIlKUgOd7X2IEe2jIeJ7kUl4NLhcXp6w3SpTYBO18Dnhv470fB1wCzPK2zwKGRaKCIiJSAq/J+XDWEW7mfq7gJVIZRW5wzqsm5yhTqj5gM4szs3RgB7DAOfc50Nw5tw3A+92siH0nmFmamaVlZGSEqdoiIhL4fO+mnBacw0c8yM1M5F+8yyCq4fLL6vneqFOqqSidczlAspk1BF41s86lPYFzbjowHSAlJcWVUFxEREojoMl5Nw3pyTKOEM9LjGAEc4PLahnBqHRUc0E75/aY2SJgMLDdzBKdc9vMLBFfdiwiIpHmNTm7nBwMaMQe/s7t9OcD2vJ9cFkNtopaJTZBm1lTL/PFzGoD5wFfA68DY71iY4F5EaqjiIhAUJPzhpwTOZNPWEg/AMbzZHDw1fO9Ua80GXAiMMvM4vAF7Bedc2+a2afAi2b2K2AzcHkE6ykiUrUFNDm/wnDG8TQAP1O3cFk1OceEEgOwc24l0D3E9kwouHikiIiEndfkfCinOrfyCNP4Hb34gjlcGTyrFajJOYZoJiwRkWhVYBWjF7iKafyOm3iQxZylKSVj3FENwhIRkXLin8s5K4sNtOEkNnANz3Aq6ziDz4LLqsk5JikDFhGJNl6T886s2gzjNXqyjJ9oTjVc4eCrVYxiljJgEZFoETCX8/v0ZwzPkkkC/+A2mhV80lNTSsY8ZcAiIhUtoK83N3MXtzGVgSygAXv5jNO5iYeDZ7XSlJKVggKwiEhF8j9e5K1gVA3Hj5zAeJ4gjRS6kx5cXk3OlYaaoEVEKorX15ubk8s0fkd/PqAzq5nJtcSRG1xWTc6VjjJgEZHyFtDk/ENOIoN4l0k8wkyuBSgcfNXkXCkpAxYRKU8BjxelMoqJTCOb6kxnPL/mycLlNbFGpaUALCJSXrwmZ3JySGUUY5hNHxYzi7FaRKEKUgAWEYm0gMeL9nIcDdjH5bxEFnX4FU8FNzmrr7fKUB+wiEikBPT17snMZhxP0Y0V7OU4anKYCTwRHHzV11ulKACLiERCwONFbzOEzqziGa7hap6jFocKl9fjRVWOArCISDilpkJSEowezcGsXK5lBhfyNo3YzWeczj3cSU0O55fXIgpVlgKwiEg4BK5ctGkTADX5hW0k8if+jzRSSGFZ8D5qcq7SFIBFRMoiMPBmZrKLRkzgP2yhBdVwzGcI/8efg7NeUJOzKACLiByzgH5eB8zlUjqxmhlcx8ecDRA8hzP4Aq+anIVSBGAzO9HMFprZWjNbbWaTvO1/MbOtZpbu/VwQ+eqKiEQJ/zO9WVn8QEuG8RojmEsi21hKL0byQnD51q0VeCVIaZ4Dzgb+4Jxbbmb1gWVmtsD77EHn3H2Rq56ISJQJeKbXbwp3soCB/JNbuImHqE5Ofvk6ddTULCGVGICdc9uAbd7r/Wa2FmgR6YqJiESVAoH3KzpTjVw6sYYp3Mlt/IM2bAzeR7NZSTGOqg/YzJKA7sDn3qaJZrbSzJ42s0ZF7DPBzNLMLC0jI6NstRURKW8FBlkdpBZ3MIUeLOcWfA2ACewKDr7q55VSKHUANrN6wFzgJufcPuBxoC2QjC9Dvj/Ufs656c65FOdcStOmTcteYxGR8lJgrd4POJeurOTv3MFoZjOb0cHl9UyvHIVSBWAzi8cXfFOdc68AOOe2O+dynHO5wBNA78hVU0SknAUMsgKYx8UM4AMcxnsMYAbjSGBXfnk90ytHqTSjoA14CljrnHsgYHtiQLHhwKrwV09EpJwFNDnn5uSygSQABvMO93IrX9GFAXwQvI+e6ZVjUJpR0H2AMcBXZpbubbsDGGlmyYADNgK/iUD9RETKR4FBVmn05AYeZRuJfE176nCQWynw0IcGWUkZlGYU9GLAQnz0dvirIyJSzgoE3l004k6m8B9+QzN2cB+3UJuDwfso8EoYaD1gEam6/IOsvH7eDSTRmy/YRWNu5BH+yl00YF9+ea3VK2GkACwiVZN/kFVODntoQEP2ksRGxvAsY5lFN1YGl9eEGhJmmgtaRKqWgEFWe3LqcSMPk8RGtnICBjzAHwoHXw2ykghQABaRqiEg8OZk7uY/TKAd3/AoNzCa2dTl58L7aEINiSA1QYtI5VZgkNVBanEmn5BOd87mIx5mEt1JD95Hg6ykHCgDFpHKKcQ6vQC1OcQQ5vMSI/iQvsHBVzNZSTlSABaRyqVA4N1HfSbzd1qyhS9JBuAe7mQEc4Ofr9RMVlLOFIBFpHIoEHhzqMZTjKMd3/APJnMFL9Kc7aH31SArqQDqAxaR2Ffged5cjH4sYjFncyZLeIOL6EVa4f3U1ysVSBmwiMS2gEUT1nMSDqiG40rm8BwjWcxZhYOvRjdLFFAAFpHYFNDkvDWnOb/mCdrxDW9wEQATeZSRvBDcz6vAK1FETdAiElsCHivaR33u5f94gJvJpjo38ghn8knhfdTULFFIAVhEYkOB53kd0I9FfEkPRvIcd/MnTmJD8D4KvBLFFIBFJLoFBF4HzOMShjCfmhzmHu6gCTtJYVnwPlo0QWKA+oBFJDoF9PG6zEwWcB6n8TnDeY1UfIF1MP8tHHz1PK/ECAVgEYkuBZ7nXcKZnMtCBrGA7TTnaa5jLLNC76vneSWGqAlaRKJDgT5e8PXz3sRD/MCJ/IuJjOcJanK48L7q65UYVGIGbGYnmtlCM1trZqvNbJK3vbGZLTCzb73fjSJfXRGpdApkvGvowGieZScJGPACV7Getkzk0cLBV48VSQwrTRN0NvAH51wH4HTgBjPrCEwG3nfOnQK8770XESmdAoF3PSdxDbPozCrmcQnL6QFAW76nLlnB+yrwSiVQYgB2zm1zzi33Xu8H1gItgEsgryNmFjAsQnUUkcqkQODNxbiex2jP17zE5fyB+9lAGwaxoPC+CrxSiRzVICwzSwK6A58DzZ1z28AXpIFmRewzwczSzCwtIyOjjNUVkZhVIPBm0hjwTRu5h4ZMYDrracs/+SNNyAzeV4FXKqFSB2AzqwfMBW5yzu0r7X7OuenOuRTnXErTpk2PpY4iEssKBN7vaMs4nuIEfmQd7QB4jqt5lImcwLbgfRV4pRIrVQA2s3h8wTfVOfeKt3m7mSV6nycCOyJTRRGJSQUC7zecwlhm0p6veZ6RXM/jNGQPQPB8zaDAK1VCiY8hmZkBTwFrnXMPBHz0OjAWmOr9nheRGopIbAnxONE+6tOD5eRSjUk8zK38k+NDrc2rx4mkCinNc8B9gDHAV2aW7m27A1/gfdHMfgVsBi6PSA1FJDYUCLwr6cJrDOPP/B/HsZ9URnEGn9KMEGNBFHilCioxADvnFhOihcgzILzVEZGYkpoKd94JmzaBGc45FnMWU5nM21xIPfYzllm0ZjOX8Hrh/RV4pQrTVJQicvQC+3c3bQLge5fEWSzmHD5mKb24mzvZTCtas7nw/urjFdFUlCJyFAo0M2cTxyZa05bvOZ6fOEwNpnED1zGDOhwsvL8yXpE8CsAiUrICgTeL2szgOu7jFuI5wlo6UIeDLKV36P0VeEUKURO0iBStwKNEGTThr/yZJDYykUc5gR95gJsxXOj91dQsUiQFYBEprEDg9YfXD+nLX/grvfmCjzibJZzFUN6imr+EeeM1W7dW4BUpgZqgRaTQaGacwwEf0J8HuJnT+Yz/5W6G8RpraU971hU+hpqZRY6KMmCRqizEaObDrjrPMIbufMl5vE8aKXkzVlUnp3DwVTOzyDFRBixSFYWYrcpvPE/wDGPpxCqeYhxX8xy1+KXwMZTxipSJMmCRqqRA3y5AOt34FU/yHW0BuImHeIfz+YoujGNG4eCrjFckLJQBi1QFBTLeI1TnVYbzL37HYs6mDj8zmHc4mfV0Jz30MZTxioSVArBIZRaiqfkI1enAWtZzMiexngf4Pdcyk0ZeP28hCrwiEaEALFKZhBjNDPA5vXmHwdzF34gnm4lM4xS+ZQjz8x8h8qtWDXJzfY8STZmiwCsSIQrAIpVBiEz3gKvDC1zFf/gNafSiPvuYwHQS+YmbeLjwMZTpipQrDcISiWUhBlUBLKYPJ/Aj43mSQ9RiGjewlRYk8lPhY2hQlUiFUAYsEosKZLwHqcVLXE599jOc1+jGCi7nJX7FU5zBp6HXE1XGK1KhFIBFYkERfbtrac9/+A3PcA27acwwXmU4r1GfAzzFr0MfS4FXJCqU2ARtZk+b2Q4zWxWw7S9mttXM0r2fCyJbTZEqKsRMVf7gexMP0pG1PMZvGcS7LKQfr3Bp0cdSU7NIVClNH/BMYHCI7Q8655K9n7fDWy2RKi7EYgifcAbjmc4OmgIwgPeZym1soSUvMJJ+fBjc1FzNu721MIJIVCqxCdo595GZJZVDXUSkQN/uNo7nWcbwNONYR3vqcoDLmMtg/stFvMlFvFn4GGpiFokJZRkFPdHMVnpN1I3CViORqiQ1FZKSfP261aoFjWbeSQKt2cRt3EtTMnia6/iJ4xnMf0MfS03MIjHlWAPw40BbIBnYBtxfVEEzm2BmaWaWlpGRcYynE6lkQvTtrnYduIV/Mo6nAGhCJo9yA+tox8ecw3XMpB4/Fz6WAq9ITDqmAOyc2+6cy3HO5QJPAL2LKTvdOZfinEtp2rTpsdZTJLYVkelu43ge4Pf0YBmdWc3DTOIA9cj1enPH8yTt+Db4WOrbFakUjikAm1liwNvhwKqiyopUaSEy3QOuDr9QA4CZXMsfeIDqZPMwN7KVFrzIlYWnh4T8TDcnxzcSeuNGBV6RGFbiICwzex7oBzQxsy3AXUA/M0sGHLAR+E3kqigSgwoMpsomjvc4j2cZw2sMYwbXcQUv8Wue5FJe4VS+KfpYGlQlUimVZhT0yBCbn4pAXURiVxETZRyiJpOZyvOMZAfNacQuruEZOrAWgKbspCk7Qx9TgVekUtNMWCJlEWIRhDWuPd/QjmHMoya/sIh+nM3HjGY2Q5hPTQ6HPpZWIRKpUhSARY5GEZnut5zMHK5kDleyii40JpOhvEl1clhGT+LILfqYynRFqiSthiRSksARzGPGFJoS8l5upR3f8r/cTQP28gi/YxWdqU4OQOHgq1HMIoIyYJHCishyAX5wLXiRK3iRK/gnt3IOH3M+/yWeI1zOS7Rka9HHVaYrIgEUgEX8QvTn4hwHqcW/+R9eZgSf0AeAHizjELUA6MZKurGy6OMq8IpICGqClqqriMkxHLCGDrzHAADiOcIU7uQA9bibO/mGU1hGCoNYEPq4amIWkVJQBixVT4hM1znHMnryCpfyCpeyjva05Tu+5RSqk8M6TiWBXcUfV5muiBwFZcBS+RWR6eYE/Oc/kWn0Io17+SMn8gOPcT0fcU7e8n5Bwde8rXFxvt/KdEXkGCgDlsqnmEFUB1wdFjCQN7iINxnKR5xDe9ZxNc/Rmy8YypvFZ7rKckUkTBSApfIoYhAVwHpOYiLT+ID+HKYmDdnNEOaT62XBffiEPnxS+JiaHENEIkQBWGJXEZluLkYaKbzBRbTjG8YwmybsZBOtmcg0LuIN+rCEeLKLPrYyXRGJMAVgiR3FNC3jHG9xAa8ynLe4kJ9IpBo53MCjjGE2DdjHGjoVfWxluiJSzhSAJXoVE3Cdc6ylA1/SnVE8B8AD3EwaKQzmHS7iDYYwXyOXRSRqKQBL9CmiL3c/9XifAcxnCO8wmM20pho5DOVNGrCPWYylGTuowZGij61MV0SihB5DkopXzIQYK+nCXo4D4GnGMZzXeJ6R9GQZ/2ECG2hDA/YB0JKthYNvwUkxtJi9iEQJZcBS/oppWt7umvI+A3iP8/gv5/MjLXiW0YwmlSt4kW6s4Ew+KT7LBTUti0jUUwCWyCvhudw9NKQlW/mRRFrwIwCN2MUA3mcI8xnoTfmYyE8k8lPoc6hpWURiTIkB2MyeBoYCO5xznb1tjYE5QBKwEbjCObc7ctWUmFJMwM121VhKL97jPBYwkE85g0t5hTlcxQlsYxo3cBqf050vi19DVwFXRGJcafqAZwKDC2ybDLzvnDsFeN97L1VVqD5cb81c5xxbaJFXtC8fciafchd/JYs6/IH7mci0vM9v4DFSWFbyGrrqyxWRGFdiBuyc+8jMkgpsvgTo572eBSwCbgtnxSSKlfB40De0YxH9+JC+LKIf+6nPLhoTTza/50F+z4Ocy8KSHxEC9eWKSKV1rH3AzZ1z2wCcc9vMrFlRBc1sAjABoFWrVsd4OqlQpXgetw0bqM0hpjKZO/g7AIn8SF8+pD8fkE114slmBHOLP5ealkWkioj4ICzn3HRgOkBKSoorobhUNH+w3bwZGjeGQ4fg55/zPs51sIZOeRnuh/Qlg2b8l0EMYgEX8QZNyaAvH3Iy3+WtJlQkBVwRqaKONQBvN7NEL/tNBHaEs1JSjoqb3jEzk2ziWEEP6rOfdnxLOsn0ZDkArdjEYN6hH4voxgoAOrOazqwu+nwKuCIiwLEH4NeBscBU7/e8sNVIIquE+ZQdsICBLKEPS+jDZ5zOz9RjIv/iX9xIN1Ywk7H05UOS2FTy+RRwRURCMueKbxU2s+fxDbhqAmwH7gJeA14EWgGbgcudcyWOqElJSXFpaWllq7EcneICLrCZE1nMWRymBtcyCwe0ZAs/cTxdWclZLKYPSzibj/Oe0S2WAq6ISB4zW+acSwn1WWlGQY8s4qMBZaqVhFdgoI2L8z2mEyLDBZjDFbzKcJbQhy2cCEBHVnMtszBgPkNowwbqc6Dk8yrgiogcE82EFauKy2xzcnDAD64ln3Man3MaK+nKfIYQRy4LOZcl9MnLbvuwhC58lbd714DXhSjgioiEhQJwrCih7/YAdanBYWpwhOcYyS3cxzZOAKAGv9CD5eykCc3ZwcNM4t9cX/I569aFWrVg1y5o1UoBV0QkjBSAo00pmpJznLGGTnnZ7eecxmo68R7ncS6LaMFW+vNB3qfdWEFNDuftH/g6iLJbEZFyowBc0UpoSs6hGuvcqSyjJx1ZQ0+Ws5we9GYp4Fu04DQ+51Je4UR+AKAvH9GXj0o+twKuiEiFUQAuL0cxSOoXanAb/2AZPfmS7vxMPQBu4Z/0ZDldWckzjOF0PivdZBd+CrgiIlFDATjcShNovcz2G9eO5fRgGT1JI4WOrOHfXE8NDjOPSziBHxnH06SQRk+W0Z6vAV8T8hhmF10Hf6D1n18BV0Qk6igAl1UJTcjgW/N2JV3ZQTOGeXOWnMGnLKU3ALXJIpl0WnsTWxjwPScpsxURqcQUgEurFJmtcy4vaM7hCl5mBOkks562OKrRiF1cQgIG3MRDZFOdZNLpyBqqkxN0umKDrwKuiEjMUwAuqJRNyL9Qg69de1bQjXSSWUE3vqILG2hDXbJIJ5l0kunGCsYyi2TS8+ZLBria50uui5qSRUQqraoZgAuu+AOQmVlkX+1615ZVdGYVnfk1T3IC25jOBG7kXwDU4iBd+IpLmMfP1KUuWdzDHfydO46uXspsRUSqjModgEuTzWZmAuTNHFWXn0lgF2n0ZALTWUsHDlEbACOX3nzBCWzjAt6mKRl0YwWn8O2xNSErsxURqbIqRwAuZbMxkLdtP/WYwXV5me0qOrOf43iU3/JbHqcRu2nGDvrzQV6JDqylLlkAtOV72vJ98fVSoBURkSLEZgAuxchj5xzrOJW1dGAtHfia9qylA5cwjz8xBYBJPEJjMunCV1zDM3RmFf35APAF2HcYcnT1UhOyiIiUUuwF4NRUmDABsnyZ6H5XNyjQNmc7k3gEgNP5jL00BKAlP9CBtRzPTwDU5wDbaUZTMkr/uI+fMlsRESmj2AvAd96ZF3wH8i7vMTDvo+ocYRivMYlHMOA5rqYZOziVdSGX1mtGRvHnUqAVEZEIib0AvHlz3stLmEd/PqA9X9OBtbRlPfFk531+AfNLPp5/xZ/MTAVaEREpN7EXgFu18vX9AhN5tPT7KZsVEZEoUq0sO5vZRjP7yszSzSwtXJUq1pQpUKdO0Z9X8y4pLs73u3VrmD3bF3Sdg+xs3++NGxV8RUSkwoQjAz7XObczDMcpHX/QLPjYkTJaERGJIbHXBA2+IKtAKyIiMaxMTdD4JpB618yWmdmEUAXMbIKZpZlZWkZGCaOORUREqoiyBuA+zrkewBDgBjM7p2AB59x051yKcy6ladOmZTydiIhI5VCmAOyc+9H7vQN4FbwFbkVERKRYxxyAzayumdX3vwYGAavCVTEREZHKrCyDsJoDr5qZ/zjPOefeCUutREREKjlzgQsZRPpkZhnApjAesglQfo9ARZauJfpUlusAXUu0qizXUlmuA8J/La2dcyEHQJVrAA43M0tzzqVUdD3CQdcSfSrLdYCuJVpVlmupLNcB5XstZR0FLSIiIsdAAVhERKQCxHoAnl7RFQgjXUv0qSzXAbqWaFVZrqWyXAeU47XEdB+wiIhIrIr1DFhERCQmKQCLiIhUgAoLwGZ2opktNLO1ZrbazCZ52xub2QIz+9b73Shgn9vN7DszW2dm5wds7+mtS/ydmT1i3uwgIc4Zcv9ouRYzq2Nmb5nZ195xphZxviQzO+itw5xuZv+Otmvxti/ytvnr2ayIc0b791I/4BrSzWynmT0U4nwR+V6O9jrMLMErf8DMphU4VkzdK0VdSyzeKyV8LzF1rxTzvcTavTLQfIsJfeX97h9wrMjfK865CvkBEoEe3uv6wDdAR+BeYLK3fTLwD+91R2AFUBNoA6wH4rzPvgDOAAyYDwwJcb4i94+WawHq4FtfGaAG8HER15IErIqB72URkFLC+aL+ewlx3GXAOeX1vRzDddQFzgL+B5hW4Fixdq+EvJYYvVeK+15i7V4p8lpi7F7pDpzgve4MbC3Pe6XCMmDn3Dbn3HLv9X5gLdACuASY5RWbBQzzXl8CvOCc+8U5twH4DuhtZonAcc65T53vL/JMwD6BQu4fTdfinMtyzi30jnMYWA60DEcdSytc13IUp4z67yXwmGZ2CtAM3//wy8XRXodz7mfn3GLgUOBxYvFeKepaYvFeKepajkLUfy+BYuRe+dJ5iwoBq4FaZlazvO6VqOgDNrMkfP8S+Rxo7pzbBr4/Jr4vEHx/xB8CdtvibWvhvS64vaCi9g+rMl5L4HEaAhcB7xdxqjZm9qWZfWhmZ4ftAoLrkETZr2WG18T0v0U04cTU9wKMBOZ4N2UoEf1eSnkdRYnFe6U0x2lIbNwrJYmle6U0Yu1euQz40jn3C+V0r5RlMYawMLN6wFzgJufcviKa2cHXDFCQK2Z7afcPmzBci/841YHngUecc9+HKLsNaOWcyzSznsBrZtbJObevbFcQUMHwXMso59xW862aNRcYg+9fkqXdPyzC9b14rsJ3HaFE9Hs5iuso8hAhtkX7vVLScWLpXilOrN0rpREz94qZdQL+gW9VPyine6VCM2Azi8f3R0p1zr3ibd7upf/+JrMd3vYtwIkBu7cEfvS2twyxvaCi9g+LMF2L33TgW+fcQ6HO5TV3ZHqvl+Hrd2gXpksJ27U457Z6v/cDzxG6aSZmvhcz6wZU9/7mhUTyeznK6yhKLN4rJYmle6VIMXivlHSsmLlXzKwlvvXsr3HOrfc2l8u9UpGjoA14CljrnHsg4KPXgbHe67HAvIDtV3nt822AU4AvvOaE/WZ2unfMawL2CRRy/2i6Fu9YdwMNgJuKOV9TM4vzXp/k7R/qX/8Vdi1mVt3MmnjHjAeGEnq96Jj4Xjwj8WVbRZ0vIt/LMVxHSDF6rxR3rFi7V4o6TizeKyWJiXvFfN0XbwG3O+eW+AuX273iIjA6sDQ/+EbQOWAlkO79XAAk4OvL+db73Thgnzvx/UtpHQEj0oAUfP/BrgemQd4MXxcDfytp/2i5Fnz/enL4Bg74j/PrgteCr69iNb7Rd8uBi6LwWuriGwG50qvrw+SPjo6p7yXgs++B9gW2Rfx7Ocbr2AjsAg7g+1d6xxi+VwpdC7F7r4S6lli9V0L+NxZL9wrwJ+DngLLpQLPyulc0FaWIiEgFiIpR0CIiIlWNArCIiEgFUAAWERGpAArAIiIiFUABWEREpAIoAIuIiFQABWAREZEK8P9sD5I61BWNTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_true = 6 * np.exp(0.1 * (t - 2000))\n",
    "fig = plt.figure(figsize=(8,4))\n",
    "plt.plot(t,yn,'ro:',label='Runge Kutta')\n",
    "plt.plot(t,y_true,'b--',label='Exact')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "InLe737QF4H4"
   },
   "source": [
    "As one might expect, because of the fourth order accuracy of the RK4 algorithm, both solutions are incredibly well matched when viewed this way.\n",
    "\n",
    "Rather than zoom in to any part of the graph, it can be easier to assess how accurate our scheme is using a table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SVzOXpE0qTRs"
   },
   "source": [
    "### Tablulating our results\n",
    "The table below shows the time, the Runge Kutta numerical approximation, $y$,  the exact solution, $y_{\\rm{true}}$, and the exact error $|y_{\\rm{true}}(t_n)-y_n|$ for the linear population equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 363
    },
    "executionInfo": {
     "elapsed": 9,
     "status": "ok",
     "timestamp": 1646995514778,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "QaOg9JvvqTRs",
    "outputId": "a111aa7c-b4f2-4b52-dde0-cb07a7f67365"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time t_n</th>\n",
       "      <th>Runge Kutta (y_n)</th>\n",
       "      <th>Exact (y_true)</th>\n",
       "      <th>Exact Error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2000.0</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2000.1</td>\n",
       "      <td>6.060301</td>\n",
       "      <td>6.060301</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2000.2</td>\n",
       "      <td>6.121208</td>\n",
       "      <td>6.121208</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2000.3</td>\n",
       "      <td>6.182727</td>\n",
       "      <td>6.182727</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2000.4</td>\n",
       "      <td>6.244865</td>\n",
       "      <td>6.244865</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2000.5</td>\n",
       "      <td>6.307627</td>\n",
       "      <td>6.307627</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2000.6</td>\n",
       "      <td>6.371019</td>\n",
       "      <td>6.371019</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2000.7</td>\n",
       "      <td>6.435049</td>\n",
       "      <td>6.435049</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2000.8</td>\n",
       "      <td>6.499722</td>\n",
       "      <td>6.499722</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2000.9</td>\n",
       "      <td>6.565046</td>\n",
       "      <td>6.565046</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   time t_n  Runge Kutta (y_n)  Exact (y_true)  Exact Error\n",
       "0    2000.0           6.000000        6.000000          0.0\n",
       "1    2000.1           6.060301        6.060301          0.0\n",
       "2    2000.2           6.121208        6.121208          0.0\n",
       "3    2000.3           6.182727        6.182727          0.0\n",
       "4    2000.4           6.244865        6.244865          0.0\n",
       "5    2000.5           6.307627        6.307627          0.0\n",
       "6    2000.6           6.371019        6.371019          0.0\n",
       "7    2000.7           6.435049        6.435049          0.0\n",
       "8    2000.8           6.499722        6.499722          0.0\n",
       "9    2000.9           6.565046        6.565046          0.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'time t_n': t[0:10],    'Runge Kutta (y_n)':yn[0:10],'Exact (y_true)':y_true[0:10],'Exact Error':np.abs(np.round(y_true[0:10]-yn[0:10],10))}\n",
    "df = pd.DataFrame(data=d)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wdkSen2QGzWW"
   },
   "source": [
    "You'll hopefully note that the error in the final column is zero. This is because the RK4 algorithm is highly accurate and we've used a very complicated algorithm for a very simple ODE (we've taken a sledgehammer to a problem that really only needed a tiny push!).\n",
    "\n",
    "Lets now explore two other Differential Equations where we will slowly crank up the complexity to see if that makes any noticeable difference to our iterative RK4 answers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fxlfUcFyqTRv"
   },
   "source": [
    "## 2. Non-Linear Differential Equation\n",
    "Our second equation will be a non-linear differential equation (effectively still an ordinary differential equation but with a higher order power of the dependent variable $y$).\n",
    "\n",
    "Since our first equation was an equation describing population evolution over time, we'll choose another differential equation describing population behaviour:\n",
    "\n",
    "$$ \n",
    "y^{'}=0.2y-0.01y^2, \\ \\ (2000 \\leq t \\leq 2020), \n",
    "$$\n",
    "\n",
    "with the same initial condition as before, namely that\n",
    "\n",
    "$$\n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "### Implementing RK4 \n",
    "\n",
    "To write the specific 4th Order Runge Kutta method we need\n",
    "\n",
    "$$\n",
    "f_2(x,y)=0.2y-0.01y^2.\n",
    "$$\n",
    "\n",
    "\n",
    "Again we could, if we choose to, calculate $k_1, k_2, k_3, k_4$ by hand, e.g.\n",
    "\n",
    "$$\n",
    "k_1=f_2(x_n,y_n)=0.2y_n-0.01y_n^2.\n",
    "$$\n",
    "\n",
    "However, we can again take advantage of the computer, and define a function to return the right hand side of the population equation, and feed in the correct values of $t$ and $y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "4Pbq4CxMqTRw"
   },
   "outputs": [],
   "source": [
    "def nl_fun(x,y):\n",
    "    ftw = 0.2 * y - 0.01 * y * y\n",
    "    return ftw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "iJjUCyoHqTRy"
   },
   "outputs": [],
   "source": [
    "yn2 = np.zeros(N + 1)\n",
    "yn2[0] = 6.0\n",
    "## 4th Order Runge Kutta\n",
    "for n in range (0, N):\n",
    "    k1 = nl_fun(t[n],yn2[n])\n",
    "    k2 = nl_fun(t[n] + 0.5 * h , yn2[n] + 0.5 * h * k1)\n",
    "    k3 = nl_fun(t[n] + 0.5 * h , yn2[n] + 0.5 * h * k2)\n",
    "    k4 = nl_fun(t[n] + h, yn2[n] + h * k3)\n",
    "    yn2[n+1] = yn2[n] + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nSVhMIVhqTR0"
   },
   "source": [
    "### Results and Graph\n",
    "\n",
    "While this ODE is more complicated, it also has an analytical solution; by again separating and integrating and making use of the standard integral result that \n",
    "\n",
    "$$\n",
    "\\int{\\frac{\\rm{d}x}{x\\left(a+bx\\right)}}=\\frac{1}{a}\\ln{\\left|\\frac{a+bx}{x}\\right|}+c,\n",
    "$$\n",
    "\n",
    "we can therefore find:\n",
    "\n",
    "$$\n",
    "y(t)=\\left({5e^{\\dfrac{-(t+c)}{5}}+0.05}\\right)^{-1}, \n",
    "$$\n",
    "\n",
    "with a contant of integration \n",
    "\n",
    "$$\n",
    "c=-5\\ln{\\left(\\frac{0.14\\exp{(400)}}{6}\\right)}.\n",
    "$$\n",
    "\n",
    "The plot below again shows the RK4 numerical approximation, $y_n$ (circles) and the analytical solution to the non-linear population equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "executionInfo": {
     "elapsed": 596,
     "status": "ok",
     "timestamp": 1646995562613,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "VGO634VDqTR0",
    "outputId": "ccf19c8e-8fb1-4516-894e-90aca19febb7"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD4CAYAAAA0JjXXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAucUlEQVR4nO3deXwURf7/8Vfl4JRDAioIJEFQ7gQJqKCAIIjIobuKoiLHrniwgq7+XAF3PVZc5YsHLLosq4Bi5BAvdhUVFAQRBaIBETxQblBCuAmBHPX7YyZJz2RCQjLJHHk/H488Mump7q7K0HxSn6quNtZaREREpGJFBLoCIiIilZECsIiISAAoAIuIiASAArCIiEgAKACLiIgEQFRFnqx+/fo2Li6uIk8pIiISMCkpKfuttQ18vVehATguLo5169ZV5ClFREQCxhizvaj3lIIWEREJAAVgERGRAFAAFhERCYAKHQP2JSsri127dpGZmRnoqlRK1apVo3HjxkRHRwe6KiIilUrAA/CuXbuoVasWcXFxGGMCXZ1KxVpLeno6u3btIj4+PtDVERGpVAKegs7MzCQmJkbBNwCMMcTExCj7ICKVV3IyxMWBMRAV5foeF+faXs4CHoABBd8A0u9eRMKarwCb9z0iAm67Dba77xTKyXF9374dRo0q9yAcFAFYRESkTHwF2qICrPu7tZY06rOWJBbye57jfvIf0JuRARMmlGuVFYCByMhIEhMTadu2LQMGDODQoUMBqcfy5cvp379//s+PPPIIV199NSdPnixyn6eeeir/9aFDh3jppZfKtY4iIgFRip6stZYDnE0KF/MWv+NZ/sw6OgLwNR2oxVHOIY3OrOVGFvIAz7GPcwrOuWNHuTYp9AJw3ocQEeG3PH316tVJTU1l48aN1KtXjxdffLHMxyyriRMnsmrVKt59912qVq1aZDkFYBEJK2fYkz1FND/YFvyPa3mBsSylFwB7OY86HCaGAySRwg28xYM8y6f0BKAJO7mD/zCFMbzHQFJJ4CB1OYd9BXVp2rRcmxrwWdBnJDnZlZfPyHD9nJenB7j1Vr+c4rLLLmPDhg0A9OjRg8mTJ5OUlMT+/ftJSkpi27ZtzJ49m0WLFpGRkcHPP//M9ddfz6RJkwB45ZVXeOaZZ2jUqBEtWrSgatWqTJs2jbS0NO666y52uP+ieuGFF+jatavPOjz77LN88MEHfPTRR1SvXp3Zs2ezbt06pk2bBkD//v158MEH+fDDDzlx4gSJiYm0adOGnJwcfv75ZxITE+nduzePPvoogwYN4uDBg2RlZfHkk08yaNAgv/yeRETKJDnZleLdvh0iI11B1RiwtqCMO9CeslFsJZ4tNOcsjtGdFeQQQUu+5xeakUtk/i73MpWr+IQGpDGSmcSxzeOrLocBaMB+nufPRdevRg2YOLFcmp7PWlthXx07drTeNm3a5Lmhe3drZ81yvT51yvXznDmun5s0sdb18Xh+1a/vej8tzVV+0SLXz3v3FjqfLzVr1rTWWpudnW1vuOEGu3jxYndVutu1a9e6D51mY2NjrbXWzpo1y8bHx9tDhw7ZEydO2KZNm9odO3bY3bt329jYWJuenm5PnTplL7/8cjt69GhrrbVDhgyxK1eutNZau337dtuyZctC9Vi2bJmtW7eubd68uT18+HD+9lmzZuUfx1prr732Wrts2TKPultr7datW22bNm3yf87Kyso/Tlpamr3gggtsbm5uofMW+gxERPzh9detjY11/T8dGen6HhNjbc2ahf4fzwW7g8Z2Pe3yt93Ga7YZW2wE2flF+/G//Pf/xFT7Vx63r3Gb/YJL7T7q21xfMaKkXxERru+xsa66+wGwzhYRE0OrB7xrl+/t6ellOmxeL3Lbtm107NiR3r17F7tPr169qFOnDgCtW7dm+/bt7N+/n+7du1OvXj0AbrzxRn788UcAli5dyqZNm/L3P3LkCEePHqVWrVoex23evDkHDx7k448/5oYbbihTu6y1jB8/nhUrVhAREcHu3bv57bffOO+888p0XBERD2fQm81NP0CEe6rTawzlE3qxmVZ8T0uOUpt2bGADCQBU5wSdWMstvEELfqI5W7iQH/MP+U/GnFk9IyIgN7egjnnfY2NdvV0/ZVJLKvgC8PLlBa+joz1/btq0YAzAKS9PX7++Z/kSBpq8MeDDhw/Tv39/XnzxRcaMGUNUVBS5ubkAhe6VdY7LRkZGkp2djXX+Y/OSm5vL6tWrqV69+mnrcu6555KcnEyvXr2IiYnhyiuv9KiHr7oUJTk5mbS0NFJSUoiOjiYuLk73/IpI6Z1BoMVadtCENXTmW9qxmVZsojXpxLCHRhjgY/qwjCtpzSaGM5tWbKYd3+YfagZ3lq6e3oE2QAG2OKE1CWviRFde3smPefo6deowdepUJk+eTFZWFnFxcaSkpACwcOHCYvfv3Lkzn332GQcPHiQ7O5u33nor/70+ffrkj+ECpKamFnmcCy+8kLfffpvbbruN1NRU4uLiSE1NJTc3l507d7JmzZr8stHR0WRlZQFQq1Ytjh49mv/e4cOHOeecc4iOjmbZsmVs9/XHi4hIUZwTooqaCGUt+4lhGT2Ywhj+wMscx/X/9IuM5kYW8iSP8DUXE89WhjKHLFxL385mOLtpzBL6MJWx3M10LmdV8fWKcIeuyEjP77Gx8PrrrrpZC9nZru/btgVd8IVg7AGfTt4vcMIE1/Twpk39/ldNhw4dSEhIYN68eTz44IMMHjyYOXPm0LNnz2L3Pf/88xk/fjyXXHIJjRo1onXr1vlp6qlTpzJ69Gjat29PdnY23bp1Y/r06UUeq1OnTsyaNYuBAweybNky4uPjadeuHW3btuXiiy/OLzdq1Cjat2/PxRdfTHJyMl27dqVt27Zcc801/OUvf2HAgAEkJSWRmJhIy5Yty/4LEpHwU5KerbVkE8n3tKQJO6nDEd7hOkbzIntplF+sAft4iElcxI+MYgaDWUBrNlGdwtm3KHJOX68Q6cmWljld2tTfkpKS7Lp16zy2bd68mVatWlVYHcrbsWPHOOuss8jOzub6669n5MiRXH/99YGu1mmF22cgIsVwBlzvQOtwgLN5j0Gk0JGvuZhUEjlBDd5lEINYxFqSeJHRtONb2rOB9mzgXOdtPCUVxoHWGJNirU3y9V5o9YBDwGOPPcbSpUvJzMykT58+XHfddYGukohUdqcLuNaSRRTf0SY/0PZmCdfxHmk0YCSzOIujdOAb7uTfdCSFS/gKgE6sYzYjSlaHvCAbE+P6+cCBcslihhIFYD+bPHlyoKsgIpXdaQKutZYMalCTDE4RTXc+4xs6cJJqANTiCHFsA6AFP7GZllzIj/kzl4sVxr1Zf1MAFhEJVSUYuz1uq5NCR1ZzGV9yKV9yKZ1YyyIGUYUsmvELXVlFR1LoSArN2ZIfbCOwtOQH3+dWoC0zBWARkVBxulRyTg4W2GIv4Ecu5Fo+AKAvH/I5VwDQgh/pzRKuYmnBIbmtZOfOC7gKtH6jACwiEqyKGbsF2MIFfEhfVtCNlVzBrzSkCic5Qm2qcorxPEUuEVzCV9SnBIsWqWdbYRSARUSCTXIyjB3rucqfe7JUCh1ZyRX8kZc5m0O8yY2M5x80YQc9+ZRurKArq4jGtT7ANXx4+nOpZxswxQZgY8xMoD+wz1rb1r0tEZgOVAOygXustWuKPEiQi4yMpF27dvk/33zzzTz88MN+OXZqaip79uyhX79+fjmeiISh0/R099CQ2QznU3qymsvIoCYACaynD0sYyUxu4Q1iKeGj8xRwg0ZJesCzgWnAa45tk4DHrbWLjTH93D/38HvtKkjeUpTlITU1lXXr1ikAi0iBIgKuBb6zrfmEXiSwnh58xkHOZgJP0Z71/JGX6cYKLufz/Ptti73vVgE3aBW7FKW1dgVwwHszUNv9ug6wx8/1CrjDhw9z0UUX8cMPrhmAQ4YM4T//+Q8Ad999N0lJSbRp04ZHH300f5+1a9fSpUsXEhIS6Ny5M4cPH+Zvf/sb8+fPJzExkfnz5wekLSISYKdZ0jHXwsv8gSG8wXn8Sjs2ch9TeJ9rAWjNJvbRgPUkMoX7+D1v+w663sszei/LGKTLMVZmpR0Dvg/4yBgzGVcQ71JUQWPMKGAUQNMSPNy4R4/C2wYPhnvucT0G2FdHcvhw19f+/eD9ACHnsxmKkvc0pDzjxo3jpptuYtq0aQwfPpyxY8dy8OBB7rjjDgAmTpxIvXr1yMnJoVevXmzYsIGWLVty0003MX/+fDp16sSRI0eoUaMGTzzxhMezfEUkzOX1bnfsgHr1IDMTjh/Pf/uErcpyerCPcxjGa0RgeZqHOU5NerOEXnxCTz7NTykbXM+uLUQ925BX2gB8N3C/tfYtY8xg4BXgKl8FrbUzgBngWoqylOcrV0WloHv37s2bb77J6NGjWb9+ff72BQsWMGPGDLKzs9m7dy+bNm3CGEPDhg3p1KkTALVr1y50PBEJY74mTrlf/0wz/ssAPqQvn9GdTKrTmJ3czmsY4Au60IA0zOmOr4AbdkobgIcBY92v3wRe9k91Tt9jrVHj9O97P42wrHJzc9m8eTPVq1fnwIEDNG7cmK1btzJ58mTWrl3L2WefzfDhw8nMzMRaizGnvXxEJJycZuLUcWrwCb24hsVEk82LjOZ5/sxFfM9dTKcvH9KNFfkB9xzSij5PTAxMmaKAG4ZK+zjCPUB39+uewE/+qU5wef7552nVqhVz585l5MiRZGVlceTIEWrWrEmdOnX47bffWLx4MQAtW7Zkz549rF27FoCjR4+SnZ1d6BGBIhLikpNdf+07H81nLXs5jxncQX/+SwzpDGIRX3EJAPfxAr8Qz/e04nn+zNV87PPpQEDBWG7eGO7+/Qq+YaoktyHNxTXDub4xZhfwKHAHMMUYEwVk4h7jDVXeY8B9+/Zl5MiRvPzyy6xZs4ZatWrRrVs3nnzySR5//HE6dOhAmzZtaNasGV27dgWgSpUqzJ8/n3vvvZcTJ05QvXp1li5dypVXXsnTTz9NYmJi/tiyiISQ08xYPklVqnGSNXTiElx3YsbzC3cxnQH8l87ubU3ZWfTxlVqutPQ4QtFnIOKLjzHdLKJYQTcWMZBFDORG3mQSfyGLKJ7lAfrzP9rwncZyJZ8eRygiUpxinpF7H88zh6EcIIZqnOAqluY/li+abB7mmcLHrFkTqlXTo/fEJwVgEancfPR0T9poltCbVXTlH4wHIJNqXMNifs9b9OFjapJR9DE1cUpKICgCsGYQB05FDkGIBIUierqZVOUjrmYhN7CIgRyhDnU4xP08zzmkMZ27iz6m0spSCgEPwNWqVSM9PZ2YmBgF4QpmrSU9PZ1q1aoFuioi5c/XmK6NJJsoqpPJQm5gKK9Tj3RuYCE3sJBefEIV90MNfFJPV8og4AG4cePG7Nq1i7S009wHJ+WmWrVqNG7cONDVECk/XoHXAl9yKcncygIG8zBP82eeZyCL+Ig+XMkyosn2fSz1dMWPAh6Ao6OjiY+PD3Q1RCRcnOa2ocd4jDkMZSvNqMYJBvBfOpICQG2O0oclvo+pnq6Ug9IuxCEiElx8LJCxx57Hm7gWiDfAai6jOVuYzTB+41wWcBPdWVH4WFoMQypAwHvAIiJl4pViPkkVFjGQWYzgI67GYOnJp8RwgA/oRxQ5RR9LPV2pQOoBi0ho8fVoP3fw/ZjeNGIPg3mTDbRnHP9gE62JcT9RtVDwVU9XAkg9YBEJDT5mMe+39XiDW2jFZnqzlNZsojdLGMEsrmIpkeT6PpZ6uhIEFIBFJLh5Bd4cIviYPsxkJO8xiCyqcC9T6c1SGrObeQwp+lgKvBJEFIBFJDj5er4ucA2LWUIf6pPGn5jGCGbRjo2+j6HbhiSIaQxYRIKDj7Fdm57OcrozgpmcwLVgzChmMJ/B7OZ8nuMB38E3JsY1ppuT47oNads2BV8JOuoBi0hg+ejpHrB1eY3bmc5d/EBLzuYA9/ASnVjHDbxV9LGUYpYQogAsIoFRRIp5K3G0ZhOZVOdSVjObYQxmge8H2CvFLCFMAVhEKpZX4M0iioXcwK+cx/28QBzbGMc/GMgiElnv+xjq6UoY0BiwiFQM50pV6ensJ4aJjCeObdzCXF7jdnIxGOBv/N138M0b29X9uhIGFIBFpHx5BV6AOdxGE3byCBNpy0bepx8pdCSCIh6PqcArYUgBWETKhyPw5qQfZBED2EA7ADqSwjBe5Tta8xF96cdiz+CrFaqkElAAFhH/cgTejPQMXuJuLuIHBrGIl7gHgNZsZjp305rNnvvq9iGpRBSARcQ/vFLNLzCWWLYzmpeoz34WcCPT+JPvfZVilkqo2ABsjJlpjNlnjNnotf1eY8wPxpjvjDGTyq+KIhK0nItnDB3KzvTq5GIASKMBl/AVn9GN1VzGjSws/DAEBV6pxErSA54N9HVuMMZcCQwC2ltr2wCT/V81EQlaXs/e3UA7htpXiWcr73MtAH/nr/yPAXRjpTskOyjwihQfgK21K8D9LK8CdwNPW2tPusvsK4e6iUiwcQTevGUi+/E+CWzgHa5nDFPpwDcAvmc0K/CK5CvtQhwXAlcYYyYCmcCD1tq1vgoaY0YBowCaNm1aytOJSED5WLUqmyiG8SonqM6TTOBu/kU9DvreXwtniBRS2klYUcDZwKXA/wMWGGMKZZkArLUzrLVJ1tqkBg0alPJ0IhIQXj3e9+nHIN4lk6pEk837XMt2YpnAU76Dr3q8IkUqbQDeBbxtXdYAuUB9/1VLRALKEXhz0w/wLoNIYh39eZ8NtOcXmgHQlu8812jO+ztc9++KFKu0AfhdoCeAMeZCoAqw3091EpFA8bqVKI36dOAbruddjlCbmYzgRy4sfP8uuHq7c+bo/l2REirJbUhzgdXARcaYXcaYPwAzgWbuW5PmAcOstUWsISciQc9r1ar1tAegPvtpy0bmcBubacUIZhNNtue+SjOLlEqxk7CstUOKeOs2P9dFRCqaY3KVBd7md/yNJ9hOLNuIoz7pJBd1qWtilUiZaCUskcrIa3LVYvqS5H7YfS4RzGIE9QrdfeimHq+IX+h5wCKViY/biX7gIvqxmDi2Mpth3MbrRJJbeF/1eEX8SgFYpDLwCrxrSWIF3XiA52jJDyymLz35lCpkFd5XgVekXCgFLRLOvGY1f0tbruMdOrOWSTzEEWoB0JePCgdfpZpFypUCsEg48gq8ezmPobxGAutZxpU8wV/ZQnNqc7Twvgq8IhVCKWiRcOJjjBcglwg+oB8PMYmHmFT0qlVKNYtUGAVgkXDgFXhPUoVp/IlVdOUtfs/57GEnTajBicL7KvCKBIRS0CKhzCvVnIvhDYbQku95kGc5QXWOcRZA4eCrVLNIQCkAi4Sq5GQYNSq/1/szzejMGm7lDc7mIEu4isX0oxbHPPdT4BUJCkpBi4Si5GQYNgxycjhJFapyivP4lSqc4nVuZQhzCz+PV6lmkaCiHrBIKHGknPfn1GU000gklVNEU5MMvqArt/KGZ/BVj1ckKCkAi4QCR+DNSj/MFMbQgp/4N3dyFUs5SdXC+yjwigQ1paBFgpnX7OZdnE9vlvA9rejNxzzP/bRhk+c+SjWLhAQFYJFg5BV4j1ODmmTQiD0ksJ5JPER//odx7hMZCa++qsArEiKUghYJJl63FR2iDvfzHM34hXTqEYFlHkMY4B18a9RQ8BUJMQrAIsHCcVtRDhH8m1G04CemMJZBvFf0fjExMGOGgq9IiFEKWiQYOG4rOkZNerCcFJLoxme8wH10ILXwPhrrFQlp6gGLBJIj5XwyJxKAszjOZazmDYawnB6Fg69mN4uEBQVgkUBwBN7c9APM4A7i2MZmWgLwT8YwhHme47wKvCJhRQFYpCJ5TbJaSxKX8iV3MoML+ZEIcgvvo8ArEpYUgEUqimOSlQXuZSqX8BU7aUIyt7CcHlzEjwXlIyMVeEXCWLEB2Bgz0xizzxiz0cd7DxpjrDGmfvlUTyRMuCdZ2YwMAAxQi6Pcxwv8wEXcwlzdViRSyZSkBzwb6Ou90RjTBOgN7PBznUTChyPlvCnnQrrzGUvpBcBTTOA5HqA2Rz330W1FIpVCsQHYWrsCOODjreeBh8D7kSsi4gy8menH+CtPkEgq39GGI9T2vY/GekUqlVLdB2yMGQjsttauN8YUV3YUMAqgadOmpTmdSOjwWkJyGT0YxQy20IKhvMazPEAD9nvuo/t5RSqlMw7AxpgawASgT0nKW2tnADMAkpKS1FuW8JU3yco9zgvwvfu2oqX0ohefepbX2s0ilVppZkFfAMQD640x24DGwNfGmPP8WTGRkOKYZDWL4czlZgDu5N9soH3h4KtJViKV3hkHYGvtt9bac6y1cdbaOGAXcLG19le/104k2DnGen/IuYArWcZIZjHPHYAjsFQn03MfTbISEUp2G9JcYDVwkTFmlzHmD+VfLZEg5wi8WemHmch42rOB9SQwgzt4h+sL76NJViLiUOwYsLV2SDHvx/mtNiKhwGusdyVX8AgTGcx8pjCW8/jNs7wmWYmID3oaksiZcI/1nsiJZhW9uIpP6Mky1pJEEimeZTXJSkROQ0tRihQnORni4sAYGDqUlTmXkUgq/fiAPTQEKBx8NclKRIqhACxSFOeDE7Zv5yhnMdr+k26s5BRV+IB+NGJv4f00yUpESkApaBFfvMZ5M6lKIqlsJZ6xvMBEJlCTDM99NNYrImdAAVjEm3ucl5wcTlCN6mRSjZPcz/NczNd0YbVneY31ikgpKAUtkseZcs7J4X360ZwtLOEqAP7Ei4WDr8Z6RaSUFIBFnIE3PZ3D1GYkr9Cf94khnQakeZbPW/88NlZjvSJSakpBS+XmNda7lF6MYBZ7aMR4JvI3nqAqpwrKa5xXRPxEAVgqL8dYb54tNOcsjrGay+jM2oKyGucVET9TCloqH6+x3mX0YCG/B1wPT/iGDp7BV+O8IlIOFIClcslLOaenc5wa3MtUerKMSTyEBQxQjZMF5XVPr4iUEwVgqTzyUs4ZGXxOVxJYz4uM5j6eZzk9MM6yenCCiJQzjQFL+EtOhrFjIT0dgM20pBsriGcry+lBN1YWlNVYr4hUEPWAJXx53V60nxgAWvE9rzKM9SR4Bl+N9YpIBVIAlvDkGOvNJpKJjCeW7aSSAMBQXucsjheU11iviFQwpaAl/DhuL/qZZgxlDqvpwk3Moyk7PMsq5SwiAaIesIQPr9uLXmEkCaxnE615gyHMYwj1OFhQXilnEQkg9YAlPHitaAWwlXgu4StmM5wm7PIsrxWtRCTAFIAl9DlSzosYQC2OciXLeYzHiCCXCGxBWQVeEQkSSkFL6HKknI/mVOeP/IdBLOI5/gxAFDkFwTcyUvf1ikhQUQCW0ON1e9EqupBIKrMYwTie4i33spL5NNYrIkGo2ABsjJlpjNlnjNno2PZ/xpjvjTEbjDHvGGPqlmstRfI4bi8C+JJL6MYKLIbP6M5TTKAKWQXldXuRiASpkvSAZwN9vbYtAdpaa9sDPwLj/FwvkcIcS0mepAoAl/AVz/IA60ngclYVlFXKWUSCXLEB2Fq7Ajjgte1ja222+8cvgcblUDcRF0fK2ebkMJthNGcL22mKAe5jCrU4VlBeKWcRCQH+GAMeCSwu6k1jzChjzDpjzLq0tDQ/nE4qFUfK+SB1uZl5jGA2F/AzkeQULq+Us4iEiDIFYGPMBCAbSC6qjLV2hrU2yVqb1KBBg7KcTiobR8p5BVeQwHre5nc8xTg+oReN2V1QVk8vEpEQU+r7gI0xw4D+QC9rrS2uvEiJeT29CGA2w6nKSb6gC51YV1BWS0mKSIgqVQA2xvQF/gJ0t9ZmFFdepMQcK1r9TDNOUYVWfM9UxgB4PkChRg2lm0UkZJXkNqS5wGrgImPMLmPMH4BpQC1giTEm1RgzvZzrKZWBO+VsMzJ4ldtJJJVRzABcgVdPLxKRcFJsD9haO8TH5lfKoS5SWTlSzoeow11MZz43043PmMNQz7JKOYtImNBa0BJYjpTzFi6gF5+wm/OZyHj+wjNEkltQVilnEQkjCsASOI6HKAA0ZQeX8BUPMpnOrPUsq4coiEiY0VrQUvEcC2vszGnIEN7gIHWpQhYLuMkz+GpFKxEJUwrAUrEcC2u8x0ASWM//6M96EgqX1YpWIhLGFICl4rhTziczshnDFK7jPeLZytdcTA8+8yyrWc4iEuYUgKX8OR8fmJPDAzzLPxnDfTzPF3ShBVsKyirlLCKVhCZhSflyzHI+QTWqk8kEJtKXD+nP+55lNctZRCoR9YCl/LhTzscyDLfzKtfyPjlE0JBfCwdfpZxFpJJRABb/c6Scv8lpx8V8TTK30o0Vhcsq5SwilZQCsPiP87m96en8kz9xKV9ynJp8Qi8e4/HCC2tolrOIVFIKwOIfjtuLAI5xFs/yAL1ZwnoSNMtZRMSLJmFJ2TlWtFpLEu3ZQC2O8QVdaMhejLOs1nIWEQHUA5aycKScc3IsTzKBS/mSSTwEQCPv4KuUs4hIPvWApXQctxftoSG38TrL6MktJDOWKYXLay1nEREPCsBy5hwp5+V050beJIMazGI4w3jVs9erwCsi4pMCsJSc47m9eRqQRnO2MIsRtOSHgrIa6xUROS2NAUvJOGY5/0I8TzEOgDZs4gu6eAZfjfWKiBRLAViKl5dyzshgHjeRSCr/x/9jJ40BCqecdXuRiEixFIClaI5ZzsdzqvJH/sMQ5tGOb0klkSbsKiirFa1ERM6IxoDFN8csZwtczUd8QRfGM5HHeZQocgrK6iEKIiJnTAFYPCUnw4QJsH071r3JAA/zNNXI5Co+8SyvWc4iIqVSbAraGDPTGLPPGLPRsa2eMWaJMeYn9/ezy7eaUu6cz+zdvp1D1GEwC5jCWAD6875n8FXKWUSkTEoyBjwb6Ou17WHgE2ttC+AT988SqrzWcf6SS0gklXe5Dus5xcpFs5xFRMqs2ABsrV0BHPDaPAh41f36VeA6/1ZLKoxjhnMuhmd4iMv5HINlJVdwPy94ltcsZxERvyjtLOhzrbV7AdzfzymqoDFmlDFmnTFmXVpaWilPJ37nTDnnuCZUfUMHxvEPrucdvqEDl/JVQfnYWKWcRUT8qNwnYVlrZwAzAJKSkmwxxaUiOGY4A2wljni20ZGvWUNnOpJSkHjWDGcRkXJR2h7wb8aYhgDu7/v8VyUpV46UcxZRjOMpWvATK7kcgCRn8FW6WUSk3JQ2AC8ChrlfDwPe8091pNx4pZy3EUt3PuNpxjGSmXQkpaCsZjiLiJS7YlPQxpi5QA+gvjFmF/Ao8DSwwBjzB2AHcGN5VlLKyCvl/A7XMZKZ5BLBPG7iJhYUlFXKWUSkQhQbgK21Q4p4q5ef6yLlwfHowDzbiaUFPzGPm2nG1oKyWlRDRKTCaC3ocOWVct5MSz6iDwBjmcIquhYEX6WcRUQqnAJwOHIsrGGBVxhJR1IYw1RyiMAA0WS7ympRDRGRgFAADjeOWc6Hqc0tvMEfeYUufMFyehBJbkFZzXIWEQkYPYwhXCQnw9ix+ctJHuBsOrGW7cTyFON4iEkFwTcyUr1eEZEAUwAOB16znAHqcZAbeZOBLKILqwvKapaziEhQUAo61DlSzvtowO94i420AeBpxnkGX6WcRUSChgJwqPKa5fwJPUlgPR/Qj+/cATifZjmLiAQdBeBQ5JjlnEUUE3iS3izhbA6yhs6FF9bQeK+ISNBRAA41jpQzwEvcw1NMYCQzWUsn2vNtQVmlnEVEgpYmYYUKr1nOh6hDXQ5zF9O5gJ/pz/sFZTXLWUQk6KkHHAocKecTVONuXiKRVA5Rh6qc8gy+SjmLiIQEBeBg50g5b6Yll/AV07mbwSygBhmeZZVyFhEJGUpBBytHytkCMxnJvfyTszjGYvrSl48KyirlLCISchSAg5HXwhoWwxvcQhe+YA5DacivBWW1sIaISEhSAA42jscHrqIL8WylEXt5m99xFscKr+WsxweKiIQkjQEHC8fCGtk58BiP0o0VPMKTANThiOdazlpYQ0QkpKkHHAwcKedtxHIryXxBV27nVV7gPs+ySjmLiIQFBeBAc6ScP6cr17pvKUrmFm5hrmdZpZxFRMKGUtCB4rWWM0BbNtKXD0kl0TP4KuUsIhJ2FIADwbGwxho6cSMLOEkV6nKY+dxMPNsKymphDRGRsKQAXNHcKeecjEz+wcN0ZRVr6MwOmhYuq4U1RETCVpkCsDHmfmPMd8aYjcaYucaYav6qWNhxpJx35ZzHVSxlPP/gd7zNehJowZaCsko5i4iEvVIHYGPM+cAYIMla2xaIBG72V8XChnOs1/0ghdt4nbV0YiYjmMfN1OVwQXmlnEVEKoWyzoKOAqobY7KAGsCeslcpjDhuLzpGTSyGWhzjX9xNFNmevV7QLGcRkUqk1D1ga+1uYDKwA9gLHLbWfuxdzhgzyhizzhizLi0trfQ1DTWOhyh8RWc68A1jmApAK75XyllEpJIrSwr6bGAQEA80AmoaY27zLmetnWGtTbLWJjVo0KD0NQ0VPla06soqTlGFEcwqXF4pZxGRSqksk7CuArZaa9OstVnA20AX/1QrRDluL9pKHJfzOY/zGLfwBhtoTzdWepbXLGcRkUqrLGPAO4BLjTE1gBNAL2CdX2oVihwrWgFEksM+zmEBN3IjCz3LaqxXRKTSK8sY8FfAQuBr4Fv3sWb4qV6hw5Fy/i0nhr/zCBZoyk5+5ELP4KuxXhERcSvTLGhr7aPAo36qS+hxzHL+L/35A69whNoMZBEJbCCKnIKyeoiCiIg4aCWs0nKnnI9lGEbxbwbyX85nNyl0JIENnmU11isiIl70NKQzlZwMY8fmL6oxkEUspwd/4Wke51GqcqqgbGSkZjiLiIhPCsAl5Qi8p4gGoqlCFo/xGBZDd1Z4llfKWURETkMp6JJw3F60nvZ0Yi1/568AdGNl4eCrlLOIiBRDAbg47rHe7IyTTGQ8nVjLb5xLZ9YULqtZziIiUkIKwEVx3F70Y04zurKKR5jI73ib72jDAP7nWV4rWomIyBnQGLA3r0lWACeozk6aMJ/BDObNwvtoYQ0RETlDCsBOjvt6fyGed7ieB3iOBDawlXjPGc6gwCsiIqWmFHQe91ivzchgOnfSng08wd/YQ0OAwrcXaaxXRETKQAHYMda7M6chffmQu5lOF75gI21pxF7P8hrrFRERP6jcKWhHyjmLKK5gJfupz7+4izv5N8a7vFLOIiLiJ5UzACcnw4QJsH07u2lEQ04QTTYvcQ8X8QMX8ItneQVeERHxs8qVgnakm3O37+BFd8B9mT8C0I/FnsFXY70iIlJOKkcAdgRe0tP5ieb0YDl/4kUu53Ou5qPC+2isV0REylF4B2CvwAswkxG0ZwPf0o5ZDGcx1xDLDs/9tJSkiIiUs/AdA3ZMsHJqwk6u5iNe4p7CM5xjY2HiRAVeEREpd+EZgN339JKTQxZRTOIhsoniUZ6gN0vpzVLP8npykYiIVLDwSkE7U845OXxFZzqSwiNM5CdaYH3to3SziIgEQHgEYK+x3qOcxRimcBmrOUA93uE6Xmeo5329MTGa4SwiIgET+gHY8azePDtoygxGMZoX2URrruO9gvIKvCIiEgRCewzYMda7h4a8xe+5l2m0YRO/0MxzklVkpG4rEhGRoFGmHrAxpq4xZqEx5ntjzGZjzGX+qthpORfUyMllOnfSis08xCS2EQvgGXx1T6+IiASZsqagpwAfWmtbAgnA5rJXqRiOlPN3tOYKVnI300liHRtoTxzbPctrkpWIiAShUqegjTG1gW7AcABr7SnwfmBuOZgwATIyyKQqPfmUbKKYzTBu57XCk6y0frOIiASpsowBNwPSgFnGmAQgBRhrrT3uLGSMGQWMAmjatGkZTue2w7VqVTVOMo+bactGGrC/4H2N9YqISAgoSwo6CrgY+Je1tgNwHHjYu5C1doa1Nslam9SgQYMynM7NEcSvZLln8NVYr4iIhIiyBOBdwC5r7VfunxfiCsjla+JEV6D1prFeEREJIaUOwNbaX4GdxpiL3Jt6AZv8UqvTufVWV6CNjQVjXN91X6+IiISYst4HfC+QbIypAvwCjCh7lUrg1lsVbEVEJKSVKQBba1OBJP9URUREpPII/aUoRUREQpACsIiISAAoAIuIiASAArCIiEgAGGt9Pqa+fE5mTBp4L9ZcJvXBuRJHSFNbgk+4tAPUlmAVLm0Jl3aA/9sSa631uQpVhQZgfzPGrLPWhsUsbLUl+IRLO0BtCVbh0pZwaQdUbFuUghYREQkABWAREZEACPUAPCPQFfAjtSX4hEs7QG0JVuHSlnBpB1RgW0J6DFhERCRUhXoPWEREJCQpAIuIiARAwAKwMaaJMWaZMWazMeY7Y8xY9/Z6xpglxpif3N/PduwzzhizxRjzgzHmasf2jsaYb93vTTXGmCLO6XP/YGmLMaaGMeZ9Y8z37uM8XcT54owxJ4wxqe6v6cHWFvf25e5tefU8p4hzBvvnUsvRhlRjzH5jzAs+zlcun8uZtsMYE+Muf8wYM83rWCF1rRTVllC8Vor5XELqWjnN5xJq10pvY0yK+5pIMcb0dByr/K8Va21AvoCGwMXu17WAH4HWwCTgYff2h4Fn3K9bA+uBqkA88DMQ6X5vDXAZYIDFwDU+zlfk/sHSFqAGcKW7TBVgZRFtiQM2hsDnshxIKuZ8Qf+5+DhuCtCtoj6XUrSjJnA5cBcwzetYoXat+GxLiF4rp/tcQu1aKbItIXatdAAauV+3BXZX5LUSsB6wtXavtfZr9+ujwGbgfGAQ8Kq72KvAde7Xg4B51tqT1tqtwBagszGmIVDbWrvaun4jrzn2cfK5fzC1xVqbYa1d5j7OKeBroLE/6lhS/mrLGZwy6D8X5zGNMS2Ac3D9h18hzrQd1trj1trPgUzncULxWimqLaF4rRTVljMQ9J+LU4hcK99Ya/e4t38HVDPGVK2oayUoxoCNMXG4/hL5CjjXWrsXXL9MXB8guH6JOx277XJvO9/92nu7t6L296sytsV5nLrAAOCTIk4Vb4z5xhjzmTHmCr81wLMOcZS9LbPcKaa/FpHCCanPBRgCzHdflL6U6+dSwnYUJRSvlZIcpy6hca0UJ5SulZIItWvl98A31tqTVNC1ElXSguXFGHMW8BZwn7X2SBFpdnClAbzZ02wv6f5+44e25B0nCpgLTLXW/uKj7F6gqbU23RjTEXjXGNPGWnukbC1wVNA/bbnVWrvbGFPLfayhuP6SLOn+fuGvz8XtZlzt8KVcP5czaEeRh/CxLdivleKOE0rXyumE2rVSEiFzrRhj2gDPAH3yNvko5vdrJaA9YGNMNK5fUrK19m335t/c3f+8lNk+9/ZdQBPH7o2BPe7tjX1s91bU/n7hp7bkmQH8ZK19wde53OmOdPfrFFzjDhf6qSl+a4u1drf7+1HgDXynZkLmczHGJABR7t95IeX5uZxhO4oSitdKcULpWilSCF4rxR0rZK4VY0xj4B3gdmvtz+7NFXKtBHIWtAFeATZba59zvLUIGOZ+PQx4z7H9Znd+Ph5oAaxxpxOOGmMudR/zdsc+Tj73D6a2uI/1JFAHuO8052tgjIl0v27m3t/XX/8Ba4sxJsoYU999zGigP7DRxylD4nNxG4Krt1XU+crlcylFO3wK0WvldMcKtWulqOOE4rVSnJC4Voxr+OJ9YJy1dlVe4Qq7Vmw5zA4syReuGXQW2ACkur/6ATG4xnJ+cn+v59hnAq6/lH7AMSMNSML1D/ZnYBrkr/A1EHiiuP2DpS24/nqyuCYO5B3nj95twTVW8R2u2XdfAwOCsC01cc2A3OCu6xQKZkeH1OfieO8XoKXXtnL/XErZjm3AAeAYrr/SW4fwtVKoLYTuteKrLaF6rfj8NxZK1wrwCHDcUTYVOKeirhUtRSkiIhIAQTELWkREpLJRABYREQkABWAREZEAUAAWEREJAAVgERGRAFAAFhERCQAFYBERkQD4/8uD+3ihpc4/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c1 = -5.0 * np.log(0.14 * np.exp(400) / 6.0)\n",
    "y_nltrue = 0.2 / (np.exp(-0.2 * (t + c1)) + 0.01)\n",
    "fig2 = plt.figure(figsize=(8,4))\n",
    "plt.plot(t,yn2,'ro:',label='Runge Kutta')\n",
    "plt.plot(t,y_nltrue,'b--',label='Exact')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "25wHP_N_X84Z"
   },
   "source": [
    "Once more the plots are indistinguishable, again suggesting that our numerical scheme has found an approximate solution which closely matches the true solution. Once more, a table should again highlight exactly how accurate the approximation is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XsDAJtseqTR2"
   },
   "source": [
    "### Results Table\n",
    "The table below shows the time and the Runge Kutta numerical approximation, $w$,  for the non-linear population equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 363
    },
    "executionInfo": {
     "elapsed": 264,
     "status": "ok",
     "timestamp": 1646995609489,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "2nvEBzeWqTR2",
    "outputId": "beb53a8b-e079-4762-d7ad-9bf7c8114e1c"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time t_n</th>\n",
       "      <th>Runge Kutta (y_n)</th>\n",
       "      <th>Exact (y_true)</th>\n",
       "      <th>Exact Error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2000.0</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2000.1</td>\n",
       "      <td>6.084335</td>\n",
       "      <td>6.084335</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2000.2</td>\n",
       "      <td>6.169332</td>\n",
       "      <td>6.169332</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2000.3</td>\n",
       "      <td>6.254983</td>\n",
       "      <td>6.254983</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2000.4</td>\n",
       "      <td>6.341278</td>\n",
       "      <td>6.341278</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2000.5</td>\n",
       "      <td>6.428207</td>\n",
       "      <td>6.428207</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2000.6</td>\n",
       "      <td>6.515759</td>\n",
       "      <td>6.515759</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2000.7</td>\n",
       "      <td>6.603924</td>\n",
       "      <td>6.603924</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2000.8</td>\n",
       "      <td>6.692689</td>\n",
       "      <td>6.692689</td>\n",
       "      <td>1.000000e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2000.9</td>\n",
       "      <td>6.782043</td>\n",
       "      <td>6.782043</td>\n",
       "      <td>2.000000e-10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   time t_n  Runge Kutta (y_n)  Exact (y_true)   Exact Error\n",
       "0    2000.0           6.000000        6.000000  0.000000e+00\n",
       "1    2000.1           6.084335        6.084335  0.000000e+00\n",
       "2    2000.2           6.169332        6.169332  0.000000e+00\n",
       "3    2000.3           6.254983        6.254983  1.000000e-10\n",
       "4    2000.4           6.341278        6.341278  1.000000e-10\n",
       "5    2000.5           6.428207        6.428207  1.000000e-10\n",
       "6    2000.6           6.515759        6.515759  1.000000e-10\n",
       "7    2000.7           6.603924        6.603924  1.000000e-10\n",
       "8    2000.8           6.692689        6.692689  1.000000e-10\n",
       "9    2000.9           6.782043        6.782043  2.000000e-10"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'time t_n': t[0:10],'Runge Kutta (y_n)':yn2[0:10],'Exact (y_true)':y_nltrue[0:10],\n",
    "     'Exact Error':np.abs(np.round(y_nltrue[0:10]-yn2[0:10],10))}\n",
    "df = pd.DataFrame(data=d)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7-YP7NTlYnW-"
   },
   "source": [
    "The table illustrates that for this equation, a very small (but nevertheless non-zero) error creeps in after a few timesteps, and appears to grow every few steps.\n",
    "\n",
    "The first few steps in the table are likely to be zero due to rounding: this worksheet displays numbers to a preset number of decimal places."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "alLm7_PfqTR4"
   },
   "source": [
    "## 3. Non-Linear Population Equation with an oscillation \n",
    "In our final differential equation, lets yet again use a more complex population model and evaluate how well RK4 does. The population model is essentially the previous model with an additional sine term (representing shorter fluctuations in population than the long term behaviour):\n",
    "\n",
    "$$ \n",
    "y^{'}=0.2y-0.01y^2+\\sin(2\\pi t), \\ \\ (2000 \\leq t \\leq 2020), \n",
    "$$\n",
    "\n",
    "with the initial condition,\n",
    "\n",
    "$$\n",
    "y(2000)=6.\n",
    "$$\n",
    "\n",
    "### Implementing RK4\n",
    "As with the previous cases, lets express the right hand side of the ODE as a function\n",
    "\n",
    "$$\n",
    "f_3(x,y)=0.2y-0.01y^2+\\sin(2\\pi x),\n",
    "$$\n",
    "\n",
    "Once again, we could mathematically calculate expressions for $k_1\\rightarrow k_4$. \n",
    "This time though, you can see that $f_3(x,y)$ now contains terms that depend on $y$ **and** $t$!\n",
    "\n",
    "Now that we're more familiar with how to implement RK4 numerically, we'll simply move on to define our function representing $f_3(x,y)$ and let the code evaluate the coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "id": "KcPh153LqTR4"
   },
   "outputs": [],
   "source": [
    "def nlo_fun(x, y):\n",
    "    ftw = 0.2 * y - 0.01 * y * y + np.sin(2 * np.math.pi * x)\n",
    "    return ftw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "id": "XqdqPm6ZqTR6"
   },
   "outputs": [],
   "source": [
    "yn3 = np.zeros(N + 1)\n",
    "yn3[0] = 6.0\n",
    "## 4th Order Runge Kutta\n",
    "for n in range (0,N):\n",
    "    k1 = nlo_fun(t[n],yn3[n])\n",
    "    k2 = nlo_fun(t[n] + 0.5 * h, yn3[n] + 0.5 * h * k1)\n",
    "    k3 = nlo_fun(t[n] + 0.5 * h, yn3[n] + 0.5 * h * k2)\n",
    "    k4 = nlo_fun(t[n] + h, yn3[n] + h * k3)\n",
    "    yn3[n + 1] = yn3[n] + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EqJ4WI_VqTR8"
   },
   "source": [
    "### Results Graph\n",
    "This final ODE is quite a lot more complicated than the two previous ones (note for example, that it's not separable). We could employ more advanced mathematical techniques to derive an analytical solution, but that isn't the focus of our work here, so for now we'll just focus on the RK4 solution that we obtain from our numerical technique \n",
    "\n",
    "The plot below shows the 4th order Runge Kutta numerical approximation, $y$ (circles) for the non-linear population equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1646995590717,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "iHNXCSpHqTR8",
    "outputId": "b7ab46e9-602b-4e10-cf8d-33706e2a9ddf"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD4CAYAAAA0JjXXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbCUlEQVR4nO3dfYwc9X3H8c/XZyfmCATYOwLB9p1TBZpghadTGiqaJ0pLqwCRWlVB54QWWqtISQ1S2oLcJiqqK0LogxOSRk4hEHyFoDQPJFGQEUpKUxHSszHEQBISsOGAxGcb8mQebPPtHzNjr9ezO7O7MzsP+35Jp92b3Z35/W4Zvv795jffr7m7AADAYC0ougEAAAwjAjAAAAUgAAMAUAACMAAABSAAAwBQgIWDPNjY2JhPTk4O8pAAABRm06ZNO919PO61gQbgyclJzc7ODvKQAAAUxsy2t3uNKWgAAApAAAYAoAAEYAAACjDQa8Bx9u7dq7m5Ob344otFNyVzixcv1pIlS7Ro0aKimwIAKJnCA/Dc3JyOOuooTU5OysyKbk5m3F27du3S3Nycli9fXnRzAAAlU/gU9IsvvqhGo1Gr4CtJZqZGo1HLkT0AVMXMjDQ5KZlJCxcGj5OTwfaiFR6AJdUu+Ebq2i8AKLuZGWlsTFq5Utoe3gi0f3/wuH17sH1srNhAXPgUNAAAWZmZkVavlnbtSn7vrl3SqlXB8+npfNsVpxQj4KKNjIzo9NNP14oVK3TBBRfo+eeflyRt27ZNK1asOPC+z372szrzzDP13HPPHdh2/fXXy8y0c+fOQTcbAGqnnynjmZkgoKYJvpE9e6RLLilmJFy5ABx9OQsWZDePf8QRR2jLli3aunWrjjvuOH3qU5867D233nqrPvnJT2rjxo069thjJUlPPfWU7r77bi1btqz/RgDAEOt3ynhmJgike/Z0f+z9+4P9j4wM9hpxpQJw9K+b7dsl9+Bx1aps/1Bnn322nn766UO23XHHHbr22mu1ceNGjY2NHdh+5ZVX6rrrruNaLwD0qDnwJo1cd+2KD8RRbIgCdq9eeSV4zCO2xCldAH7nO6Wbbw6e790b/L5hQ/D71Vcf/q+bPXukK64Inu/cGbz/a18Lfv/pT7s79v79+3XPPffowgsvPLBt+/bt+uAHP6iNGzfqhBNOOLD9zjvv1EknnaTTTjutu4MAQE30u8K4lylj6WAgjkasK1f2NvLtZM8eac2abPfZqnQBuJO5ufjt3X55rV544QWdfvrpajQa2r17t84777wDr42Pj2vZsmW64447Dmzbs2eP1q5dq2uuuaa/AwNABWWxwrifKeNINGJN0mgEAzn34HFkJN3nnnyy97al4u4D+znrrLO81SOPPHLYtnYmJtyDP+GhPxMTqXcR68gjj3R39+eff97POeccX7dunbu7P/HEE37qqaf6jh07/OSTT/YNGza4u/tDDz3k4+PjPjEx4RMTEz4yMuJLly71Z599tq/+AUCZbdjg3mjE/3+43U+jEXyun330+jMycuixm9swOpr8+X5ji7u7pFlvExMrFYDj/mijo/F/4G5EAdjdffPmzb506VJ/+eWXDwRgd/fHH3/cly1b5nfddddhn5+YmPD5+fnYfROAAdRB2qDV7mfBguDRbDDBNyk2JP1DIIvY4t45AFdqCnp6Wlq/XpqYCOb9JyaC37O8f+uMM87Qaaedpttvv/2Q7cuXL9edd96pSy+9VPfff392BwSAAen1mm2W08Xuye9tnTJuNLo71shIcmyYng7WDW3YEMSS6HNSPrEljnmav0ZGpqamfHZ29pBtjz76qN70pjcNrA2DVvf+ASi/NMkpGg1p3brDg060UCrrRU5xRkakW26JD3xpE2yMjg4meKZlZpvcfSrutUqNgAEA3Um70jjuFp8sRr5pjY62D75S5xHroEeuWSEAA0BN9RJAo0Ac3d6TdG9tv9PF0T7SBs7paWnbtuB4+/Yd+rhtW3WCr1SSADzIafBBqmu/AJRb821C/SanaGdkJAi4O3ceDHrNo9ROgTjKXTQxcfg+hknhAXjx4sXatWtX7YKVe1APePHixUU3BUBFxS2aSlo81Wtyi270M108MSHdems1R6xZK3wR1t69ezU3N1fLurmLFy/WkiVLtGjRoqKbAqBCuqno07x4KppyzmvUK3VeKIXDdVqEVXgABgAEugm8vWoOoN0er2wrjKuAVdAAUHJFTB2nvWYrdbdQCukkBmAzu8nMdpjZ1qZtp5vZd81si5nNmtlb820mANTXIG736RRAk67ZDvNCqTylGQHfLOn8lm3XSfoHdz9d0kfC3wFgaPWTZSqLUnrtxK1WbifuFp9hXyiVp8QA7O73StrdulnS0eHz10p6JuN2AUAlJFUGev/72wflbka+C8L/W3dTfjxptTKK1es14CskfdzMnpJ0vaSr273RzFaF09Sz8/PzPR4OAMonzXXbaJ1ra7m+bhNd7N8f7OuVV7hmWxe9BuDLJV3p7kslXSnpxnZvdPf17j7l7lPj4+M9Hg4AyiXv67adpo65ZlsPqW5DMrNJSV939xXh7z+XdIy7u5mZpJ+7+9Gd9iFxGxKAesi7QAG3+9RHHrchPSPpHeHzd0t6rMf9AEClDGLkS/AdDmluQ7pN0n2STjGzOTO7TNJfSPpnM3tQ0j9JWpVvMwEgW72sWk67YrmXBVMSi6aGTZpV0Be7+4nuvsjdl7j7je7+HXc/y91Pc/ffcvdNg2gsAPQradXyypXBKLSXFcvRddtowdSttx68RpsUjFk0NXxIRQlgaOR57TbNdduZGWnNmiDQj4wEgXpiQlq7lsBbV52uAS8cdGMAoAh5FipIe912eppAi4PIBQ2g9vLMNsV1W/SKETCAWivDyBeIwwgYQG0x8kWZEYABVEp0+9CCBZ1vG+pm1bL7oVmlWLGMQSAAA6iE1tuH3NvfNpRm5BtXGzeqBBTlWybNI/LEbUgASm1mRlq9OttC9SMjTB9jMLgNCUAl5XHfLnmWURYEYACllMfqZVYto0y4BgygdPJYvcyqZZQNI2AApcLIF8OCETCA0kg78m00pMsvT3fbECNflBUBGEApdHPf7s6d0qc/ne62IUa+KCsCMIBcxNXbjau723x/bzf37bZqvo93377gcds2gi/Ki/uAAWSK+3aBg7gPGMBAcN8ukB4BGEAmWL0MdIdrwAD6xn27QPcIwAD6tmZNttPOVBvCMCAAA+hZtNJ5+/bO71sQ/p8mTZk/qg1hWCQGYDO7ycx2mNnWlu0fMrMfmtnDZnZdfk0EUEbRtHOn4Bvdt7t/P2X+gFZpFmHdLOkGSZ+PNpjZuyRdJOkt7v6SmR2fT/MAlFGaBVftVi9PTxNkASnFCNjd75W0u2Xz5ZKudfeXwvfsyKFtAAYsLnlGL0kzyEAFJEuViMPMJiV93d1XhL9vkfRVSedLelHSh939/9p8dpWkVZK0bNmys7YnXSwCMHBZJs+YmAgyUAHonIij10VYCyUdK+ltkv5a0h1m8csr3H29u0+5+9T4+HiPhwOQh+YRbRbBd3RUWru2//0Aw6DXADwn6Use+J6kVySNZdcsAHmLFlFllTKSpBlAd3oNwF+R9G5JMrOTJb1K0s6M2gQgZ2kqD3WDpBlA99LchnSbpPsknWJmc2Z2maSbJL0hvDXpdkmX+CCrOgDoWdZZq0iaAfQm8TYkd7+4zUsrM24LgAHIKmtVoyGtW0fgBXpFJixgiMzMJGetirJRuZM0A8gT1ZCAIRFNPbcTV3OXpBlAfhgBA0MgadEVi6iAwSMAAzWWNnMVi6iAwSMAAxUUlzKy9XHBgnQJNiYmCL5AEbgGDFRMdC03mk6ORratj2luDCRzFVAcAjBQIWmqEKVF5iqgWExBAxWRZQINFl0BxSMAAxWQZepIMlcB5cAUNFByWY18yVwFlAsjYKDk0qaOjDJVtT6SuQooJ0bAQIklpY4cHWU6GagqRsBASaVJHUnwBaqLAAyUEKkjgfojAAMlk2bRFSNfoPoIwMAApUkhuXJl50VXpI4E6oFFWMAAzMxIq1cfmpe5XQrJTkgdCdQHARjIWWvu5l6x6AqoF6aggRxllcGKRVdA/RCAgZxklcGKkS9QTwRgIAeMfAEkSQzAZnaTme0ws60xr33YzNzMxvJpHlA9aUe+C8Kzr1MKSUa+QH2lGQHfLOn81o1mtlTSeZKezLhNQKUl5W4eGQlyM+/fL7lL+/bFP27bRvAF6iwxALv7vZJ2x7z0r5L+RpJn3SigqtLkbmZKGYDU4zVgM7tQ0tPu/mCK964ys1kzm52fn+/lcEAlkLsZQDe6DsBmNippjaSPpHm/u6939yl3nxofH+/2cEAlkLsZQLd6GQH/hqTlkh40s22SlkjabGYnZNkwoCrI3QygF11nwnL370s6Pvo9DMJT7r4zw3YBlRCNfDsFX3I3A4iT5jak2yTdJ+kUM5szs8vybxZQvLjCCWNjwY9ZcBvRypWdgy+5mwG0kzgCdveLE16fzKw1QAl0KpzQvM0T1v+z6ApAJ2TCAkIzM8HoduXKQwNtL1h0BSAJ1ZAAZVexSGLkCyAdRsCAkrNXpcXIF0BaBGAMvaTsVWk1Gox8AaRHAMZQS8peFYkKJzQawY90aNGEDRuknTsJvgDS4xowhlrS1HOjIa1bR2AFkD0CMIZW0tTzhg0EXgD5YQoaQylp6pnsVQDyxggYQycpfSTZqwAMAiNgDBUKJwAoC0bAGBoUTgBQJoyAURtxxROiRwonACgbRsCovE7FE6JHCicAKBsCMCotixzOo6MEXwCDRwBGZaW5ppuEkS+AonANGJWUZjVzEgonACgSARiV1G/1IgonACgaARiVk6Z6UVQ8ISqYQOEEAGXDNWBUSlIKyZERppUBVAMjYFRKp6lnrukCqBICMCojaeqZa7oAqiQxAJvZTWa2w8y2Nm37uJn9wMweMrMvm9kxubYSQ4/qRQDqJs0I+GZJ57dsu1vSCnd/i6QfSbo643YBB0T3+3aaeiaFJICqSQzA7n6vpN0t2za6+77w1+9KWpJD2wCqFwGorSyuAV8q6ZvtXjSzVWY2a2az8/PzGRwOwyJp5Csx9QyguvoKwGa2RtI+STPt3uPu6919yt2nxsfH+zkcaiquihHViwDUXc/3AZvZJZLeI+lc96RaM8DhOlUxonoRgLrrKQCb2fmS/lbSO9y9j4SAGEZxgbcbVC8CUAdpbkO6TdJ9kk4xszkzu0zSDZKOknS3mW0xs8/k3E7URLSoqtfgy8gXQF0kjoDd/eKYzTfm0BbUXL/lAxn5AqgTMmFhIPotH0j1IgB1QwBG7tLcThRprWJE9SIAdUU1JOQq7ci30ZDWrSPIAhgeBGDkqlP1IonygQCGF1PQyE1S9SLKBwIYZgRg5CKpehG3EwEYdgRg5KLT1DMjXwAgACMHSVPPjHwBgACMjCVNPVO9CAACrIJGZpIyXVG9CAAOYgSM1OLKBnZTPpCpZwA4iBEwEnUqG5i2fCBTzwBwKAIwOoqu6aZJI9kOU88AcDgCMNrqt3qRxP2+ANAO14ARq9/qRRL3+wJAJwRgxErK4ZyE8oEA0BkBGIeIVjp3SqQhHV42kPKBANAdrgHjgDQLrqheBADZYASMA5KmnbmmCwDZIQBDUnL+5okJrukCQJYSp6DN7CZJ75G0w91XhNuOk/QFSZOStkn6E3d/Lr9mIk9p8jdv2zaw5gDAUEgzAr5Z0vkt266SdI+7v1HSPeHvqKik0oEk0QCA7CUGYHe/V9Luls0XSbolfH6LpPdm2ywMCqUDAaAYvV4Dfp27PytJ4ePx7d5oZqvMbNbMZufn53s8HPJA6UAAKE7ui7Dcfb27T7n71Pj4eN6HQ0pRmkmmngGgGL0G4J+Z2YmSFD7uyK5JyAKlAwGg3HoNwHdKuiR8fomkr2bTHPRrZkYaGwsCbHRtl9KBAFA+iQHYzG6TdJ+kU8xszswuk3StpPPM7DFJ54W/o2DRNd3mur3dYuoZAAYj8T5gd7+4zUvnZtwW9IHSgQBQLWTCqgFKBwJA9RCAa4DSgQBQPQTgiktKpCFROhAAyohyhBWWlEiD0oEAUF6MgCssKYczwRcAyosAXFHkcAaAaiMAVxA5nAGg+rgGXDFJ9/uSSAMAqoERcIWkud+XqWcAqAZGwBWRJtMVU88AUB2MgEuuubhCp+DL1DMAVAsBuCTiygdGZQOTiiuQwxkAqocp6ILNzEirVx8aZNOWDZSCkS/BFwCqhwBcoGhRVa95nBn5AkB1MQVdkGhRVa/Bl0xXAFBtBOAC9Fs+kOpFAFB9TEEPWJrbidppNKR16wi8AFAHjIAHKO3It7V8IGUDAaB+GAEPUKfqRRLlAwFgmDACHpCk6kUsqgKA4UIAHoCk6kXcTgQAw6evAGxmV5rZw2a21cxuM7PFWTWsTjpNPTPyBYDh1HMANrOTJP2VpCl3XyFpRNL7smpYXSRNPTPyBYDh1O8U9EJJR5jZQkmjkp7pv0n1kTT1TPUiABhePQdgd39a0vWSnpT0rKSfu/vG1veZ2SozmzWz2fn5+d5bWjFJma6oXgQAw62fKehjJV0kabmk10s60sxWtr7P3de7+5S7T42Pj/fe0pKJq17UWsWo0/2+TD0DwHDrZwr6dyU94e7z7r5X0pck/XY2zSqv5vq80bXdKNCmrWLE1DMAoJ9EHE9KepuZjUp6QdK5kmYzaVVJ9Vu9SGLqGQAQ6DkAu/v9ZvZFSZsl7ZP0gKT1WTWsbPrJ4Rzhfl8AQKSvVJTu/lFJH82oLaXVb/UiKRj5EnwBABEyYaWQlMM5CeUDAQCtCMAJkhJpSIdXL6KKEQAgCdWQOkiTw5k0kgCAXjACbiNNIg2CLwCgVwTgGGkWXXFNFwDQDwJwjKRFVyTSAAD0iwDcImnRFYk0AABZIAA3SbPoiqlnAEAWCMAhFl0BAAZpKAJwp8pFY2PSa15D9SIAwGDV+j7gmRlp9Wpp166D21orFzW/1g6LrgAAWattAM6icpHEoisAQD5qGYCzqFwksegKAJCf2l0DzqJykcSiKwBAvmoVgJNWMqdF9SIAQN5qE4DTjnxbKxc1GsGPGdWLAACDU5trwEnpI6lcBAAok1qMgNOkjyT4AgDKpPIBmPSRAIAqqnQAJn0kAKCqKhuAqdkLAKiyvgKwmR1jZl80sx+Y2aNmdnZWDeskze1GpI8EAJRZv6ug10m6y93/2MxeJWk0gzZ1lGbkS/pIAEDZ9RyAzexoSW+X9KeS5O4vS3o5m2a1l+Z2I6aeAQBl188U9BskzUv6nJk9YGb/YWZHtr7JzFaZ2ayZzc7Pz/dxuMCTT7Z/jUVXAICq6CcAL5R0pqR/d/czJP1a0lWtb3L39e4+5e5T4+PjfRwusGxZ/HZGvgCAKuknAM9JmnP3+8Pfv6ggIOdq7dpgpNuMkS8AoGp6DsDu/lNJT5nZKeGmcyU9kkmrOpieDka6ExMH8zcz8gUAVE2/q6A/JGkmXAH9uKQ/679JyaanCbgAgGrrKwC7+xZJU9k0BQCA4VHZTFgAAFQZARgAgAIQgAEAKAABGACAApi7D+5gZvOStme4yzFJOzPcX5HoS/nUpR8SfSmruvSlLv2Qsu/LhLvHZqEaaADOmpnNunstVmHTl/KpSz8k+lJWdelLXfohDbYvTEEDAFAAAjAAAAWoegBeX3QDMkRfyqcu/ZDoS1nVpS916Yc0wL5U+howAABVVfURMAAAlUQABgCgAIUFYDNbambfMrNHzexhM1sdbj/OzO42s8fCx2ObPnO1mf3YzH5oZr/ftP0sM/t++NonzMzaHDP282Xpi5mNmtk3zOwH4X6ubXO8STN7wcy2hD+fKVtfwu3fDrdF7Ty+zTHL/r0c1dSHLWa208z+LeZ4uXwv3fbDzBrh+39lZje07KtS50q7vlTxXEn4Xip1rnT4Xqp2rpxnZpvCc2KTmb27aV/5nyvuXsiPpBMlnRk+P0rSjyS9WdJ1kq4Kt18l6WPh8zdLelDSqyUtl/QTSSPha9+TdLYkk/RNSX8Qc7y2ny9LXySNSnpX+J5XSfqfNn2ZlLS1At/LtyVNJRyv9N9LzH43SXr7oL6XHvpxpKRzJP2lpBta9lW1cyW2LxU9Vzp9L1U7V9r2pWLnyhmSXh8+XyHp6UGeK4WNgN39WXffHD7/paRHJZ0k6SJJt4Rvu0XSe8PnF0m63d1fcvcnJP1Y0lvN7ERJR7v7fR78RT7f9JlmsZ8vU1/cfY+7fyvcz8uSNktakkUb08qqL10csvTfS/M+zeyNko5X8D/8gei2H+7+a3f/jqQXm/dTxXOlXV+qeK6060sXSv+9NKvIufKAuz8Tbn9Y0mIze/WgzpVSXAM2s0kF/xK5X9Lr3P1ZKfhjKvgCpeCP+FTTx+bCbSeFz1u3t2r3+Uz12Zfm/Rwj6QJJ97Q51HIze8DM/tvMfiezDhzahkn135fPhVNMf99mCqdS34ukiyV9ITwp4+T6vaTsRztVPFfS7OcYVeNcSVKlcyWNqp0rfyTpAXd/SQM6VxamfWNezOw1kv5L0hXu/os20+xSMA3QyjtsT/v5zGTQl2g/CyXdJukT7v54zHuflbTM3XeZ2VmSvmJmp7r7L/rrQVMDs+nLtLs/bWZHhft6v4J/Sab9fCay+l5C71PQjzi5fi9d9KPtLmK2lf1cSdpPlc6VTqp2rqRRmXPFzE6V9DFJvxdtinlb5udKoSNgM1uk4I804+5fCjf/LBz+R1NmO8Ltc5KWNn18iaRnwu1LYra3avf5TGTUl8h6SY+5+7/FHSuc7tgVPt+k4LrDyRl1JbO+uPvT4eMvJf2n4qdmKvO9mNlpkhaGf/PD5Pm9dNmPdqp4riSp0rnSVgXPlaR9VeZcMbMlkr4s6QPu/pNw80DOlSJXQZukGyU96u7/0vTSnZIuCZ9fIumrTdvfF87PL5f0RknfC6cTfmlmbwv3+YGmzzSL/XyZ+hLu6x8lvVbSFR2ON25mI+HzN4Sfj/vXf2F9MbOFZjYW7nORpPdI2hpzyEp8L6GLFYy22h0vl++lh37Equi50mlfVTtX2u2niudKkkqcKxZcvviGpKvd/X+jNw/sXPEcVgem+VGwgs4lPSRpS/jzh5IaCq7lPBY+Htf0mTUK/qX0QzWtSJM0peA/2J9IukE6kOHrQknXJH2+LH1R8K8nV7BwINrPn7f2RcG1iocVrL7bLOmCEvblSAUrIB8K27pOB1dHV+p7aXrtcUm/2bIt9++lx35sk7Rb0q8U/Cv9zRU+Vw7ri6p7rsT1parnSux/Y1U6VyT9naRfN713i6TjB3WukIoSAIAClGIVNAAAw4YADABAAQjAAAAUgAAMAEABCMAAABSAAAwAQAEIwAAAFOD/AWuolsNCaTpQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig3 = plt.figure(figsize=(8,4))\n",
    "plt.plot(t,yn3,'bo:',label='RK4')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Tq5oFE2zbJCm"
   },
   "source": [
    "The graph shows the behaviour of the solution; as we might expect the overal population increase that we witnessed in case 2 has been modified by a shorter oscillatory pattern.\n",
    "\n",
    "This emphasises the power of RK4: we have used it to solve an equation which has **no simple analytical solution**. Given our previous test case, we can be confident of a high degree of accuracy, despite having no solution to compare it too."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tzflYVlOqTR9"
   },
   "source": [
    "### Results Table\n",
    "The table below shows the time and RK4 solution, $yn3$,  for the non-linear population equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 363
    },
    "executionInfo": {
     "elapsed": 7,
     "status": "ok",
     "timestamp": 1646995515418,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "Mj0ZTItPqTR-",
    "outputId": "e9a175c3-c66c-44a5-f0f3-fbf053cc5a9d"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time t_n</th>\n",
       "      <th>Runge Kutta</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2000.0</td>\n",
       "      <td>6.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2000.1</td>\n",
       "      <td>6.114813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2000.2</td>\n",
       "      <td>6.279909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2000.3</td>\n",
       "      <td>6.465115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2000.4</td>\n",
       "      <td>6.632805</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2000.5</td>\n",
       "      <td>6.752294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2000.6</td>\n",
       "      <td>6.811571</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2000.7</td>\n",
       "      <td>6.821875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2000.8</td>\n",
       "      <td>6.813379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2000.9</td>\n",
       "      <td>6.823656</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   time t_n  Runge Kutta\n",
       "0    2000.0     6.000000\n",
       "1    2000.1     6.114813\n",
       "2    2000.2     6.279909\n",
       "3    2000.3     6.465115\n",
       "4    2000.4     6.632805\n",
       "5    2000.5     6.752294\n",
       "6    2000.6     6.811571\n",
       "7    2000.7     6.821875\n",
       "8    2000.8     6.813379\n",
       "9    2000.9     6.823656"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'time t_n': t[0:10], 'Runge Kutta':yn3[0:10]}\n",
    "df = pd.DataFrame(data=d)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cH-1Bl9rMOxR"
   },
   "source": [
    "## Python ODE solvers\n",
    "\n",
    "Python is a clever bit of kit, with entire libraries dedicated to solving numerical problems like these. Some of these libraries even come with built in solvers to carry out the tasks we have just performed. If you feed a particular function, initial/final values, and even explicit grid-size to these solvers, we can obtain automated solutions to our problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "id": "7J6X6wCBMjzT"
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "sol = solve_ivp(nlo_fun, [t_start, t_end], [6.0], t_eval=t)\n",
    "sol2 = solve_ivp(nlo_fun, [t_start, t_end], [6.0],  method='DOP853', t_eval=t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "executionInfo": {
     "elapsed": 345,
     "status": "ok",
     "timestamp": 1646995616115,
     "user": {
      "displayName": "James Threlfall",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh5AMmwJn-QgnyAtyVfjDwEfQL-5FfwAhQ1uZjJ=s64",
      "userId": "12893452708265254015"
     },
     "user_tz": 0
    },
    "id": "eu2IpzNtN5kW",
    "outputId": "65b1915d-96d7-4670-e311-b50fa56568db"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD4CAYAAAA0JjXXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyAElEQVR4nO3dd3xVRf7/8dckEEKvoRk6ShEhQERBRRQpKog/dRVUEEUQFF2woX5d22JZRVkFd1lcBRUFO6KigoALqJSAAUVUWoAASkJvIW1+f8wNN4Q0kptbkvfz8cgjZ849d85Mbk4+mXOmGGstIiIi4l9hgS6AiIhIWaQALCIiEgAKwCIiIgGgACwiIhIACsAiIiIBUM6fJ6tTp45t2rSpP08pIiISMKtWrUq21kbl9ppfA3DTpk2Ji4vz5ylFREQCxhizNa/XdAtaREQkABSARUREAkABWEREJAD8+gw4N2lpaSQmJpKSkhLoopR6kZGRREdHU758+UAXRUSkzAt4AE5MTKRq1ao0bdoUY0ygi1NqWWvZs2cPiYmJNGvWLNDFEREp8wJ+CzolJYXatWsr+JYwYwy1a9fWnQYRkSAR8AAMKPj6iX7OIiLBIygCsIiISEnYuhWmTIFdSzbCRx9hk5IDXaQTFICBbt26Ffm94eHhxMTE0K5dO/r378/+/fsBSEhIoF27dieOe+211+jUqRP79u07sW/ChAkYY0hOTj7xnooVKxITE0NMTAwjR44scrlERMqivXvhqadg69Lt0K8fmxZuZdQosF/MhYEDmf/pUTp2hPXrA11SBWAAvv/++yK/t2LFisTHx/Pzzz9Tq1YtXn311VOOefvtt5k0aRLz5s2jZs2aAGzfvp358+fTuHHjk45t0aIF8fHxxMfHM2XKlCKXS0QkVGVmwsGDcDpdVqx13/ftg7//HdYuOworV9KtcSLbt0O9h26FdevIjG5M06Zw1sSR8PTTrFoFu3aVSDUKpAAMVKlSBYAbbriBuXPnntg/dOhQPvroI6ZPn86AAQPo27cvrVq14sknn8w1n65du7Jjx46T9r3//vs899xzzJs3jzp16pzYP3bsWJ5//nk9lxURAX791d0uBki56nr+W/1e3nnHpdPT4bPPXGDOzejRcOedQEYGLTbNY+tW6H9/K0hIILLnBURHQ3iNqnDWWfTtC598mEH4kUPYw0e47TYYMABITIS0NH9U9YSgC8A9esD06W47Lc2lZ8xw6aNHXfq991z6wAGX/vhjl05OdunPPnPpP/44vXMPHDiQ9zyZp6amsmDBAq644goAVqxYwTvvvEN8fDwffPDBKXNaZ2RksGDBAq666qoT+7Zu3cro0aOZN28e9evXP7F/zpw5nHHGGXTo0OGUMmzZsoWOHTty8cUXs2TJktOrgIhICDq0P4NxHb7iH/9w6XIx7bi4O3TpAljLzktu4t9XzeXrr93rGzbAc89531+pElSuDHbqa9CnDw13ev4+V6yY+wnDw+GddzBPj+ejj+CNUSuhUSOYP7/E6pibgI8DDiaXX34599xzD8ePH+err76ie/fuVPR8gL169aJ27doAXHPNNSxdupTY2FiOHTtGTEwMCQkJdO7cmV69ep3ILyoqilq1avH+++8zduxYAI4ePcrTTz/NvHnzTjl/gwYN2LZtG7Vr12bVqlVcffXVrFu3jmrVqvmh9iIiRZeSAs8+C336QGG61bz2GsTHw6uvQtWZU/k09U729l8GnEfE+MfonHVg8h4a7Y3n6bsupUMfIDOTD943PP6EYfhwqF3L8vwTx1wUPn4r1KgOsbGFK3RYGC1bAqYWTJ0KMTFFqHnRBV0A/vZb73b58ienK1U6OV29+snpOnVOTmdrdBZKZGQkPXr04Ouvv+a9995j0KBBJ17Leas4K531DPjAgQP069ePV199lXvuucdT3kp8+eWXXHjhhdStW5ebbrqJTZs2sWXLlhOt38TERDp16sSKFSuoX78+FSpUAKBz5860aNGC33//ndjC/jKJiPjJ8ePwl2sz+UvthQweUZEK3S5gwgSoVC6Vbt0iyMyERYuga1f3t3vvXvj+e7jySjAGtm2D9b9Y0tIM5YcNg3r1qNW3y6knqlMHs2YNHY1x92xnvMsDs15kzJ6VVKpWDobcAj/9BMuWQYUKkO3vdqG1aOG+/CzobkEH2sCBA5k2bRpLliyhT58+J/bPnz+fvXv3cuzYMWbPns0FF1xw0vuqV6/OK6+8woQJE0jL9hwhKiqKr776ikceeYSvv/6ac845h927d5OQkEBCQgLR0dGsXr2a+vXrk5SUREZGBgCbN29mw4YNNG/e3D8VFxEpwPHjLoiCi3Vn1Uziphl9Yc0ajIF9q7fw4MQG8NVXfPYZXHYZrFjhjp83D/r3h5UrXfqJC+azsFxvyh8/DBERcM01LjLnplw5d9sYoGpVylev7IIvwLXXwh13gOdvZyhRAM6hd+/eLF68mMsuu4yIiIgT+y+88EIGDx5MTEwM1157ba6t0o4dO9KhQwdmzZp10v5mzZoxZ84cbrvtNpYvX57nuRcvXkz79u3p0KED1113HVOmTKFWrVq+q5yIlHkZGZCU5E1v2QLHjhXuvePGQf+eRzk07UMAJrxdj7BFC+G22wCICM/A9OkN7drRpw/Me3MXHRvvAaBnT1iyBDp1cnmFh1nXcSfb0MxCGTAAli49OT1ypGtmhxhjs/pu+0FsbKzN2Xlp/fr1tGnTxm9lKIrp06cTFxfH5MmTA12UYguFn7eIlJwePdz3rMd1HTq4x3cLFrj00aPeWLZ1Kzz5JDzyCLRsCRs3ghn/d5q//QTm11/hzDPzP9n/+3+wfDls3+5asOPHQ1iYyxDcfwNZLdtSyhizylqb63PEoHsGLCIivnP0KLz7rmukhu3+g7vvqkt6Zhj8+CN88w3/GDeMStHuTtvx41C1KkyeDKNGuU7Es2fDJZe4ANyyJTB1HNzaveDgC25A7u+/e4PswoUnd3Qq5cG3IGoBlzH6eYuULTNnwo03wtqnP+Ocx69x4zT79oUPP4QbboDNm6FJE8jM5PDRMCZPho4dXW9m8DRSF86HF16ATz/Ne2iP5EotYBGRELZoEbz0Erz8MjRv7iakCMujB09GBvznPxAVBX/5C1x/PTRrBufEXg7ln4HWrd2B113nuiKfcYZLDx9OlbAwHnrtNW9m1hIeblzTePdu15U563gpNnXCEhEJQjt3uukYwd1G3rABGjRw6YkT3fCerBuY8+d7JzAKC4Np09ytY1avJrxvL85vtc/1JH7gAWja1HuSrGBqLTRsCHXrel9r0ACeecZt9+sHq1Yp+PqYWsAiIkEmOdk1VEeOhOefhysuTeHKZSkQWR0wNGzo5jnIGrUzY4brVDV0qNs3bx7UqAEs2Ot6Uu3dC5556HNljHtem92IEScH3DL+vLYkqAUsIhIErHXzIYPrlTx53HZG3uGauOalF6FWLdcsBgYNtHzyife9L78Mv/ziTdfcs9EF58sug3XrijbJxJNPuoguJUYBGN8sR3j22WfToUMHXnrpJTKzzRi+dOlSunTpQuvWrWndujVTp0498doTTzzBGWeccWI5wzlz5gCwbds2LrnkEjp27Ej79u1PWiAi63wxMTEnzTs9bNgwOnToQPv27bnuuus4fPhwkeskIr7x559uIYHCeOEFNyRoyxZgyRKG/P1Mmm//n3uxVy944w1vi3TUKM/qA06NGm4uZAA++gjatnW9nMFNKSjByVqb7xfwBrAb+DnbvhhgGRAPxAFdCsrHWkvnzp1tTr/88ssp+0JJ5cqVT2z/+eeftmfPnvaxxx6z1lq7a9cu26hRI7tq1SprrbVJSUm2U6dO9vPPP7fWWvv444/bF154wVrrfg61a9e2GRkZdvjw4fZf//qXtdbadevW2SZNmuR6vuwOHDhwYnvs2LH22WefzfW4UP95iwSzhQutvf9+b7pbN/eVJSHB2uPH3fbBg9a+9JK1K1e69I4d1k6aZG1amrX22DFrx4xxO3PKzLT2gQesHTfOu+/KK6198023nZTkXktL82ndpGiAOJtHTCxMC3g60DfHvueBJ621McBjnnTI8tVyhHXr1mXq1KlMnjwZay2vvvoqQ4cOpZNn6pc6derw/PPP81z2ZTw82rRpQ7ly5UhOTsYYw0FP74sDBw7QsGHDAuuQtWCDtZZjx45pmUORAPj8c/jiC++yeX/9Kzz6qNu21rVwPeuyEBYGDz8MX33l0g3rZTA69SXKpadAZKTraZXbtW+MezCc/e/Inj2upxa4+9fPPec6XUlQKzAAW2sXA3tz7gayluipDuz0WYkCuB5hcZYjzNK8eXMyMzPZvXs369ato3Pnzie9Hhsby7p160553/LlywkLCyMqKoonnniCGTNmEB0dzRVXXMGkSZNOHJeSkkJsbCznn38+s2fPPimPW2+9lfr16/Prr79y9913n1bdRaRovvwSfv7ZbY8fD3Fx3iFC118Pl1/utjMz4ZVX4OabXbpyZTcKKCtAs3Qp3H+/9+/Z6fjhB9djS0JKUZ8BjwFeMMZsByYAD+d1oDFmhDEmzhgTl5R9AtIgdPnll7Nw4UKOHz/Ol19+metyhBUrVjyxHGFerGdsgLU215Zo9n0TJ04kJiaG+++/n/feew9jDDNnzmTo0KEkJiYyd+5cBg8efOK58rZt24iLi+Pdd99lzJgxbNq06URe06ZNY+fOnbRp0+bEPxIiUnIOH4ZbbuHEOrYVK0Klz9+HK67wNoOPHwdcJ+IhQ9zwoSx16+IdS3Txxe657Y03+q8CElBFDcCjgLHW2kbAWOD1vA601k611sZaa2OjoqIKzjmrLz141yPM+pcxaz3CG25w6az1CK+5xqWz1iPs39+lT3M9wpzLEQ4cOPDEa3ktR5jT5s2bCQ8Pp27dupx99tmntJRXrVpF27ZtT6THjh1LfHw8S5Ys4aKLLgLg9ddf5/rrrwega9eupKSkkJycDHDidnTz5s3p0aMHP2Z1tPAIDw/nhhtu4KOPPjqtuouUZRkZ3jG3BUlPd/2crIUqVeCbb+C/j2z29rY6fNitdpCa6tIPPujmcsxttZ71693atVnXsWeZUikbihqAbwGy7pN8AOSyiGNoKupyhABJSUmMHDmS0aNHY4zhrrvuYvr06cTHxwOwZ88exo0bx4MPPphvGRo3bswCz8zo69evJyUlhaioKPbt28dxz3/TycnJfPfdd7Rt2xZrLRs3bgRcq/uzzz6jddZsNyKSL2vd1IvZl5H96qsTI35OHJPVoP3wQzeJ1Pz5Lt0+ahcVunaCv/3N7bjtNrfmXmSkS8fEQJcu3nG0Dz8Mr3vaLPXru8yPHCmx+knwKupT+p3AxcC3wKXABl8VKNB69+7NkCFDuOqqq3JdjnDjxo3ceOONJ5YjPHbsGDExMaSlpVGuXDkGDx7MvffeC0CDBg2YMWMGw4cP59ChQ1hrGTNmDP2zWuh5ePHFFxk+fDgTJ07EGMP06dMxxrB+/XruuOMOwsLCyMzM5KGHHqJt27ZkZmZyyy23cPDgQay1dOjQgX//+98l90MSCXHWukeuF13k+jTdfbd3gqhjx9xNtTvvhAkTXNeTmjXhxRdh9GgXfKtVcyODADdj1JNPwtVX536yW2/1bqeluYfG/fq5dM2aboYpdZoskwpcjMEYMxPoAdQB/gQeB34DXsYF8BTgTmvtqoJOFqqLMWg5QpHSZdYs1+KdNy9bIPWw1jupxdlnu/Sjj0K3bnDlldkO/Owzd0Dz5n4rt4SeYi3GYK0dlMdLnfPYLyISlFJTISICrr3W3QW+9NJTjzHGxdXs6aefznFQSgrccYdrQqvDoxSRZsIqhKFDh5aK1q9IWTZpEnTu7B63li/vHtWG79jmIvFez0jL48fdbeKCREa6B8XZVw4SOU1BEYALug0uvqGfs5RFWb/2bdu6/lDZZoqFNWvg9tvht99c+sMPoV49t/RQTikpbn2/rP4V7du7h8EiRRTwABwZGcmePXsUHEqYtZY9e/YQmdUzUyQEJSfDBx9402+8AY8/7g2y2aWnw4ABbh1dgJ494e23oeqeBO/0U/36we+/u6FAAGed5cbhZi1e8PzzMGyY246MdCfat69E6iZlT8DnKouOjiYxMZFgn6SjNIiMjCQ6OjrQxRApsvHj3WxSF17oOh/Hx8Py5a4TMsDAga6hOnu2m4mxcuUcLV6AMWNgxQrYtMnNnHHmmd7Xzj3XfWVJS3MrKmT54AP1WBafKbAXtC/l1gtaRCQ/GRmwfz/Uru2C68aN0K7dya9nDbF98UXX8h03Lp8Mk5Nh9253T1qkhBWrF7SISCDddJNrrH7/vbsLnD34wsnrxN93Xx6ZbN0K//oXPPOMmzGvTp0SK69IYQX8GbCISH4GDYIRI4q5uM8XX8DUqZCQ4KtiiRSbbkGLSFCxFv7zH6hVy60mdJL0dDfXco0aLp2Wlv+C88eOuee8ADt2eBe0F/GT/G5BqwUsIkHFWnjnHXj//RwvHD8OrVp57zNbC82aubmVs785y5tvulmqPAuZKPhKsFEAFpESl5LiWrNff+3S6elu7ousBYJ++snNsbx3r1tLd/bsbAE4K6hWqAAPPeSmsQLX+h061C3jBy7Q1qsHCxe69LnnQt++fqidSNEoAItIiVi5EqZPd9sHDsDatW6VPoDNm12v5lmzXDozE5YscavzgXstLAzX/fnSS8GzOhjDh7u1dsHNKTl+vDfIpqS4cb2eFcNo2xamTVOHKwla6gUtIj5nrRsKlJrqGqn16sGvv3pfr1ULJk70zn/RoQMkJubyODcjAw4dcs99CxId7WbmEAkR6oQlIj6TmenuDFeo4OavCAuDqKgiZJSW5ro9G3PyQF+REKNxwCJS4jIz3XPcSpXclI/16hUjoyFD3DRWr72m4CullgKwiPhEWBh06eICcLEYAy1bQpUqmvZRSjXdghaRYvniCzfCJyammBkdPuweBLdu7YtiiQQFjQMWkRKRkgJ33gl/+1suL27d6uaR/P57l/79dzc0aPFil875z/+AAe4rPb1EyywSLBSAReQkv/8Of/+7W/QA3GRSW7a4Hs3ghhPde6+Lk5GRMG+eW0b3FJUqwbJlbuEDcNG6ShXvw+FPP3U9tPbudemnnoLXXy/mnJMioUMBWEROkpICjz3mGrAAS5e6CaVWr3bpTZvcVJFr17p0q1au1zPgWrWffOK+R0W5he6vvtq91r49LFrk3gBQv757bf9+l77gArfOoEgZoWfAIsKGDW6ui5EjXfroUde6DQuD7dth/nw3x0Xdut6WcERELhl9/LGbqerTT+Gqq/xWfpFgld8zYAVgEWH0aDf146+/ukkyisxamDPHBV/1YBYpXicsY8wbxpjdxpifc+y/2xjzmzFmnTHmeV8VVkT8J+v/75deghUrihF8v/4a/vjDBd0BAxR8RQqhMM+ApwMnzWhujLkEGAC0t9aeDUzwfdFEpCR9+CH06OFG/0REQNOmRczoyBG4+WbXM0tECq3A7obW2sXGmKY5do8CnrPWHvccs7sEyiYiJShrpsdij/qpXBm+/LIYEVykbCpqL+izgIuMMcuNMf8zxpyb14HGmBHGmDhjTFxS1lIoIhIQx4+7kUHgOiAvWuRd2/60LVrkHX8UG6tVh0ROU1EDcDmgJnA+8ADwvjG5P/Sx1k611sZaa2OjijQru4j4ygMPQM+e3qG5J121GRnwv//Btm0ufegQPP+8d0BwThMnuvFKx46VaJlFSquiBuBE4GPrrAAyAf37KxIAq1d7J5sCN2nGkSNu+9gxeOstt/4uuPXsZ81yw4lOsX49XHYZvPOOS2/b5tYUXLPGpX/6CZo0cd/BrbW7ciVUrFgi9RIp7YoagGcDlwIYY84CIoBkH5VJRArp2DHo3Rvee8+lrYV27dxMVlmGDYNJk9x2w4bQv3+OTLImwmjXzt2fHjXKpc8+G5KT4fLLvZl37ep9aFy7tnv+KyJFUmAnLGPMTKAHUMcYkwg8DrwBvOEZmpQK3GL9OaBYpIw7cACqV3eNzzlzoEULtz8z081S1by5S1es6Bq2zZrlkdHy5dCnj4vgffpA584nv167tne7fXvXfBYRnyhML+hBebx0s4/LIiKFsG0bnH++WwBh1Cjo1s37Wni4GxGUXcuW+WR21lmuN1anTiVRVBHJh2Y9FwkxDRq4iaYuvrgYmRw86BZGqFkTpk/3VdFE5DQoAIuEiKNH3fdKlWDKlGJklJ7uJnZu1Mjb4UpE/E4BWCREDBkCO3e65XSLtWJfeDhcd51brUhEAkYBWCRE3HyzC8DFCr779rnbzvfc47NyiUjRaD1gkSCWmekddnv11XDnncXIbOpUaNPGOyhYRAJKAVgkiD35JJx3HiQk5HNQRoZ3O7fRgFmvX3yxW6mocWNfFlFEikgBWMTP9u2DvXvzft1aNwskuBbvP//pJqDK1bhxcG62qdhvv93NZpXl6qvdw2OAVq3cIOFi3cMWEV/RlSjiZ8OGuQmndu506ddfd9t/+5tL9+8PFSrARx9BvXowYkSODNLSvEsZtW/v9qWnu32xsa53c5bzztNUkSJByvhzAqvY2FgbFxfnt/OJBANr4ZVXYOBAF1DXrIG1a2HwYPf6kCGwfbtbXAjgxRehVi249dZcMjt+3E0NedVVMGaMv6ogIkVkjFllrY3N7TW1gEVK2JYt8Mgj7rbyo49Chw7uK8tbb518/H335ZOZMRAd7Xoyi0hIUwAWKWHNm0NcnHsEWyyZmRARcWrEFpGQpE5YIiXkxx/dQgngRv+EFedq++ADuOACSErySdlEJPAUgEVKyDPPwOjRPlqvvlw5qFbNLYEkIqWCOmGJlJCUFLdy0VlnFSMTa91z35zbIhIS8uuEpRawiI99/bXrrBwZWczgm5rqFk349FOXVvAVKVUUgEV8KCHBxcynn/ZBZocOwZ49cPiwDzITkWCjXtAiPtS0KcyeDT16FCOT3bvd897ateG779zqRSJS6qgFLFJM1sIbb7hlAgGuvBIqV87nDS+/DK++6k0vWwbLl7vtlBQ4+2wYO9alFXxFSi0FYJE8ZF/jICEB5s6FI0e8+7L6L6akwHPPwWuv5ZNZerp3e9kyN89kljFj3Awd4B4cP/oo3HVXMUsvIsFOAVgkFwsWQIMG8NtvLj13rmvZHjjg0q++6m43p6W5qZYXLYLp0/PI7McfoXVriI936RkzYOFC7+vTpsGkSd70X/8K7dr5tkIiEnT0DFgkF23bQvfubvIpgEGDICbGBWWAli3dogn790NUFJxxRj6ZNWjgFkiIiHDpnLeV27TxcelFJBQUOA7YGPMG0A/Yba1tl+O1+4EXgChrbXJBJ9M4YAl2u3ZB/fo+GvFz6BBUreqDjEQkVBV3HPB0oG8umTYCegHbilU6kSCRnAydOnmXBSyWI0fgwgvhoYd8kJmIlEYFBmBr7WIgt+XDJwIPAv6bSkukBNWsCXffDTfc4IPMKlRwD40vucQHmYlIaVSkZ8DGmKuAHdbaNaaAe3XGmBHACIDGjRsX5XQiJS4jwz2afeQRH2RmrZu7+ZlnfJCZiJRWp90L2hhTCfg/4LHCHG+tnWqtjbXWxkZFRZ3u6URK3MqVrtPx+vU+yGztWujSxdt9WkQkD0UZhtQCaAasMcYkANHAamNMfV8WTMRfMjJcT+b6vvgNPnLEjfmtVs0HmYlIaXbat6CttT8BdbPSniAcW5he0CLB6PzzvbNYFVl6urvt3LUrrF6thRNEpEAFtoCNMTOBH4BWxphEY8ywki+WSMkbMwYmTPBBRnv3wgUXuAk1QMFXRAqlwBawtXZQAa839VlpRPwkPR127nSN1iLL6rlVtSo0bAi1avmsfCJS+mkqSilTkpLcV7lyMGsW/OMf+Ry8ZYt3Qui33nLzNG/c6NLTprmHxgcOQPny8MknMGBAiZdfREoPBWApVZ56Cvr08aYXL3bzOIObVjI2Fq6/3o0UCgvLZ7Ghb76BFi3g889d+pxz4J57oF49lz7zTBg82M12JSJSBJoLWkLepk3QvLl79Fqnjne+ZoCJE93woiuucAF38mQXVwt8THvJJe4B8bnnunTHju4ry4UXui8RkSIqcC5oX9Jc0OJrq1bBRRe5W8l3333q63v2uLvEzZsXMsMpU1wTWc9zRcQHijsXtEjQ6tgRxo93MTM3tWufRvDduNEtBThlis/KJyKSF92ClpD0zTfQoYObQOPee32UacuWEBen5QFFxC/UApaQc+iQWzDh/vt9lOHBg/Ddd277nHOKOTZJRKRw9JdGQk7VqvDll9CkiY8yfOopmDQJNm+GM87wUaYiIvlTC1hCRkaGu0MMbr2DrBFBxfbEE/Dhhwq+IuJXCsASMiZNgvPOcwsO+cSWLW5KrCpVoH9/H2UqIlI4ugUtIeP2212sbN/eB5kdOQIXXww9erhZrkRE/EwBWILesWMQEeGC7+23+yjTypXh6aehbVsfZSgicnp0C1qC3l//Cpde6u4WF4u1MG8erFzp0oMHQ+fOxS6fiEhRqAUsQe+ii6Bx40KODtqxw63Hm/VMd9EiSEiAW2916QcfdBmtWOHmphQRCRBNRSmly/Dh8N57bmwvuHvWs2dDcrJLb9oE0dFQoULAiigiZYemopSQ9OyzLpYWaMUK+PNPt/3YY24JpKx/LF9+2Y3vzdKihYKviAQFBWAJSmlpruG6cGEBBx45Av36wejRLt2oEcTEeJc7qlwZqlUrwZKKiBSNngFLUCpfHr7/HlJSCjiwcmX44ANo3dov5RIR8RW1gCXo/PCDG3oUHu7ia64OH3a3msGN5/XZtFgiIv6hACxB5cAB6NsX7ryzgAOfegouuwy2bfNLuUREfK3AW9DGmDeAfsBua207z74XgP5AKrAJuNVau78EyyllRPXq8OmnrqNyvh57DLp3d+OTRERCUGFawNOBvjn2zQfaWWvbA78DD/u4XFIGZU200aOHW5o3V3v3eudv7tfPX0UTEfG5AgOwtXYxsDfHvnnW2qx5iZYBBbVXRPJ15Ah06ADTp+dzUGYmXH21C7x+HL8uIlISfNEL+jagMKM1RfJ09Ci0auWG6eYpLAxGjXJDjLKGGYmIhKhiBWBjzP8B6cA7+RwzAhgB0FjP6yQPUVHw8cf5HHDggHtAPGiQ38okIlKSitwL2hhzC65z1k02n/ksrbVTrbWx1trYqKioop5OSqn0dHjkEe9EVrn68kto1gw0jamIlCJFCsDGmL7AOOAqa+1R3xZJypK4OHjxRe+Q3lx17AhXXplPzywRkdBTYAA2xswEfgBaGWMSjTHDgMlAVWC+MSbeGDOlhMsppdT558OGDXBdu1/d+KOsrtCffAJ/+YvrbFW/Prz9NtSoEdCyioj4UmF6QQ+y1jaw1pa31kZba1+31ra01jay1sZ4vkb6o7BSeth9+9n1wEuQmkrjxmA+/8z1cE5NdQds2gS7d7tnvyIipZBmwpKAWPbYXBpMuI+1//SstjB0KKxaBZGRLn3//fC//6nVKyKllhZjEP86fhwqVKDdMzfynxrduf0+zxDyqCj3JSJSRqgFLP7z++/YVq1I//EnqlaFO/4eTXh4oAslIhIYCsDiP5mZ7CrfmOuGVWffvkAXRkQksBSAxX9at2bZPxZTvkVjPdoVkTJPAVhK3gcfwMMPQ2oq11zjkppJUkTKOgVgKXF2+Qo2vb6Ijz7Rr5uISBb9RZQSd+ixF7it6SK2bFenexGRLPqLKCXnzTexF3WnWvNmzFtSkYiIQBdIRCR4qAUsJePAAVLuuo9FvZ8hIwMqVNBzXxGR7NQClpJRvToz7lnBV6ui6JaGxvuKiOSgFrD41N4FP7L1kf8AcPszzXlvbtUTs0uKiIiXArD41P9GzCD8hWdJ338YUMtXRCQvCsByWlbP38Nn7R/hwDcrAZj3r42sL9eOg7PmAnDGy+M4PHcJ5WpUCWQxRUSCngKwFGhXYgYpGxMBSNoTRt+fXuDId/EAREbXITW6BekRlQDo0q8urXs1ClRRRURChrHW+u1ksbGxNi4uzm/nk+I7cADmNhhGbM2NnJn4revKnJ4O5dR/T0SkIMaYVdba2Nxe019RyVf16lDl5qupWmmddxyRgq+ISLHpL6nkavVqqHZkFy0vakD/qf2B/oEukohIqaJnwHKKzEx465rZRF/cHPv9D4EujohIqaQALKcIC4Mx73fj6C13YmI7B7o4IiKlkgKwnOTXXzLBWpp2qUutaS+iCZxFREpGgQHYGPOGMWa3MebnbPtqGWPmG2M2eL7XLNliij/89BNMP2cCmzpeBykpgS6OiEipVpgW8HSgb459DwELrLVnAgs8aQlxbdpA7ysjiG5Zwa2eICIiJabAXtDW2sXGmKY5dg8Aeni23wS+Bcb5smDiX6mp7m7zpXPGgLVaukhEpIQV9RlwPWvtLgDP97p5HWiMGWGMiTPGxCUlJRXxdFKS3n/PMqvePez68Du3Q8FXRKTElXgnLGvtVGttrLU2NioqqqRPJ0XQpOZBeqXNJer37wJdFBGRMqOoE3H8aYxpYK3dZYxpAOz2ZaHEPxI2Z9KkqeG83tVhyw9Qu3agiyQiUmYUtQU8B7jFs30L8KlviiO+lpnp3V6xAv74w20vmZ/C/had2XjdQ+6Zb1SUGwAsIiJ+UZhhSDOBH4BWxphEY8ww4DmglzFmA9DLk5YgkpoK93f7nrmXTjixb0XXv7K3+wAAzu8RycFufanXNyZAJRQRKdsK0wt6UB4v9fRxWcSHIiLgin0zOG/358D9APQc0ZyoQ66DVfny0P27ZwNYQhGRsk3LEZYihw7BswPXMPL/atO4WzTs3QuRkVCpUqCLJiJSJmk5wjJi0/pURn09gPKJZ8GaeVCrVqCLJCIieVAALkViukSQumAmEU0bBrooIiJSAHV7LQVmf5zJRw8sw1qIuLgrNGkS6CKJiEgBFIBLgb3PTuHaCV3J+GFFoIsiIiKFpFvQpcCt3w7l8BsVqdL13EAXRURECkkt4BC2dMFx9u1Ow1SuRJW7b9UcziIiIUQBOEQdOwbr+o8jufUFWrtXRCQE6RZ0iKpYEXo91Z0am6u4sb4iIhJSFIBD0PZtlkaNDc3vvwa4JtDFERGRItAt6BCzbFYC+5rEsPLhjwNdFBERKQYF4BDTvnsNKjapS9u+jQNdFBERKQYF4CC0dUMqM1/eze7EVAB+/jyBTxuO4vCe41RqWIMzE+ZT+eJcpxYVEZEQoQAcBJYtg1tjfyJ+4V4Akj5ZyqAx9Uj+eDEAx9Zu4LI/Z5A0Pz6ApRQREV9SAA4CrStv57+rYqg2/RUA2t0cw54nJ9PiilYAnDv2Qirv20GzgecFspgiIuJDCsABsnQp3HILpKdDjXMaEf7uDJr/8x4AIhvWovZjd1GhZSN3cMWKUK1aAEsrIiK+pgAcIGtXp3PTx9dycKFnfeRBg7R8oIhIGaIAHCB3XvMHl9VbS61d6wJdFBERCQAFYD/avx969oS4OCA6mrCf1rr70CIiUuYoAPvRH39A6pYd1Pnvs+7hb8WKgS6SiIgEiAKwH7VuDd+OmkXTt8fDtm2BLo6IiARQsQKwMWasMWadMeZnY8xMY4xWBcjF3Lnw9NOQmQnhD9wHP/0EzZsHulgiIhJARQ7AxpgzgHuAWGttOyAcGOirgpUmX3wBa99eQ+pGT6tXwVdEpMwr7mpI5YCKxpg0oBKws/hFKn0mv5JJZochhA8q53pgGRPoIomISIAVuQVsrd0BTAC2AbuAA9baeTmPM8aMMMbEGWPikpKSil7SEJOZCePHw65dYMLDCP/4A/jvfxV8RUQEKN4t6JrAAKAZ0BCobIy5Oedx1tqp1tpYa21sVFRU0UsaYn77DSY/c5D4cTPdjrPOgo4dA1soEREJGsXphHUZsMVam2StTQM+Brr5pliha80ayMiANm1g3UNvc/nMIbBpU6CLJSIiQaY4AXgbcL4xppIxxgA9gfW+KVZoSNiYzhdN7mTd+E8AmD/7CAdjLuKXB6cDUHvsENcDq0WLAJZSRESCUXGeAS8HPgRWAz958prqo3IFrV0Jx9n27lIAKlUrR8ddcwnb9DsA3ftWomnrijRrFeEOrloVevcOVFFFRCSIGWut304WGxtr4+Li/HY+X7MWPq0zjPMOzaf+0S2YcuGBLpKIiAQxY8wqa21sbq8VdxhSmWIMNPn3Q6Rv6angKyIixaIAXAir4zJJnvEVvSdeTsfrzwTODHSRREQkxGku6EJYPnYWvV++kmNz5ge6KCIiUkqoBVwIt8+/gaRp5Ym6qlegiyIiIqWEWsD5mPPWfg4mHqR8ZDhRo/6iWaxERMRnFIDzsH07pAwdSUqHLpCaGujiiIhIKaNb0Hlo1Aj2T7mH6mm/QEREoIsjIiKljAJwDhkZ8NvqI7Q9tzLnjOiGZtcUEZGSoFvQOfzrwQSqdzmLHc+9HeiiiIhIKaYAnMPNDzQg5dzunHH1uYEuioiIlGJlJgCnp3u3f11veXliJnv2uPQ3r2/li+ajydyzj5r1K9BixUxo3TowBRURkTKhTATguG/2806Nu9j23LsA/Pj1bu6+txx26msApCXtp2fC6+z7clkgiykiImVImeiE1eKcSkRnziXFXgjADYMjSEn8P2pe2hGAPg92wNyxi8iaNQJYShERKUtKbQBO2pHK8sGT6fvF3dSsFwH7fzsxnCisdk0qTfj7iWPDwgAFXxER8aNSewt6zbNz6bfoPra/9pXbobG8IiISREptC/iyyVfz5zU/0+zSswNdFBERkVOUuhbw4me/Y+OcXwCop+ArIiJBqlS1gNNSLXWeuIvMihGwb7kWTxARkaBVqgJw+QhDneVzsfsPKPiKiEhQKzW3oNe8HofNtNSNaUi9Hm0CXRwREZF8lYoAvPpfyzjn9i4sG/56oIsiIiJSKMUKwMaYGsaYD40xvxpj1htjuvqqYKcjZkQXvr/xVTq/dFMgTi8iInLaitsCfhn4ylrbGugArC9+kQpv5ZNz2bl8O2HlwrjwnVFEVK/oz9OLiIgUWZEDsDGmGtAdeB3AWptqrd3vo3IV6MDW/Zz55E1suvkxf51SRETEZ4rTAm4OJAHTjDE/GmP+a4ypnPMgY8wIY0ycMSYuKSmpGKc7WfUmNdg5bR6df3jVZ3mKiIj4S3ECcDmgE/Bva21H4AjwUM6DrLVTrbWx1trYqKioYpzuVG1vOZdKdSr5NE8RERF/KE4ATgQSrbXLPekPcQFZREREClDkAGyt/QPYboxp5dnVE/jFJ6USEREp5Yo7E9bdwDvGmAhgM3Br8YskIiJS+hUrAFtr44FY3xRFRESk7CgVM2GJiIiEGgVgERGRAFAAFhERCQAFYBERkQAw1lr/ncyYJGCrD7OsAyT7ML9AUl2CT2mpB6guwaq01KW01AN8X5cm1tpcZ6HyawD2NWNMnLW2VPTCVl2CT2mpB6guwaq01KW01AP8WxfdghYREQkABWAREZEACPUAPDXQBfAh1SX4lJZ6gOoSrEpLXUpLPcCPdQnpZ8AiIiKhKtRbwCIiIiFJAVhERCQAAhaAjTGNjDGLjDHrjTHrjDF/9eyvZYyZb4zZ4PleM9t7HjbGbDTG/GaM6ZNtf2djzE+e114xxpg8zpnr+4OlLsaYSsaYL4wxv3ryeS6P8zU1xhwzxsR7vqYEW108+7/17MsqZ908zhnsn0vVbHWIN8YkG2P+mcv5SuRzOd16GGNqe44/bIyZnCOvkLpW8qpLKF4rBXwuIXWt5PO5hNq10ssYs8pzTawyxlyaLa+Sv1astQH5AhoAnTzbVYHfgbbA88BDnv0PAf/wbLcF1gAVgGbAJiDc89oKoCtggC+By3M5X57vD5a6AJWASzzHRABL8qhLU+DnEPhcvgViCzhf0H8uueS7Cujur8+lCPWoDFwIjAQm58gr1K6VXOsSotdKfp9LqF0redYlxK6VjkBDz3Y7YIc/r5WAtYCttbustas924eA9cAZwADgTc9hbwJXe7YHALOstcettVuAjUAXY0wDoJq19gfrfiJvZXtPdrm+P5jqYq09aq1d5MknFVgNRPuijIXlq7qcximD/nPJnqcx5kygLu4Pvl+cbj2stUestUuBlOz5hOK1klddQvFayasupyHoP5fsQuRa+dFau9Ozfx0QaYyp4K9rJSieARtjmuL+E1kO1LPW7gL3w8R9gOB+iNuzvS3Rs+8Mz3bO/Tnl9X6fKmZdsudTA+gPLMjjVM2MMT8aY/5njLnIZxU4uQxNKX5dpnluMf0tj1s4IfW5AIOA9zwXZW5K9HMpZD3yEorXSmHyqUFoXCsFCaVrpTBC7Vq5FvjRWnscP10r5Qp7YEkxxlQBPgLGWGsP5nGbHdxtgJxsPvsL+36f8UFdsvIpB8wEXrHWbs7l2F1AY2vtHmNMZ2C2MeZsa+3B4tUgWwF9U5ebrLU7jDFVPXkNxv0nWdj3+4SvPhePgbh65KZEP5fTqEeeWeSyL9ivlYLyCaVrJT+hdq0URshcK8aYs4F/AL2zduVymM+vlYC2gI0x5XE/pHestR97dv/paf5n3TLb7dmfCDTK9vZoYKdnf3Qu+3PK6/0+4aO6ZJkKbLDW/jO3c3lud+zxbK/CPXc4y0dV8VldrLU7PN8PAe+S+62ZkPlcjDEdgHKen/kpSvJzOc165CUUr5WChNK1kqcQvFYKyitkrhVjTDTwCTDEWrvJs9sv10oge0Eb4HVgvbX2pWwvzQFu8WzfAnyabf9Az/35ZsCZwArP7YRDxpjzPXkOyfae7HJ9fzDVxZPXeKA6MCaf80UZY8I9280978/tv/+A1cUYU84YU8eTZ3mgH/BzLqcMic/FYxCutZXX+UrkcylCPXIVotdKfnmF2rWSVz6heK0UJCSuFeMeX3wBPGyt/S7rYL9dK7YEegcW5gvXg84Ca4F4z9cVQG3cs5wNnu+1sr3n/3D/Kf1Gth5pQCzuF3YTMBlOzPB1FfBUQe8Plrrg/nuyuI4DWfncnrMuuGcV63C971YD/YOwLpVxPSDXesr6Mt7e0SH1uWR7bTPQOse+Ev9ciliPBGAvcBj3X3rbEL5WTqkLoXut5FaXUL1Wcv0dC6VrBXgUOJLt2Higrr+uFU1FKSIiEgBB0QtaRESkrFEAFhERCQAFYBERkQBQABYREQkABWAREZEAUAAWEREJAAVgERGRAPj/0TAFT6pygL4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig3 = plt.figure(figsize=(8,4))\n",
    "plt.plot(sol.t,sol.y[0],'b:',label='ivpRK45')\n",
    "plt.plot(sol2.t,sol2.y[0],'r:',label='ivpDOP853')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HA1HMapeQM5E"
   },
   "source": [
    "Note that the native RK45 routine in this package finds a slightly different solution to the 8th order accurate *DOP853* method. If you opt to use intrinsic functions, you need to take care that they are set up correctly, preferably testing them on some simple examples first!\n",
    "\n",
    "Investigate the documentation of the *solve_ivp* function from the *scipy.integrate* package, to explore this process in detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bxqFs-xhb76b"
   },
   "source": [
    "## Further work\n",
    "Now its over to you:\n",
    "\n",
    "Try choosing a worked example seen in the lectures, or indeed a tutorial question **once you have worked through the mathematical steps** and see if you can create a new implementation of RK4 to verify the solution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5UolyYCMV1uQ"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "nlr1Ru5KqTRX",
    "5EQeGB9yFXyS"
   ],
   "name": "MAT301_DifferentialEquationSolving2.ipynb",
   "provenance": [
    {
     "file_id": "https://github.com/john-s-butler-dit/Numerical-Analysis-Python/blob/master/Chapter%2003%20-%20Runge%20Kutta/301_2nd%20Order%20Runge%20Kutta%20Population%20Equations.ipynb",
     "timestamp": 1646748928131
    }
   ]
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}