EENP1_tasks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f6666086b79f74a7",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# <center>Introductory numerical physics & Numerically solving the equations for radioactive decay and for oscillatory motion\n",
    "\n",
    "\n",
    " <center>Department of Physics, University of Surrey module: Energy, Entropy and Numerical Physics (PHY2063)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents:\n",
    "\n",
    "* [Numerical Physics part of Energy, Entropy and Numerical Physics](#EENP-intro)\n",
    "* [Background to numerically solving first and second order ODEs](#ode-background)\n",
    "* [Solving a single first-order ODE on a computer: Euler Method](#euler-method)\n",
    "* [Task 1](#Task-1)\n",
    "* [Solving ODEs using a Python ODE solver](#ode-solver)\n",
    "* [Task 2](#Task-2)\n",
    "* [$F = ma$ as an example of a second-order differential equation](#2nd-order-ode)\n",
    "* [Task 3](#Task-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-fcd26299e06875e0",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "#### Learning Objectives:\n",
    "> Revision of technique to solve first and second order ODEs, also revision of taylor series\n",
    "\n",
    "> Introduce Euler method as a means of computationally solving these equations\n",
    "\n",
    "> Use Euler method to accurately demonstrate the decay of a radioactive sample\n",
    "\n",
    "> Introduce concept of utilising python packages for prebuilt solutions\n",
    "\n",
    "> Use scipy's prebuild ODEint function to solve equations more efficiently\n",
    "\n",
    "\n",
    "## 1) Introduction to the Numerical Physics part of Energy, Entropy and Numerical Physics <a class=\"anchor\" id=\"EENP-intro\"></a>\n",
    "\n",
    "This numerical physics course is part of the second-year\n",
    "Energy, Entropy and Numerical Physics module (PHY2063),\n",
    "and is online at the EENP module on SurreyLearn.\n",
    "See there for assignments, deadlines etc.\n",
    "\n",
    "Very few models in physics can be solved analytically with just pen and paper.\n",
    "Therefore most predictions\n",
    "and modelling done in physics requires numerically solving the equations of the model.\n",
    "This is true for electromagnetism, quantum physics (including nuclear and particle physics),\n",
    "astrophysics, soft matter physics, etc.\n",
    "This course introduces numerical physics through\n",
    "some simple common examples, and discusses general\n",
    "points about numerical physics, in particular the key problem of estimating\n",
    "how accurate the numerical solution is.\n",
    "\n",
    "It is taught using Python, which is a good language for numerical physics;\n",
    "One of the main strengths of Python is its modularity and so it often makes it trivial to solve \n",
    "otherwise complex problems by simply importing a module that has a premade solution. Additionally\n",
    "it is widely used throughout scientific communities meaning there are many very useful and efficient\n",
    "scientific modules that you will become familiar with during this course. Additionally Python is\n",
    "widely used outside of academia and so it is a very useful skill to have regardless of what career\n",
    "path you intend to follow. However almost all high performance computing (e.g., solving large PDEs,\n",
    "quantum matrix problems, molecular dynamics, fluid mechanics, etc) is done using Fortran or C. So Python\n",
    "is by no means the best at every aspect of computing, however for the size of the problems you will be \n",
    "handling python will be more than sufficient. The numerical physics skills you will learn using Python \n",
    "will be transferable to number crunching applications in other computer languages. And the coding \n",
    "skills are transferable to a huge variety of other tasks people use computers to do.\n",
    "\n",
    "The course looks at numerically solving ODEs (ordinary differential equations) and\n",
    "PDEs (partial differential equations), and introduces the numerical modelling\n",
    "of random processes, using both \n",
    "a very widely used numerical\n",
    "technique called Monte Carlo (MC), and the use of Bayesian statistics.\n",
    "MC is used a lot in fields from statistical physics,\n",
    "to nuclear and particle physics. \n",
    "Bayesian statistics is used from particle physics to predictive text on mobile phones.\n",
    "A lot of modelling in physics is based on PDEs,\n",
    "from heat and mass flows in statistical physics, a lot of electromagnetism and\n",
    "quantum physics, etc. And the numerical techniques we will learn\n",
    "are also used extensively outside\n",
    "physics, for example to model share prices in stock markets.\n",
    "Increasing numbers of graduates in physics (and related subjects) are going\n",
    "into data science, which uses many models of random processes such as Monte Carlo modelling\n",
    "and Bayesian statistics.\n",
    "\n",
    "Examples of ODEs are the equations for radioactive decay and for a mass attached\n",
    "to a spring.\n",
    "PDEs are even more common in physics than ODEs, for example, Schrödinger's\n",
    "time-independent equation in three dimensions is a PDE, Maxwell's equations\n",
    "which govern electromagnetism (including light), and the wave equation,\n",
    "are all PDEs. In many cases (all except very simple cases) we cannot\n",
    "solve them analytically, and so can only solve them numerically.\n",
    "\n",
    "### 1.1) Course structure\n",
    "\n",
    "The autumn semester computing is composed of 3 sections, ODEs, PDEs and Data Analysis.\n",
    "Of these data analysis is the largest and made up of several subsections, Probability and Monte Carlo\n",
    "methods, uncertainty estimates, fitting.\n",
    "\n",
    "Each topic will have an assignment, the first (ODEs) is formative. The second assignment is mainly based on the quality of the\n",
    "code that you produce and then the final assignment has you write a small report alongside some code \n",
    "in order to demonstrate greater understanding.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f621c2de542938ad",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## 2) Background to numerically solving first and second order ODEs <a class=\"anchor\" id=\"ode-background\"></a>\n",
    "\n",
    "We will start with a very simple example, the ODE for the mass of a radioactive isotope remaining\n",
    "after a time $t$, $M(t)$.\n",
    "If the rate of\n",
    "decay of the isotope is $r$ then the ODE is\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}M(t)}{{\\rm d}t}=-rM(t)\n",
    "$$\n",
    "\n",
    "and if the decay rate $r=0.3$s$^{-1}$ then\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}M(t)}{{\\rm d}t}=-0.3M(t)\n",
    "$$\n",
    "\n",
    "Now as this a first-order ODE we know to obtain a particular solution we\n",
    "need ${\\bf exactly~one}$ boundary condition (BC).\n",
    "Usually this is the value of the mass of the isotope at time $t=0$.\n",
    "If for example, $M(t=0)=0.1$kg then the particular solution is\n",
    "\n",
    "$$\n",
    "M(t)=M(t=0)\\exp(-rt)=0.1\\exp(-0.3t)\n",
    "$$\n",
    "\n",
    "this is the familiar exponential decay of the amount of remaining\n",
    "isotope.\n",
    "\n",
    "This is just one example,\n",
    "ODEs are common in physics, e.g., the time-independent\n",
    "one-dimensional Schrodinger equation is a second order ODE\n",
    "\n",
    "$$\n",
    "-\\frac{\\hbar^2}{2m}\\frac{{\\rm d}^2\\psi(x)}{{\\rm d}x^2}+\n",
    "V(x)\\psi(x)=E\\psi(x)\n",
    "$$\n",
    "\n",
    "which we need to solve for the wavefunction $\\psi(x)$ and energy $E$.\n",
    "Here $V(x)$ is the potential energy function.\n",
    "Another ODE is the equation for a mass $m$ attached to a spring\n",
    "of spring constant $k$\n",
    "\n",
    "$$\n",
    "m\\frac{{\\rm d}^2x}{{\\rm d}t^2}=-kx\n",
    "$$\n",
    "\n",
    "where we want to solve for the position as a function of time, $x(t)$. \n",
    "\n",
    "\n",
    "\n",
    "### 2.1) Taylor Series Revision\n",
    "To get the formula for the Euler method we will need\n",
    "the Taylor series expansion expression for the value of $M(t)$\n",
    "at the point $t_f=t_s+h$ in terms of $M$ and its derivatives at the nearby point\n",
    "$t=t_s$. Nearby means that the distance $h$ between the two points is small.\n",
    "You did the Taylor series in first year maths.\n",
    "The series expansion for the point $t=t_s$ is\n",
    "\n",
    "<center>$M(t_f)=M(t_s+h)=M(t_s)+h\\left(\\frac{{\\rm d}M}{{\\rm d}t}\\right)_{t_s}\n",
    "+\\frac{1}{2}h^2\\left(\\frac{{\\rm d}^2M}{{\\rm d}t^2}\\right)_{t_s}+ .... $\n",
    "\n",
    "where the subscripts $t_s$ indicate that the derivatives are evaluated at\n",
    "the point $t=t_s$.\n",
    "\n",
    "We can also write this as\n",
    "\n",
    "<center>$M_f=M_s+hM'_s+\\frac{1}{2}h^2M_s''+ ....$\n",
    "\n",
    "where I use the standard notation of a dash to indicate the first derivative, and\n",
    "two dashes to indicate the second derivative. Here $M_f=M(t_f)$, $M_s=M(t_s)$, $M_s'=M'(t_s)$\n",
    "and $M_s''=M''(t_s)$.\n",
    "If $h$ is small we can drop all but the first two terms\n",
    "\n",
    "<center>$M_f\\simeq M_s+hM'_s$\n",
    "\n",
    "Note that $M(t)$ denotes a function $M$ of the variable $t$ but $M_s$ and $M'_s$\n",
    "are not functions they are just numbers, the values of the functions $M$ and $M'$\n",
    "at the point $t_s$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c68e3fe9d7f8bd1a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## 3) Solving a single first-order ODE on a computer: Euler Method <a class=\"anchor\" id=\"euler-method\"></a>\n",
    "\n",
    "We will introduce the simple Euler method\n",
    "for numerically integrating an ODE, and then also introduce\n",
    "the more accurate and slightly more complex modified Euler method.\n",
    "For all methods of solution it is important to remember that\n",
    "for a first-order ODE\n",
    "the solution of a\n",
    "first order ordinary differential equation requires a single boundary\n",
    "condition to be specified (this is usually an initial condition). The\n",
    "specification of the value of the required function at a single point\n",
    " is thus\n",
    "necessary to generate the required solution.\n",
    "\n",
    "So we would like to solve the\n",
    "first-order ODE\n",
    "\n",
    "<center>$\\frac{{\\rm d}M(t)}{{\\rm d}t}=-rM(t)$ &emsp;&emsp;&emsp; $r=0.3~{\\rm s}^{-1}$</center><br>\n",
    "\n",
    "with the BC $M(t=0)=0.1$, and\n",
    "over some range of values of $t$, say from $t=0$\n",
    "to $t=30$~s.\n",
    "\n",
    "The above ODE is just an example of a first order ODE, in general\n",
    "the derivative is some function of the variable, here $t$.\n",
    "The key point to take from the above equation is that\n",
    "at any point along the $t$ axis the equation does not tell us\n",
    "what we ultimately want to know, the value of the function\n",
    "$M(t)$, but it does tell us the slope of the function, i.e., how fast\n",
    "it is changing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-6550812111ee4c70",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# useful imports for arrays and plotting capabilities\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# constants and initial boundary conditions\n",
    "r_decay = 0.3\n",
    "t_start = 0\n",
    "t_final = 30\n",
    "mass_initial = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ac604e429cf0db08",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 3.1) The Euler Method\n",
    "\n",
    "The Euler method  works by dividing the range between $t=0$ and $t=30$~s\n",
    "into many small\n",
    "intervals each of length $h$: this\n",
    "is called the step length. For example, if $h=0.1$, then\n",
    "Euler's method of solving an ODE starts from $t=0$ and\n",
    "integrates from there to $t=30$~s, in $30/0.1=300$ little steps, each $0.1$ in size..\n",
    "The method works one step of length $h$ at a time.\n",
    "\n",
    "We start at the point $t=0$. At that point we use the known value of\n",
    "$M(t=0)$ to calculate\n",
    "the value of the function a distance $h$ away along the time axis,\n",
    "i.e., we calculate $M(t=h)$. Then we use the now known value\n",
    "of $M(t=h)$ to calculate the slope at $t=h$, and so calculate\n",
    "$M(t=2h)$, and by repeating this we integrate along the $t$ axis.\n",
    "\n",
    "\n",
    "#### 3.1.1) Single Euler Step\n",
    "Euler's algorithm for a single step starts with a point along the $t$-axis,\n",
    "call it $t_s$ where we know the value of the function $M(t_s)=M_s$\n",
    "and its derivative, $M'(t_s)=M'_s$. Note that as $M'$ is a known function\n",
    "of $t$ and $M$ (this is the ODE) whenever we know $t$ and $M$ we can calculate $M'$.\n",
    "We need the expression we got from truncating the Taylor series in the box above.\n",
    "This is\n",
    "\n",
    "<center>$M_f\\simeq M_s+hM'_s$<br></center>\n",
    "\n",
    "this tells us how to calculate the value of the function $M(t)$ at the end\n",
    "of the step, $M_f$, in terms of the value of the function at the start\n",
    "of the step, $M_s$ {\\em and} its derivative at the start of the step, $M'_s$.\n",
    "You can start writing a program to solve an ODE, by writing one to do a single\n",
    "step like this.\n",
    "\n",
    "In terms of what we need to calculate in a single step, it looks like\n",
    "\n",
    "<center>$M_s'=-rM_s$</center><br>\n",
    "<center>$t_f =t_s+h$</center><br>\n",
    "<center>$M_f =M_s+hM_s'$</center><br>\n",
    "\n",
    "where I have picked sensible variable names. You should always use sensible\n",
    "variable names as this makes debugging much easier. With poor variable names,\n",
    "if you leave a program for a few days, it can be difficult for you to remember how it works.\n",
    "These lines implement one step, and at the end of this step, $t_f$ \n",
    "(${\\tt t\\_final}$) and\n",
    "$M_f$ (${\\tt mass\\_final}$) are the values of $t$ a distance $h$ along the $t$-axis, and\n",
    "our estimate for the $M$ value there, respectively.\n",
    "\n",
    "#### 3.1.2) Obtaining the complete function (array)\n",
    "\n",
    "To obtain the complete function $M(t)$ from $t=0$ to $t=30$~s,\n",
    "at the discrete set of points $t=0,h,2h,3h, ... , 30$, we need\n",
    "to appreciate that once we discretise the time axis into\n",
    "a grid of spacing $h$, the function $M$ becomes an array\n",
    "\n",
    "${M(0), M(h), M(2h), ..., M(30.0)} → \\{ {\\tt mass(0)}, {\\tt mass(1)}, {\\tt mass(2)}, ... , {\\tt mass(300)}\\}$\n",
    "\n",
    "where on the left hand side we have the function of time $t$, and\n",
    "on the right hand side we have the array with its elements denoted by\n",
    "$i$, which is an integer.\n",
    "This is\n",
    "if $h=0.1$~s, and so we need array elements 0 to 300 to have enough\n",
    "elements to go from $t=0$ to $t=30$s in steps of 0.1.\n",
    "\n",
    "Then both the time and mass arrays can be calculated in increments of one step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-331ec09ce84a3821",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# step size in seconds\n",
    "h = 0.1\n",
    "\n",
    "# calculating values 1 step forwards in time\n",
    "massdash = -r_decay*mass_initial\n",
    "t_final = t_start + h\n",
    "mass_final = mass_initial + h*massdash"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d59144b747c7528a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "in the first line of the loop we calculate the derivative (which we don't need\n",
    "to be an array), then in the next two lines we calculate\n",
    "the next values of the arrays ${\\tt time}$ and ${\\tt mass}$.\n",
    "\n",
    "### Task 1 <a class=\"anchor\" id=\"Task-1\"></a>\n",
    "\n",
    "Write a Python program to evaluate the solution of the first-order ordinary\n",
    "differential equation for radioactive decay at a rate $r=0.3$s$^{-1}$\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}M(t)}{{\\rm d}t}=-0.3M(t)\n",
    "$$\n",
    "\n",
    "with the single BC $M(t=0)=0.1$kg. Obtain the mass as a function of time\n",
    "from $t=0$s to $t=30$s, in steps of $0.1$s (300 steps).\n",
    "Note that here we know the solution, it is\n",
    "\n",
    "$$\n",
    "M(t)=M(t=0)\\exp(-rt)=0.1\\exp(-0.3t)\n",
    "$$\n",
    "\n",
    "When you write any program, it is always best to at first get it to calculate\n",
    "something you know already, so you can compare solution\n",
    "to what you know it should be to check and debug\n",
    "the program. Then later you and modify the program you know works to do something new.\n",
    "\n",
    "When you have the exact solution this is not only helpful to debug the program,\n",
    "you can use to see how accurate the solution is, for example as you vary $h$.\n",
    "As $h\\to0$ then the solution\n",
    "is exact but then you calculate an infinite number of points,\n",
    "which is not possible. So for $h>0$ the array you calculate\n",
    "is approximate. Basically all numerical physics is approximate,\n",
    "at least to some degree.\n",
    "\n",
    "Check your program by plotting the values it gives you and the true values on the same graph. If it works the lines should be basically on top of one another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c36c1e40153c87e9",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0010369555489464187\n",
      "0.0011108996538242435\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8FdX9//HX595sBELYF1lkFwVX4oq7VbFWrVjXttrW79f2W21r7fKtfmtr7ab218W2amvVVq2K1pVWK651AUUCoggIshMEAdkhIcv9/P6Yyc0lZLlgbiY3eT8fj3lkzpkzdz7jlXwyc2bOMXdHRESkKbGoAxARkbZPyUJERJqlZCEiIs1SshARkWYpWYiISLOULEREpFlKFiIi0iwlCxERaZaShYiINCsn6gBaSq9evXzIkCFRhyEiklVmzpy53t17N9eu3SSLIUOGUFpaGnUYIiJZxcyWp9NOt6FERKRZShYiItIsJQsREWmWkoWIiDRLyUJERJqV0WRhZhPMbIGZLTKzHzSw/Xgzm2Vm1Wb2uXrbLjOzD8LlskzGKSIiTcvYo7NmFgduA04FyoAZZjbZ3eelNFsBfAn4br19ewA/BkoAB2aG+25syRjdnXVbd/L+mq1s3FHJOYcMaMmPFxFpNzL5nsURwCJ3XwJgZpOAc4BksnD3ZeG2RL19Tweed/cN4fbngQnAQy0Z4OrNFRxz00sAFHfK5exR3bHCwpY8hIhIu5DJ21ADgJUp5bKwrsX2NbMrzKzUzErXrVu3xwH2X76QIq8CYHN5FWu++b09/gwRkY4gk8nCGqjzltzX3e909xJ3L+ndu9m31Xc/yNq17F+2IFl+v6xF73KJiLQbmUwWZcCglPJA4MNW2Dd9Bx7IfuuWJYvztzvU1LT4YUREsl0mk8UMYKSZDTWzPOAiYHKa+04BTjOz7mbWHTgtrGtZffsyesf6ZPH97gNh8eIWP4yISLbLWLJw92rgKoJf8vOBR9x9rpndaGZnA5jZ4WZWBpwP/NnM5ob7bgB+SpBwZgA31nZ2tygzRvfITxYX9B4Cc+a0+GFERLJdRkeddfdngGfq1f0oZX0GwS2mhva9B7gnk/EB7Desb3J9cY+B7Hx3NvnnnZfpw4qIZJUO/wZ3lwMPYPDG1QBUx3NYvLAs4ohERNqeDp8sOPBARqd0cr+/dnt0sYiItFFKFmPGMHp93dwf73sn2LEjwoBERNoeJYvOnRltdclhfu8hMG9e4+1FRDogJQtgdN8uyfUFvfaFd9+NMBoRkbZHyQLYd/QQCqoqAFhb1JOP57wfcUQiIm2LkgUQP+hA9lu3IllesGzPx5kSEWnPlCwgGPZj/bJkcf6mquhiERFpg5QsAEaMYPTGuvcr3j/sOKhSwhARqaVkARCPM/qPtySL7w8ZA7m5EQYkItK2KFmERo+uG+R24UdbqUmkO5q6iEj7p2QR6tE5j75dg0EFd1YnWLp+W8QRiYi0HUoWKcbuU5xcn7Nqc4SRiIi0LUoWKcYMqEsW763aEmEkIiJti5JFirH7dE2uv7f4I9jQ8lNoiIhkIyWLFGPfeD65Pm/ZOhKPPhZhNCIibYeSRYr+Hy6jx46gr2JrfmdWvLMg4ohERNoGJYsUNu4wxnxUNwf3eyt0G0pEBJQsdnXooYxdk5IsyuNQXR1hQCIibYOSRar+/RlbsT5ZnNtzMCzQrSgRESWLVGaM7dc5WXyv73B81tsRBiQi0jYoWdQzeMwIiiqCt7c3Fhazavb8iCMSEYmekkU9dtihjE3t5F66NsJoRETaBiWL+g7dNVnM3ergGlRQRDo2JYv6hg1j7JbVyeJ73QbC0qURBiQiEj0li/rMGNOrIFl8r+9weFud3CLSsSlZNGDoAUMorCwHYF2XHnw0c07EEYmIREvJogHxww7bpd9i9qKPIoxGRCR6ShYNOfxwDl69MFmcvTkBiUSEAYmIREvJoiEjR3JIed0js7OP+7SShYh0aBlNFmY2wcwWmNkiM/tBA9vzzezhcPt0MxsS1uea2b1mNsfM5pvZtZmMczexGIc8PSlZfLcyn5pYvFVDEBFpSzKWLMwsDtwGnAEcAFxsZgfUa3Y5sNHdRwC/BW4O688H8t39QGAc8NXaRNJa9uleSO+iYE7u7ZU1LF6nOblFpOPK5JXFEcAid1/i7pXAJOCcem3OAe4N1x8FTjEzAxzobGY5QCegEmjVeU7NjEMGdUuWZ6/Y1JqHFxFpUzKZLAYAK1PKZWFdg23cvRrYDPQkSBzbgdXACuD/uXurTy6RmizeXqlkISIdVyaThTVQV3/cjMbaHAHUAPsAQ4HvmNmw3Q5gdoWZlZpZ6bp16z5pvLvZ5cpCyUJEOrBMJosyYFBKeSDwYWNtwltOxcAG4BLgWXevcve1wFSgpP4B3P1Ody9x95LevXu3+AkcNLAYC9PZwtWb2fHq1BY/hohINshkspgBjDSzoWaWB1wETK7XZjJwWbj+OeAld3eCW08nW6AzcBTwfgZjbVDR05MZsWEVADUY7933eGuHICLSJmQsWYR9EFcBU4D5wCPuPtfMbjSzs8NmdwM9zWwRcA1Q+3jtbUAX4D2CpPNXd383U7E2qmdPDlk5L1mcvWZ7q4cgItIW5GTyw939GeCZenU/SlmvIHhMtv5+2xqqb3XjxnHI6lv4x0GnAjA7pzts3gzFxREHJiLSuvQGd1OKijiksDpZfKf/KCgtjTAgEZFoKFk0Y7/RgymoqgBgVXEf1r45K+KIRERan5JFM3KOOJyD1ixKlmfOXxVhNCIi0VCyaM5RR1FSVtfJPWNjQtOsikiHo2TRnLFjOXzDsmSxtNdQWLAgunhERCLQ7NNQZlYCHEfwNnU5weOsL0Qx/EYk4nEOG9wN8wRuMeb2Hc7216bRefToqCMTEWk1jV5ZmNmXzGwWcC3BYH4LgLXAscDz4RDig1snzGgVH1XCfuuWA1ATizO7VFcWItKxNHVl0RkY7+7lDW00s0OAkQRvW7dvxx5LyWv/4P0+QwGYsXoH4yMOSUSkNTV6ZeHutzWWKMLts939xcyE1cYccQSHf1g32khpYV9Yu7aJHURE2pd0+iz+wO6jxW4GSt39qYxE1dZ07kxJt7q8Omuf0VS/PpWciedGGJSISOtJ52mofOAQ4INwOQjoAVxuZr/LYGxtyoDDD2SfLcHVxI68TswvnR9xRCIirSedZDECONnd/+DufwA+BewPnAuclsng2pRPf5qSvJ3J4oxPXxRhMCIirSudZDGAoLO7VmdgH3evAXY2vEs7dOqpHH7hhGSxdPnGCIMREWld6Yw6ewsw28z+QzCz3fHAL8J5Jl7IYGxtTsmQHsn1Gcs24u6YNTTZn4hI+9JssnD3u83sGYKpTg24zt1rZ7z7XiaDa2tG9S2iqCCHrRXVrNu6k+Uf72BIr87N7ygikuWavQ1lZje6+2p3f8rdnwTWmNkDrRBbmxOPGSX7dk+W31zycYTRiIi0nnT6LAab2bUAZpYPPEnwVFSHdMzwXsn1aYuVLESkY0gnWXwZODBMGP8EXnb3GzIaVRt29IC6207TZizE33wzwmhERFpHU2NDHWZmhwGHArcCFxJcUbwS1ndI+19/DcXlWwFYn9eFxc9PjTgiEZHMa6qD+9f1yhuBA8J6B07OVFBtWfyE4zlyyns8N+poAKbN/5AREcckIpJpjSYLdz+pNQPJGiedxDF3fieZLN6oLOTSnTshPz/iwEREMqep21BfsCZeIjCz4WZ2bGbCasNGjuSY8tXJ4hsDDiAxfXqEAYmIZF5Tt6F6EryMNxOYCawDCgiG/zgBWA/8IOMRtjVmjDxsNL22b2R95+5s6tSV+S9NZ8zxx0cdmYhIxjQ1RPmtwGHAQ0Bv4JSwvAr4oruf5+4d8hFaO+kkjloxJ1l+Y8FHEUYjIpJ5Tb7BHY7/9Hy4SK2TTuKYP/yDf+0fXE28UVPEf1VUQEFBxIGJiGRGOu9ZSH1Dh3J01fpkcfqAA6ieOi3CgEREMkvJYm+YMeTwsfTfsg6AbfmFvPtyacRBiYhkjpLFXrKTTuKY5e8ky69p6A8RacfSGUiwr5ndbWb/DssHmNnlmQ+tjTvlFE5YOitZfCWnD2zaFGFAIiKZk86Vxd+AKcA+YXkhcHWmAsoaAwdyXN4OzBMAzO4/ks1TXow4KBGRzEgnWfRy90eABIC7VwM1GY0qS3Q/5XgOWh08PZyIxXn91XcjjkhEJDPSSRbbzawnwXhQmNlRwOZ0PtzMJpjZAjNbZGa7vcBnZvlm9nC4fbqZDUnZdpCZvWFmc81sjpm1vedSTz+dE5bOTBZfGXhQhMGIiGROOsniGmAyMNzMpgL3Ad9obicziwO3AWcQDEB4sZkdUK/Z5cBGdx8B/Ba4Odw3B/g78DV3HwOcCFSlc0Kt6vjjOeGbX0wWX/Fi3D3CgEREMqPZZOHuswiG9zgG+Cowxt3Tud9yBLDI3Ze4eyUwCTinXptzgHvD9UeBU8LxqE4D3nX3d8IYPg5fEGxbCgo4eOJpFHfKBeCjLTtZ8NHWiIMSEWl5jb7BbWYTG9k0ysxw98eb+ewBwMqUchlwZGNt3L3azDYTjEk1CnAzm0Iw1Mgkd7+lgRivAK4AGDx4cDPhZEZOPMaxI3vx9LvB4IKvLFjH6H5dI4lFRCRTmhru46wmtjnQXLJoaMTa+vdoGmuTAxwLHA7sAF40s5nuvsvjRu5+J3AnQElJSWT3f04Y1bsuWSxcx1dPGB5VKCIiGdHUfBZf/oSfXQYMSikPBD5spE1Z2E9RDGwI619x9/UAZvYMwSCGbfLZ1ONH9k6uly7byPad1XTOb3LYLRGRrJLOS3k9zez3ZjbLzGaa2a3h01HNmQGMNLOhZpYHXETQUZ5qMnBZuP454CUPeoinAAeZWWGYRE4A5qV7Uq2tX3EBo/sVAVBZk2DqY20yp4mI7LV0noaaRDCXxXkEv9DXAQ83t1P4PsZVBL/45wOPuPtcM7vRzM4Om90N9DSzRQRPXf0g3Hcj8BuChDMbmOXuT+/JibWqlSs5ceGbyeKLkzUvt4i0L9bco55hX8G4enWl7l6S0cj2UElJiZeWRjSY3+bNzDxwPOddcjMAvbZvZPqVhxMfvV808YiIpCn8Hd/s7/N0rixeNrOLzCwWLhcAbfev/CgUF3PIfvvQa/tGANZ37s7sJ1+IOCgRkZaTTrL4KvAgUBkuk4BrzGyrmW3JZHDZJH72WZyy6K1k+fm5ayKMRkSkZaXzUl6Ru8fcPSdcYmFdkbvrhYJaZ53FqR9MTxafLxgAH2vYchFpH9KazyIcp+lsM5tYu2Q6sKwzZAjji2ooqKoAYHHPQSx54tmIgxIRaRnpPDp7D3APwdNQZ4XLZzIcV1bq9JkzOG7p28nyC29+EGE0IiItJ503x45y9/oDAEpDzj6bU5/6P54fdTQAL1R24YqdOyE/P+LAREQ+mXRuQ73RwGix0pBx4zh52/LkhEil/Uax4bmXIw5KROSTSydZ3EuQMBaY2bvh3BKa5achsRi9TjuJcavmA8GESC+8MKuZnURE2r50ksU9wBeBCdT1VzQ1yGDHdvbZnJ7yNvfTG2KQSEQYkIjIJ5dOsljh7pPdfam7L69dMh5Ztjr5ZM4oq+vkntpvfza+/FqEAYmIfHLpJIv3zexBM7tYj86moaCAgScdw6Gr3gegOp7Dc/96s5mdRETatnSSRSdgJ8HsdXp0Nh3nn8+Z79ddTfwr1ifCYEREPrlmH51tgXktOp4JE/j0I0/xs7A4Lb8vG7ZX0qNzXqRhiYjsrXReyiswsyvN7HYzu6d2aY3gslZBAfs8+FfG7dsdgJqEM0VjRYlIFkvnNtT9QD/gdOAVghnvtmYyqPbizAP7J9drp10VEclG6SSLEe5+PbDd3e8FzgQOzGxY7cOnU5LFtMXr+XjbzgijERHZe+kki6rw5yYzG0swT/aQjEXUjvQrLuDwIcGtqITD03N0dSEi2SmdZHGnmXUHrieYM3secHNGo2pHzjp4n+T64y/NhaqqJlqLiLRN6cxncZe7b3T3V9x9mLv3cfc/t0Zw7cFn4hvJDceKmr0VFj/1XMQRiYjsuUaThZmdZWb7ppR/ZGbvmNlkMxvaOuFlvx5/uZ2TUiZFeuKFORFGIyKyd5q6svg5sA7AzD4DfAH4CsGtqD9lPrR24tJLmfjeS8niE7F+JDZujDAgEZE911SycHffEa5PBO5295nufhfQO/OhtRPjx3NSYj3dyoPpyld17c30+/8ZcVAiInumqWRhZtbFzGLAKcCLKdsKMhtWOxKLkf/5izlr/qvJqsdnr4owIBGRPddUsvgdMBsoBea7eymAmR0K6BnQPfGFL+xyK+qZ4hGUL1wcYUAiInum0WTh7vcAJwCXA59O2bQG0HhRe2LkSA4Z0pNhH5cBsD2/kGcfeDbioERE0tfko7Puvsrd33b3RErdandfkfnQ2he79FLOe6/uTt5DZdWaFElEskY6L+VJS7jgAs5f8Co5NdUAvNV7BB/888VmdhIRaRuaes9C71K0pB496DPhZD61qO6diwefnR1hQCIi6WvqyuJRADPTn78t5YoruGR2XV/FY/mDqViloctFpO1rKlnEzOzHwCgzu6b+0loBtisnnMCx+eUM3hg8TLaloAtP3/1UxEGJiDSvqWRxEVBBMJteUQNLs8xsgpktMLNFZvaDBrbnm9nD4fbpZjak3vbBZrbNzL6b3um0cWbErvhvLnpnSrLqwRWV6ugWkTav0WlV3X0BcLOZvevu/97TDzazOHAbcCpQBswws8nuPi+l2eXARncfYWYXEYxme2HK9t8Ce3zsNu3SSzn/pzfxm+O+QHU8h5m9hrFg6Vr2G94v6shERBqVztNQ08zsN2ZWGi6/NrPiNPY7Aljk7kvcvRKYBJxTr805wL3h+qPAKWZmAGb2WWAJMDetM8kWvXrR+1tf5/Sulcmq+95dF2FAIiLNSydZ3EMwjeoF4bIF+Gsa+w0AVqaUy8K6Btu4ezWwGehpZp2B/wV+ksZxss8NN/CFS05OFh+bVcbG7ZVN7CAiEq10ksVwd/9xeIWwxN1/AgxLYz9roM7TbPMT4Lfuvq3JA5hdUXvFs25ddv11ftSwHhzQvysAFVUJHnxL7zmKSNuVTrIoN7NjawtmNh4oT2O/MmBQSnkg8GFjbcwsh2DK1g3AkcAtZrYMuBq4zsyuqn8Ad7/T3UvcvaR37+waCNfM+Mqxda+y3PfGMqpq1NEtIm1TOsnia8BtZrYs/OX9R+Craew3AxhpZkPNLI/g6arJ9dpMBi4L1z8HvOSB49x9iLsPIRjQ8Bfu/sc0jplVzjq4P7265APw0ZadPKM5ukWkjUpnWtV33P1g4CDgIHc/1N3fTWO/auAqYAowH3jE3eea2Y1mdnbY7G6CPopFwDXAbo/Xtmf5OXEuPTo5GSF33/ciXlMTYUQiIg0z9/rdCNmppKTES0tLow5jj63/zR85ZvVAKuO5APzjYDj84jMjjkpEOgozm+nuJc2100CCEes1603OTZnr4k+vLIkwGhGRhilZRO3b3+a/ZzyBhaPAv9htGPOenxZxUCIiu2o2WZjZ+WZWFK7/0MweN7PDMh9aBzFuHCPGDGPCwjeSVbc9MTPCgEREdpfOlcX17r41fHz2dII3ru/IbFgdzHXXceW0h5PFZ4qGsHjq2xEGJCKyq3SSRe3jOWcCd7j7U0Be5kLqgE47jbGDunPi4qCD3i3G7Q+9FnFQIiJ10kkWq8zszwRDfTxjZvlp7ifpMoMf/pBvvDEpWfVk4RBWvj0/wqBEROqk80v/AoJ3JSa4+yagB/C9jEbVEZ19NuO6xTlqefAKS00szm33vtTMTiIirSOdZNEfeNrdPzCzE4HzgbcyGlVHFIvBddfxzWl1Vxf/yBvEkpm6uhCR6KWTLB4DasxsBMEb10OBBzMaVUd1wQUck7eDY5a9AwRXF7+59+WIgxIRSS9ZJMKhOyYCv3P3bxNcbUhLy8mBG27ge6/em6z6V+G+vPef7HszXUTal3SSRZWZXQxcCvwrrMvNXEgd3MUXc2ivfE5Lee/i//1jRoQBiYiklyy+DBwN/Nzdl5rZUODvmQ2rA4vF4Be/4Luv3Z98q/s/RYOZvvCjiAMTkY4snVFn57n7N939obC81N1vynxoHdiZZzJqv0Gcu/mDZNXPn/uARKJ9DPooItknneE+RprZo2Y2z8yW1C6tEVyHZQbPPce3f/E18nKCr+jdss08/vaqiAMTkY4qndtQfyUY3qMaOAm4D7g/k0EJUFjIoB6FXHFc3Qy2Nz/7Ptt2VkcYlIh0VOkki07u/iLB3BfL3f0G4OTMhiW1/ufE4fQpCmbTW7d1J7e/vCjiiESkI0onWVSYWQz4wMyuMrNzgT4ZjktCnfNz+MEZo5Plu15fyopVH0cYkYh0ROkki6uBQuCbwDjgi9TNmy2t4LOHDOCQQd0AqKxO8JNbHqO9zHAoItkhnaehZrj7Nncvc/cvu/tEd3+zNYKTQGzHdn783lN1EyR1GsC/J70QcVQi0pHkNLbBzCY3taO7n93y4UiDCgo4dOqzXNJ/Ow8c+mkAfjx9PeMnbKW4e1HEwYlIR9BosiB4EW8l8BAwHbBWiUh2l5MDd9zB9085nedHHMnaop6sK+jKLb/6Bz//xVeijk5EOoCmbkP1A64DxgK3AqcC6939FXd/pTWCkxTjx1P85S9yw4t3JqseSPSl9CWNGyUimddosnD3Gnd/1t0vA44CFgH/MbNvtFp0squbbuKMLUs4ZVHdCPHfe3IeO8p3RhiUiHQETXZwm1m+mU0kGAvqSuD3wOOtEZg0oLgY++MfufG5O+iycwcASwt78stfPtzMjiIin0yjycLM7gWmAYcBP3H3w939p+6uMSeidO65DDj5GH6Ucjvq/sqevPJvPaAmIpnT1JXFF4FRwLeAaWa2JVy2mtmW1glPdmMGf/oT53/4NqemDGP+/SlL2LRha4SBiUh71lSfRczdi8Kla8pS5O5dWzNIqadfP+zOO/nllD/Sc/smAD4qKObamx7Vy3oikhHpvMEtbdHEifQ6/7P88tk/JKv+HevDfXc9E2FQItJeKVlks1tv5bSqNXxx1r+SVT9bWM07MxdGGJSItEdKFtmsuBj+/nd++MpfGbsmGI22Kp7Dlc+XsXlHVcTBiUh7omSR7Y49lvyf3sjtL99OUTyoKttUzrcfmU2NZtYTkRaS0WRhZhPMbIGZLTKzHzSwPd/MHg63TzezIWH9qWY208zmhD81f0ZTvvtdBr/5H3518bhk1Uvvr+WWKe9HGJSItCcZSxZmFgduA84ADgAuNrMD6jW7HNjo7iOA3wI3h/XrgbPc/UCC4dA1M19TYjHo358JY/vxtROGJ6v//MoSHptZFmFgItJeZPLK4ghgkbsvcfdKYBJwTr025wD3huuPAqeYmbn72+7+YVg/Fygws/wMxtpufO/0/fjU/nVzU137+BxmLlgdYUQi0h5kMlkMIBi1tlZZWNdgG3evBjYDPeu1OQ942913GwDJzK4ws1IzK123bl2LBZ7N4jHjdxcdyqi+XQCorEnw3396lSVLlDBEZO9lMlk0NKR5/R7XJtuY2RiCW1NfbegA7n6nu5e4e0nv3r33OtD2pkt+Dnedtz/dq8sB2JDfhUt//xJr126KODIRyVaZTBZlwKCU8kDgw8bamFkOUAxsCMsDgSeAS919cQbjbJcG3/Zr7nnw/+hUWQFAWUE3LvvFZLZsK484MhHJRplMFjOAkWY21MzygIuA+rPvTaZuPu/PAS+5u5tZN+Bp4Fp3n5rBGNuvH/+YQ4f15vanfkk8UQPA/IKeXP7jR9heXhlxcCKSbTKWLMI+iKuAKcB84BF3n2tmN5pZ7ZSsdwM9zWwRcA1Q+3jtVcAI4Hozmx0ufZD0FRbCP//JSYU7ueWZW5PVM3J78eXrJ7GjQglDRNJn7WXguZKSEi8t1axxu1m1CsaP555eB3Pjp65IVh9VvZ6//uwSOhXkRhiciETNzGa6e0lz7fQGd3s3YABMmcJXyqbzw5fuSla/mdOLy374EFt2aJY9EWmekkVHsN9+8PLL/NeKN7j25XuS1W/l9OTi/3uE9Ru2RRiciGQDJYuOYv/94eWX+eryqfzfS3cnq+fm9+D8G59kZdn6CIMTkbZOyaIjCRPGf5e9yS3P3EosfEpqaUF3Jv76BWa/vSjiAEWkrVKy6Gj23x+mTeOCimXc/tRN5FUHQ5mvyy/iwkcX8swcvektIrtTsuiIhgyB119nQlEl9z1yPd3KgynVd9Y4X39gFre+8AEJDW8uIimULDqqPn3g5Zc5asxAnhixnWG9Oic3/faFhXzpbzPYsF3vYohIQO9ZdHSJBMRibN5Rxf88MJNpiz9ObupfXMAfLzyYccN6RRigiGSS3rOQ9MSC/wWKC3O57ytH8PUT6+bDWL25ggv/NI0/3z5Zs+6JdHBKFpKUE4/x/QmjuedLJRR3ygGgOhbnlyviXPjd+1i6uP44kCLSUShZyG5OHt2Xp3dO5+APFyTrSvN6ccYdb/K3P/9Tnd8iHZCShTRoYK8uPPrwdXzr9QfJqakGoCInnxuWxrjoO3/j/feWRhyhiLQmJQtp2DXXkDv1db696R2evP877LduWXLTW/l9OPO+OfzkJ/ezZcv26GIUkVajZCGNKymBmTMZ+18XMfnv3+HrbzySvMqoicX5a3kPTv7RZB7527+pqUlEHKyIZJKShTQtPx9+9jPyp03l+xvf5t9//QZHL38nuXl9QVe+/36CM759P8899h/ay6PYIrIrJQtJz7hxMHMmI6//Dg/++xb+8NTN9N1a907GwsJeXDFjOxOve4Rpi9YraYi0M0oWkr7cXPjWt7CFCznr2P148e7/4erXH6Dzzh3JJm97Fy65azrn3j6N5+au0ZNTIu2EkoXsuT594C9/ocvUV7m6cD2v/vm/+MqMJ8kLR7EFmL1yE1fcP5PTfvcqj5SupKKqpokPFJG2TsN9yCf3yitw/fWU3fw7bl+Ty6OlZVTW6/DuVl3BhTnr+fwlJzH4oFERBSoi9aU73IeShbS4tVtVFzwrAAAQFklEQVQquHvqUh54cwXbdlbvss08wfFbVzJxTG9O+/zpdOrWNaIoRQSULKQN2FxexUNvreDv/5pJWbzzbtu7VO7gjKo1TDxmOEdO/BSx/LwIohTp2JQspG2orqZmyBD+U7AP9x12Jq8Ma/j/yX7bN3CabeD0w/bliHNPJrdYVxwirUHJQtqObdtg0iS44w5WLl7FkwecyONjT2ZpjwENNi+u2Mop5R/yqVE9Gf+ViRR3L2rlgEU6DiULaZtmz4YHHsAfeojZic48MeZk/rn/cWwsLG6weczgwIHdOH5kL44d0YtDB3cnL0cP8Ym0FCULadtqauDVV+GBB6h+/AneKhrAcyOP5rlRR/Fh1z6N7tYpN87Bg4oZt/w9SroZh40/iOLDDw3eARGRPaZkIdmjqgpeew0mT8YnT2bujhjPjTyKV048lznVBTT3Xt+o9Ss4rHoDY7rnMmZ4X/Yftz+dDjkQOnVqnfhFspiShWQnd5g3D/75T7jgAjb3G8S0xet5bdF6XvtgHSs3lDf7EbFEDcM2rGJM+XoO6OyM6FvEiOH9GTh2JPHR+0GXLq1wIiLZQclC2qU1myuYecffKX3hLWbtsz9z+w6jOp6T1r551ZUMpZwRh45meO/ODOvdhUE9ChnUvRO9O+dia9ZA//7JqWZFOoJ0k0V6/8pE2oh+xQWcec54zszdDG9Mpfzh3zObIub0H8HcPsOZ23cYS3oMIBGL77ZvZU4eC8hjwZzVu23LjxsD1q5g4NZ1DKrZwcC8BAOL8ujbvTN9enelT/9eFO7TF/r1C5auXcGsNU5ZpE3QlYVkN3coK4NZs+Cdd+CddyifM5f522Fun2Es7L0vi3oOYnGPgawt6vmJDtVl5w76bNtA7+0b6VO+mT5U0jsnQY8jDqX7mafTrTCX7oW5dCvMo1unXHJwKC2FHj2CpWtXdcRLm9MmbkOZ2QTgViAO3OXuN9Xbng/cB4wDPgYudPdl4bZrgcuBGuCb7j6lqWMpWcgutmyBd98N+j8WLoQFC9i8ZDlL/nw/iwp6sHjddpau30bZxnLKNpazubyqxUMoyovTbe0qupdvoVv5Vop2bqeoeiddqKGLJegSd7rEoUtejC55OXTplEuXTnl06ZxP4ZjRdPrceeTnxIjF6l3BzJgBOTlBB36nTlBQULeek6MrHtkjkScLM4sDC4FTgTJgBnCxu89LafN14CB3/5qZXQSc6+4XmtkBwEPAEcA+wAvAKHdvdOhSJQv5JLbceQ+rfvV7VlJAWXFfyor7sKprH9Z26cHaLt1Z17kHlTnRXBXk58QoyI3TKTdOQW6MgoXvU7BzBwVVlXSq3klBdSWdqnaSX72T/Joq8nFycXJJkGuEi5Pbpxe53/wGefEYufEYuXEjNycWlO++i9yVy8mNx8jLjROLx4nnxonHYsESN2LxGDkxIx6LEYvHyYnHiMVjxP/3+8Rzc7D6Ser114NkHY/vvsRidT/Ndl/Gj4du3Xb/j/Hhh7BgQV27xvZvaNvgwdCr1+6fWVEBS5fuvn+t2vWG6oqKoHfvhr+45ct3b5/Oem5uMLJzQz7+GKqrd90nFmv4vNLUFpLF0cAN7n56WL4WwN1/mdJmStjmDTPLAdYAvYEfpLZNbdfY8ZQspEVs3QorVgT/0FeuhI8+gjVr8DVr2Lx+E2u3VLC2vIa1uZ1Z27kH67p0Z+OnJrCp/2A27qhk044qNu6oZHN5Fe3kDm/aYgbxmBEzIydmxHZWEK8oJydRQyyRIO4JYp4g5k7ME5hDzBOAE3PHan+6w8iRxAoLg9+FZuHvb8PWriW2ZDGW0h6vtz/BZ5jvWo4deCAMGUKsNh8QfC6bNmEvvACAkfKlhV9gagq0ZF34c/hwOPLIXdrVJk27//4GPqvu85Of4buW6dkTO+ez4Wexy08efwLWr9tlPwo7YZdeytkH78ORw/b8Vmtb6OAeAKxMKZcBRzbWxt2rzWwz0DOsf7Pevg2PDSHSkoqKYMyYYElhQLdwGeUeJJUwkTByZNDpnaIm4WyZOp1N1/2IjeXVbKpMsLXG2JZbwLb8QrbndWJrfiHb8hpeLy8qpjy3gJ3V2TO3ecIhUeOAsxPAcqHTXl6Nba6GzVsa2JAHA/bfu88sB+Z/1PC20cfu3WcCzP6w4fqxJ+/9Z5aubLh+YAkMbKB++gr27991r5JFujKZLBq6cVr/b63G2qSzL2Z2BXAFwODBg/c0PpG9YxZ0VnftGiSKBsRjRvfjjqL7a88xtLbSHcrLg0SzbVvwM3VJ1n0MR/SH008nkXB2VieoqKqhvKqGim07KP/SV6ioTlBR7VTUBEt5wqhwoyKWQ3Ush6p4sFTGc6iK51IVy6Gq/z5UfuZsqmoSyaWyxqmqTlBVOpOqisqwfQ4Ji1ETi1FjMWpicWosRiIWo7p2vXZ7LE5NTm6zL05K9stksigDBqWUBwL1U3Btm7LwNlQxsCHNfXH3O4E7IbgN1WKRi2SCGRQWBkvfvmntEosZnfLidMqL0x2gWyd44anGd6iuDu7BV1XVLZWVwc/cXBg6tOH9hm2DDRt2bV9VGQzLkrokEruW3eGnP8XdqUk4Ne4kElCdSJB4aBI1b7xJTU0i2JZIkEgkSNQ47gkSNQncHffgL8GEg7uTAPx//xcfOCjc5nXbXnsNHpoUtHHCn3Wfkfwcgm24J8t+4UUkjjoqaOue/OvTl6+A3/w6KHuwj4f3fVJv09f/BeMYHH00XHABtZ9W29wduOaaXfbzXfZllz+Ld9k2YCBcfXXdfql3xn73O1hVFhw73OBdu8L113PE0B4Nf7ctJJN9FjkEHdynAKsIOrgvcfe5KW2uBA5M6eCe6O4XmNkY4EHqOrhfBEaqg1tEpGVF3mcR9kFcBUwheHT2Hnefa2Y3AqXuPhm4G7jfzBYRXFFcFO4718weAeYB1cCVTSUKERHJLL2UJyLSgaV7ZaFBcEREpFlKFiIi0iwlCxERaZaShYiINEvJQkREmtVunoYys3XA8mYbNq4XsL6FwolSezkP0Lm0Re3lPEDnUmtfd29kNMQ67SZZfFJmVprO42NtXXs5D9C5tEXt5TxA57KndBtKRESapWQhIiLNUrKoc2fUAbSQ9nIeoHNpi9rLeYDOZY+oz0JERJqlKwsREWlWh08WZjbBzBaY2SIz+0HU8XwSZrbMzOaY2Wwzy6pRFc3sHjNba2bvpdT1MLPnzeyD8Gf3KGNMRyPncYOZrQq/l9lm9ukoY0yXmQ0ys5fNbL6ZzTWzb4X12fi9NHYuWfXdmFmBmb1lZu+E5/GTsH6omU0Pv5OHzSyvxY/dkW9DmVmcYM6NUwkmXJoBXOzu8yINbC+Z2TKgxN2z7tlxMzse2Abc5+5jw7pbgA3uflOYyLu7+/9GGWdzGjmPG4Bt7v7/ooxtT5lZf6C/u88ysyJgJvBZ4Etk3/fS2LlcQBZ9NxZM8N3Z3beZWS7wOvAt4BrgcXefZGZ/At5x9zta8tgd/criCGCRuy9x90pgEnBOxDF1SO7+KsGcJqnOAe4N1+8l+MfdpjVyHlnJ3Ve7+6xwfSswHxhAdn4vjZ1LVvHAtrCYGy4OnAw8GtZn5Dvp6MliAJA6M3oZWfg/UAoHnjOzmeH85Nmur7uvhuAfO9An4ng+iavM7N3wNlWbv21Tn5kNAQ4FppPl30u9c4Es+27MLG5ms4G1wPPAYmCTu1eHTTLye6yjJwtroC6b78uNd/fDgDOAK8NbIhK9O4DhwCHAauDX0YazZ8ysC/AYcLW7b4k6nk+igXPJuu/G3Wvc/RBgIMHdkf0batbSx+3oyaIMGJRSHgh8GFEsn5i7fxj+XAs8QfA/Ujb7KLzXXHvPeW3E8ewVd/8o/AeeAP5CFn0v4X3xx4AH3P3xsDorv5eGziWbvxt33wT8BzgK6GZmtdNkZ+T3WEdPFjOAkeGTBHkEc4BPjjimvWJmncOOO8ysM3Aa8F7Te7V5k4HLwvXLgKcijGWv1f5iDZ1LlnwvYWfq3cB8d/9Nyqas+14aO5ds+27MrLeZdQvXOwGfIuh/eRn4XNgsI99Jh34aCiB8VO53QBy4x91/HnFIe8XMhhFcTQDkAA9m07mY2UPAiQSjZ34E/Bh4EngEGAysAM539zbdedzIeZxIcJvDgWXAV2vv+bdlZnYs8BowB0iE1dcR3OvPtu+lsXO5mCz6bszsIIIO7DjBH/uPuPuN4b//SUAP4G3gC+6+s0WP3dGThYiINK+j34YSEZE0KFmIiEizlCxERKRZShYiItIsJQsREWmWkoVkNTPrmTJi6Jp6I4hOa8U4hqSONNuazGxbI/WdzOyVcMDMxvZ9IRuGuJDo5TTfRKTtcvePCZ6Tz9rRXTPoKwQjkdY00eZ+4OtA1ryTI9HQlYW0W7V/cZvZieFf2I+Y2UIzu8nMPh/OCzDHzIaH7Xqb2WNmNiNcxjfwmXEz+1W4/V0z+2oDbYaY2WtmNitcjkmJ41Uze8LM5pnZn8wsFn7m38zsvTCeb4fth5vZs+HAkK+Z2eiwfqiZvRHG8NMm/hN8nvBNXjPrHx57dnic48I2kwleTBNpkq4spKM4mGDAtQ3AEuAudz/CgklwvgFcDdwK/NbdXzezwcAUdh+k7XJgs7sfbmb5wFQze45dB25bC5zq7hVmNhJ4CCgJtx0BHAAsB54FJgJLgQEp8190C9veCXzN3T8wsyOB2wmGor4VuMPd7zOzKxs62XD4mmHuviysugSY4u4/D29LFQK4+0YzyzeznuFVmkiDlCyko5hRO4yDmS0Gngvr5wAnheufAg4IhhECoKuZFYXzH9Q6DTjIzGrH4SkGRhJMolUrF/ijmR0C1ACjUra95e5LwjgeAo4FXgSGmdkfgKcJhpnvAhwD/CMlnvzw53jgvHD9fuDmBs63F7Ap9fyBe8LB9J5099kp29YC+wBKFtIoJQvpKFLHyUmklBPU/TuIAUe7e3kTn2PAN9x9yi6VwRwJtb5NMC7UweFnVqRsqz++jod/3R8MnA5cSTB729UEcxQc0kgczY3TUw4UpBzk1XDI+jOB+83sV+5+X7i5IGwv0ij1WYjUeQ64qrYQXhnUNwX4n/AvdMxsVDjKb6piYHU47PUXCQZ9q3VE2OcQAy4EXjezXkDM3R8DrgcOC+daWGpm54fHsTChAEwlGCEZgn6J3bj7RiBuZgXh/vsCa939LwSjrx5W+7lAP4JB9EQapWQhUuebQEnYcT0P+FoDbe4C5gGzwkdl/8zuV+i3A5eZ2ZsEt6C2p2x7A7iJYCjspQQjBQ8A/mPB7Gd/A64N234euNzM3gHmUjfl77cIJreaQZCYGvMcwW0uCEa+nW1mbxPcwro1rB8HvJkyy5pIgzTqrEgrMbMTge+6+2da6XiHAte4+xebaHMrMNndX2yNmCR76cpCpJ1y97eBl5t6KQ94T4lC0qErCxERaZauLEREpFlKFiIi0iwlCxERaZaShYiINEvJQkREmqVkISIizfr/0ULuWfrr21AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x205a3bc3ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write Euler method solution code here...\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "\n",
    "# initialise mass and time arrays + add starting boundary conditions.\n",
    "massEuler = np.empty(301)\n",
    "time = np.empty(301)\n",
    "\n",
    "massEuler[0] = mass_initial\n",
    "time[0] = t_start\n",
    "h = 0.1\n",
    "\n",
    "for i in range(300):\n",
    "    massdash = -r_decay*massEuler[i]\n",
    "    time[i+1] = time[i] + h\n",
    "    massEuler[i+1] = massEuler[i] + h*massdash\n",
    "\n",
    "# plotting the array values from the Euler method and the true solution\n",
    "true_masses = np.array([0.1*np.exp(-0.3*t) for t in time])\n",
    "\n",
    "plt.ylabel('Mass of Sample (kg)')\n",
    "plt.xlabel('Time elapsed (s)')\n",
    "plt.plot(time, massEuler, lw=4, c='r',linestyle='--')\n",
    "plt.plot(time, true_masses, lw=3)\n",
    "\n",
    "# Note the values aren't exactly the same as can be seen when you inspect a specific increment.\n",
    "print(massEuler[150])\n",
    "print(true_masses[150])\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-792f1290ba03a296",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "#### Task 1A: estimating the accuracy of a numerical solution\n",
    "\n",
    "You should write the difference between the array $M(i)$ and\n",
    "the exact solution, as a function of time, to an array, and plot it to see how big\n",
    "the error is. You should then decrease $h$ (remembering\n",
    "to increase the size of the array).\n",
    "How big is the largest difference between your numerical solution, and the\n",
    "exact solution, for 300 steps and $h=0.1$~s? How big is the largest\n",
    "difference when you reduce the number of steps to 100 (and so $h=0.3$~s)?\n",
    "\n",
    "You should\n",
    "find that the error varies approximately\n",
    " linearly with $h$. So the smaller is $h$ the more accurate is\n",
    "the solution, but of course the slower the program runs. This is not such a big problem\n",
    "for ODEs, but it is much more serious for PDEs.\n",
    "\n",
    "Once you have the error for your chosen value of $h$, it is good idea to write\n",
    "out the estimated accuracy of the solution to screen (just as you would quote error\n",
    "bars in a value you determine in an experiment), for example write to the screen\n",
    "\n",
    "${\\tt solution~for~mass~estimated~to~be~accurate~to~3~decimal~places~for~timestep~=0.01}$\n",
    "\n",
    "Note that the error on some points will be bigger than on other points, the error you write\n",
    "to the terminal should be the largest error in any point.\n",
    "\n",
    "You should also adjust the formatting of writing numbers\n",
    "such as ${\\tt mass}$, ${\\tt time}$, etc to the screen, to be around 1 decimal place\n",
    "more than the accuracy. So if the error is estimated to be in the 4th decimal place\n",
    "then a good format would be ${\\tt \\%9.5f}$, which writes numbers to 5 decimal places.\n",
    "As the accuracy will be greater for some points than for others, you will typically want \n",
    "to write numbers with a number of decimal places large enough to be larger than\n",
    "the accuracy of most or all points. However, this will vary a little from one problem to\n",
    "another so some judgement is often needed here.\n",
    "\n",
    "You should use the package pyplot from matplotlib to do the relevant plotting that was imported earlier. The package is fully documented at: https://matplotlib.org/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-5ab8359f1e7f39ab",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x205a3fea860>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecVPX1//HXm6WLqDQVQamK2FA3qDFGYomoKLEhltijRjFGY0+xYhQTaxRjAaNGwWiMqPmBsUWxAoIF1LgCKhZABUREcOH8/vjc/e6w7szeXebOndk9z8fjPmbmzp0757qyZ++nnI/MDOeccy7fmqUdgHPOucbJE4xzzrlEeIJxzjmXCE8wzjnnEuEJxjnnXCI8wTjnnEuEJxjnnHOJ8ATjnHMuEZ5gnHPOJaJ52gGkqVOnTtajR4+0w3DOuZIybdq0z82sc13HNekE06NHD6ZOnZp2GM45V1IkfRDnOG8ic845lwhPMM455xLhCcY551wiPME455xLhCcY55xzifAE45xzLhGeYJxzziWiSc+DKUlLl8K4cfDMMzB/PrRvD+XlMGwY9O2bdnTOOfd/PMGUmjvvhLPOWnPfv/4Fv/tdSDI33AAbbZRObM45l8ETTKn5xS/C3cu++0Lv3rBwITzxBIwfD9OnQ9u2aUfonHMAyMzSjiE15eXlVvSlYszgttvgyCNh3XWzHzdvHixbBltsUf05qTAxOueaFEnTzKy8ruO8k7/Y/fWvcOqpsNdesHp19uO6dVszufz613DllYWJ0TnnauFNZMXs7ber+1vOPBOaxfx7YMoUuOmmkGhat4azz04uRuecy8LvYIqVGZxyCnz7LRx3XGgii2vgQLjrrvD8nHPgsceSiNA553LyBFOsxo2D55+Hzp3huuvq//ljjoHLLguJ6sgjYebM/MfonHM5eIIpRpWV8Ic/hOcjR8L66zfsPFVDl5cuhcMOg2++yV+MzjlXB08wxehf/4KKCujTB44/vuHnkWDsWOjXL/TnXHJJ3kJ0zrm6eCd/MTrkkDC3RYLma/kjatsW7r8fLr889Mc451yBeIIpRhLsvXf+zjdgADz0UP7O55xzMXgTWbFZvDjZ81dWwr//nex3OOccnmCKy2efhTpiBx+ce1JlQ61aBbvtBvvvD089lf/zO+dcBk8wxeSee2DFijC0OO6kyvooKwvJBeCXv4SVK/P/Hc45F/EEU0zGjQuPxx2X3Hecd14oKfPee2G2v3POJcQTTLGYMwdeew3atYN99knue1q2rJ64edllYU0Z55xLgCeYYvHPf4bH/fcP9cOStO++4Xu++gouuijZ73LONVmJJhhJgyW9K6lC0gW1vN9K0vjo/Vck9ch478Jo/7uS9qnrnJLukjRH0oxoG5DkteVdVYI55JDCfN+110KLFmEi5uuvF+Y7nXNNSmLzYCSVATcDewPzgCmSJpjZrIzDTgQWmVkfScOBq4HDJfUHhgNbAV2BJyVtHn0m1znPNbMHk7qmxHzxBbzySrhz2Xffwnzn5puHmf2dOsHWWxfmO51zTUqSEy0HAhVmNhtA0jhgKJCZYIYCl0TPHwT+IknR/nFmtgKYI6kiOh8xzll6OnaEjz4KK1K2a1e47/XmMedcgpJsItsE+Cjj9bxoX63HmFklsATomOOzdZ1zpKQ3JF0nqVU+LqJgNt4Y9tsvve+fPz/Mk3HOuTxJMsHUtl5vzfWZsx1T3/0AFwL9gB8AHYDzaw1KOlnSVElTFy5cWNshhVUMS1bfeCP06hVqljnnXJ4kmWDmAd0zXncDPsl2jKTmwHrAlzk+m/WcZvapBSuAsVQ3qa3BzG4zs3IzK+/cuXMDLy2PJk8OVZNHjkwvhvbtQyn/3/8+TPR0zrk8SDLBTAH6SuopqSWh035CjWMmAMdGzw8FnjYzi/YPj0aZ9QT6Aq/mOqekjaNHAT8D3krw2vLnqafg/ffTnY/y859D//4wdy7cdlt6cTjnGpXEEkzUpzICmAS8DTxgZjMlXSbpwOiwO4GOUSf+2cAF0WdnAg8QOu8nAqeb2aps54zO9XdJbwJvAp2AK5K6trx6+unwuOee6cVQVgZXXhmejxzpC5M55/JCVgx9ACkpLy+3qVOnphfAsmWwwQahc/2LLxq+cmU+mMEPfgDTpoWZ/r/+dXqxOOeKmqRpZlZe13E+kz9Nzz8P330HO+6YbnKBsAbNxReH51dfDcuXpxuPc67keYJJUzE0j2UaMgR22CEssVwMI+yccyXNV7RM0/PPh8dBg1IN4/9IIemtt17akTjnGgFPMGm65JIwTHmnndKOpJonF+dcnuRMMJJaA0OA3Qg1wZYThv8+njF6yzXUPvskW5p/bcyYEdaLufnm5Ks7O+capawJRtIlwAHAs8ArwAKgNbA5cFWUfH5jZm8kH6YrKDM44YRQG23AADjjjLQjcs6VoKzDlCXtb2aPZ/2g1AXY1MxSHOe7dlIdpjx6dFiyeNiwUIes2PzrX3DQQdC1K8yeDa1Kq7Sbcy45az1MOVdyid5fUMrJJXVVc00+qVk9p0gMHRrK+H/yCdx7b9rROOdKUJ3DlCU9KmlCje0eSWdGzWSuvr78Et57L/RtbLtt2tHUToILovXcrr7aKy075+otzjyY2cDXwO3R9hUwn9AXc3tyoTVir74aHnfYIawqWawOPxx69AjJ8OGH047GOVdi4iSY7c3sSDN7NNqOBgaa2enADgnH1zi98kp4LKbhybVp3hzOPTc8/9Of0o3FOVdy4syD6SxpUzP7EEDSpkBVnfuViUXWmFXdwRR7ggE4/vhQZfm009KOxDlXYuIkmN8AkyW9T1jwqydwmqR1gL8lGVyj9frr4XGHErgBbNMGRo1KOwrnXAmqM8GY2b8l9SWsFingnbDbVgDXJxxf47NyZRid1bp1WEWylKxeDUuWhArQzjlXhzijyMaY2Qoze93MZgBlwL+TD62RatkSJk6EioqwDkupmD49jHg7/vi0I3HOlYg4nfwfSxoNIGkD4D+AT4xoarp2DStvPvIIzPQqQc65utWZYMzs98BXkm4FngD+bGZjE4+ssfrgg7DQWKnZcMNQPgbCvBjnnKtD1gQj6eCqDXgV2BmYDli0zzXEsGGw7rrw8stpR1J/554bmvXuuy8kSuecyyHXHcwBGdsQQnJpkfHa1deqVfDWW6GY5Oabpx1N/fXoAUccEa7D58U45+qQtdhlU1DwYpfvvRcSS7du8NFHhfvefHrrLdhmmzAK7oMPoEuXtCNyzhXYWhe7lPS7qFM/2/t7SPI7mfp4I1rZoFjrj8Wx9dZwwAFhYbJ33kk7GudcEcs1D+ZN4DFJ3wKvAQsJ68H0BQYATwJXJh5hY9IYEgzArbdChw6+EJlzLqesCcbMHgEeiSZZ7gpsTCh0eS9wspktL0yIjUhjSTBdu6YdgXOuBMSZyf8e8F4BYmn8Zs0Kj1ttlW4c+bJkCYwdG+qUtWyZdjTOuSITpxaZy5fHHw/9FqU4gqw2++0HL74Y+mN8hr9zroY4M/ldvvTpA0OGNJ6+i1NOCY+jRoU6Zc45l8ETjGu4I46A7t3DXdmjj6YdjXOuyMQpdrm5pKckvRW93lbS7+KcXNJgSe9KqpB0QS3vt5I0Pnr/FUk9Mt67MNr/rqR96nHOmyR9HSe+grr//tCM9OSTaUeSPy1awG9+E55fdVWYQOqcc5E4dzC3AxcC3wGY2RvA8Lo+JKkMuBnYF+gPHCGpf43DTgQWmVkf4Drg6uiz/aPv2AoYDNwiqayuc0oqB9aPcU2F99RTcNdd8O67aUeSXyedFIYsv/wyPP982tE454pInATT1sxerbGvMsbnBgIVZjbbzFYC44ChNY4ZSvWiZQ8Ce0pStH9ctEzAHKAiOl/Wc0bJ5xrgvBixFV5VYtlii3TjyLd11oERI8JzL4LpnMsQJ8F8Lqk3YACSDgU+jfG5TYDMeijzon21HmNmlcASoGOOz+Y65whggpnljE3SyZKmSpq6cOHCGJeRJ401wQCccQbsuSf88pdpR+KcKyJxhimfDtwG9JP0MTAHODrG51TLvpqN9NmOyba/toRokroChwGD6grKzG4jXA/l5eWF6TRYtAgWLoS2bWGTmjm2EejUqXH1LTnn8iLORMvZwF6S1gGamdnSmOeeB3TPeN0N+CTLMfMkNQfWA76s47O17d8e6ANUhBY22kqqiPp20ld197L55tCsCQzcMwPV9jeCc64pyZpgJJ2dZT8AZnZtHeeeAvSV1BP4mNBpf2SNYyYAxwIvAYcCT5uZSZoA3CfpWqArof7Zq4Q7m++d08xmAhtlxPh10SQXaNzNY5lmz4aLL4bOneHauv73cM41drnuYNZdmxObWaWkEcAkoAwYY2YzJV0GTDWzCcCdwD2SKgh3LsOjz86U9AAwizCg4HQzWwVQ2znXJs6C6NQJBg+GXXdNO5JkLVsG994LbdrAhReGROOca7J8PZhCrgfTFAwZEkri/OEPcOmlaUfjnEvAWq8Hk3GiXpIelbRQ0gJJj0jqlZ8wXaNz/vnh8aab4Ovim+/qnCucOD3O9wEPEMr1dwX+AdyfZFCNihnMmAGLF6cdSWH86Efwwx+GkXN33JF2NM65FMVJMDKze8ysMtru5fvDjV02n38O228Pm23WNEqpSNV3MddeC999l248zrnUxEkwz0i6QFIPSZtJOg94XFIHSR2SDrDkzZkTHnv1ajpDd4cMgf79Yd48mDw57WiccymJM9Hy8OjxlBr7TyDcyXh/TC6zZ4fHnj3TjaOQmjWD228Po8j69k07GudcSuJMtGxCvxkTkHkH05T88IdpR+CcS1mdCSYqIrk/0CPz+BgTLR1UJ5imdAeTyQxeeskTjnNNUJw+mEeB4whFKNfN2FwcVU1kTe0OBkJy+elPwwRT74txrsmJ0wfTzcy2TTySxqop38FIsMsuoRDm1VeHIczOuSYjzh3M/5P008QjaawmTQpbU7yDgVDKv00beOwxeOuttKNxzhVQnATzMvCwpOWSvpK0VNJXSQfWaPTpE5qJWrZMO5J0dO4MJ5wQno8alW4szrmCipNg/gzsQljZsr2ZrWtm7ROOyzUmv/kNlJXBfffB3LlpR+OcK5A4CeY94C1rylUxG+rxx+GYY+Dhh9OOJF09e8IRR8CqVXDVVWlH45wrkDid/J8Cz0r6f8CKqp0+TDmGl16Ce+4Jv2APOijtaNL129/CxIlNty/KuSYoToKZE20to83F1ZSHKNfUr18oHdOqVdqROOcKJM5Mfl/Uo6Ga8hDl2nhyca5JiTOTvzNwHrAV0Lpqv5ntkWBcjcOHH4bHzTZLN45i8s038Ne/wvLlcNFFaUfjnEtQnE7+vwPvAD2BS4G5wJQEY2ocvvsOPv00TDbs2jXtaIrH//4HZ58NV1wBCxakHY1zLkFxEkxHM7sT+M7M/mtmJwA7JxxX6fvkk1AqZaONoEWLtKMpHgMGwAEHhDuYP/857WiccwmKk2CqVoz6VNL+krYHuiUYU+NQWQn77Qd7eEvi9/z+9+Hx5pvDgmzOuUYpToK5QtJ6wG+Ac4A7gLMSjaox6N07zIO59960Iyk+P/gBDB4My5bB9denHY1zLiFqyvMny8vLberUqWmH0TRVlfBv3z7M7t9gg7Qjcs7FJGmamZXXdVzWOxhJv5DUN3ouSWMlLZH0RtRM5nL55BP44ovQD+O+b5ddYM894auv4O67047GOZeAXE1kZxJGjAEcAWxLWB75bODGZMNqBEaMgE6d4MEH046keF15JTzwQKi47JxrdHIlmEozq+rgHwLcbWZfmNmTwDrJh1biPvooPHbvnm4cxWzgQDjsMGgWpyvQOVdqcv3LXi1pY0mtgT2BJzPea5NsWI1AVYLp5gPuYpk7FxYtSjsK51we5UowfwCmEprJJpjZTABJuwOz45xc0mBJ70qqkHRBLe+3kjQ+ev8VST0y3rsw2v+upH3qOqekOyW9HvURPSipXZwYE7FyJcyfH0rUb7xxamGUjNGjYfPNw6qXzrlGI2uCMbPHgM2ALc3sFxlvTQUOr+vEksqAm4F9gf7AEZL61zjsRGCRmfUBrgOujj7bHxhOKE8zGLhFUlkd5zzLzLaLlnf+EBhRV4yJ+fjj8Ni1a0gyLrcddwyVD268ET77LO1onHN5krPx28wqzWxRjX3LzOzrGOceCFSY2WwzWwmMA4bWOGYo8Lfo+YPAnpIU7R9nZivMbA5QEZ0v6znN7CsII94ITXjpDd/y/pf6GTgQDjwwzO739WKcazSS7F3dBPgo4/W8aF+tx5hZJbAE6JjjsznPKWks8BnQD7gpHxfRIJ5g6u/yy8Pj6NHV//2ccyUt1zyYXaPHhtZYVy37at5VZDumvvvDE7Pjga7A22RpxpN0sqSpkqYuXLiwtkPW3j77wFNPwTnnJHP+xmjbbeHww0P/1RVXpB2Ncy4Pct3BVM11eamB554HZP4J3w34JNsxkpoD6wFf5vhsnec0s1XAeOCQ2oIys9vMrNzMyjt37lzPS4qpU6dQg6y8zomuLtMll4Qhy2PGwPvvpx2Nc24t5VoP5ruoyWkTSd+bWGlmv6rj3FOAvpJ6Ah8TOu2PrHHMBOBYQhI7FHjazEzSBOA+SdcS7kj6Aq8S7mC+d86o36W3mVVEzw8gLDHgSkm/fnDMMbB4cdqROOfyIFeCGQLsBewBTKvvic2sUtIIYBJQBowxs5mSLgOmmtkE4E7gHkkVhDuX4dFnZ0p6AJgFVAKnR3cmZDlnM+BvktoTktDrwC/rG3Pe/O53sHo1nHUWJHWX1Fjdfjs0j7OSt3Ou2NVZ7FLSdmb2eoHiKajEil1uuGFYTOvjj32xsbVhFhZsc84VlbUudpnhC0kPS1ogab6khyT59PRsKith4cLwi7FLl7SjKV2vvAI//jE891zakTjnGihOghlL6CvpShgS/Gi0z9Vm/vzwl3eXLt7UszYmTYLJk8NIPK9I7VxJipNgupjZ2GjSZaWZ3QV4x0I2n34aHr1EzNo5++yw3PSUKfCPf6QdjXOuAeIkmIWSjq4q1SLpaOCLpAMrWZ5g8qNdO7j00vD8wgthxYp043HO1VucBHMCMIwwQ/5TwnDiE5IMqqR5gsmfE06ALbeE2bPh1lvTjsY5V091Jhgz+9DMDjSzzmbWxcx+ZmYfFCK4krTOOrDddqE6sFs7zZtXV1i+7DKfH+NcifGVnvLtqKNgxgw4//y0I2kchgyB3XcPJWSm1Xs6lnMuRT7MyRU3Ce64A9ZdN8wvcs6VDE8w+fbVV+GXoU8QzJ8+fdKOwDnXAHU2kUk6U1J7BXdKek3STwsRXEnaaito3dpLziehshJuugleeCHtSJxzMcQaRRYt5vVTwvyX4wFfFao2q1eHFRlXrvQaZEm4+Wb41a9gxAhYtSrtaJxzdYiTYKraevYDxkZ1ybz9pzZffBH+yt5gg3AX4/LrF78Ii7jNmBFK+jvnilqcBDNN0hOEBDNJ0rrA6mTDKlE+ByZZbdvCNdeE5xdd5MOWnStycRLMicAFwA/M7BugBaGZzNXkCSZ5w4aFIpiff1490985V5TiJJhdgHfNbHFUJuZ3wJJkwypRn0SLa3qCSY4EN9wQVr78y19g1qy0I3LOZREnwYwGvpG0HXAe8AFwd6JRlSq/gymMAQNCf0xlpZeQca6IxZkHUxktYzwUuMHM7pR0bNKBlaSDD4Zu3UL9LJeskSNhhx3gxBPTjsQ5l0WcBLNU0oXA0cCPJZUR+mFcTf36hc0lr2NHOPnktKNwzuUQp4nscGAFcKKZfUZYdOyaRKNyrj4+/BBuvz3tKJxzNdR5BxMllWszXn+I98HU7k9/CvNfjj8+VFV2yfvqq9Ans3gxbLMN7Lxz2hE55yJxSsXsLGmKpK8lrZS0SpKPIqvNJZfAGWf4LPNCat8eTjklLKt86qmh4985VxTiNJH9BTgCeA9oA5wE3JxkUCXpm29g2TJo2TIUu3SF8/vfQ48e8PrrcP31aUfjnIvEWg/GzCqAMjNbZWZjgUGJRlWKFi4Mj126eCXlQmvbNtQpg5BsKirSjcc5B8RLMN9IagnMkDRK0lmAdzDUtGBBeOzSJd04mqr99oOjj4Zvv4WTTgqFR51zqYqTYH4OlAEjgGVAd+CQJIMqSZ5g0nf99eG//wsv+OqXzhWBOKPIPoieLge8+FM2nmDS17Ej3HtvWCphwIC0o3GuycuaYCS9keuDZrZt/sMpYatWQadOvqxv2vbeO+0InHORXE1kq4FVwD3AMOCAGludJA2W9K6kCkkX1PJ+K0njo/dfkdQj470Lo/3vStqnrnNK+nu0/y1JYyQVttrASSeFjv6rry7o17ocHnoIxo9POwrnmqysCcbMBhCGJ7cD7gNGAlsBH2c0m2UVlZS5GdgX6A8cIal/jcNOBBaZWR/gOuDq6LP9geHR9w0GbpFUVsc5/w70A7ahejh14fkIsuLw7LNw6KGhKOacOWlH41yTlLOT38zeMbOLzWwH4FHCDP6zYp57IFBhZrPNbCUwDhha45ihwN+i5w8Ce0pStH+cma0wszlARXS+rOc0s39bBHgV6BYzTtcY7b57KD66dCkcc4xPfnUuBTkTjKRNJP1G0mRCscuzCOX749gE+Cjj9bxoX63HmFklYZ2Zjjk+W+c5o6axnwMTs1zTyZKmSpq6sGruSj786EfQs6evT1IsJPjrX2GjjWDy5OqVMJ1zBZM1wUj6L+GupQVwHHAs8DjQUlKHGOeura3IYh5T3/2ZbgGeM7PnawvKzG4zs3IzK+/cuXNthzTM3Llha9cuf+d0a6dTJxg7Njz/wx9g+vR043Guicl1B7MZsAFwCvAEMDXapkWPdZlHmDNTpRvwSbZjJDUH1gO+zPHZnOeUdDHQGTg7Rnz5Y1Y9TDmfScutvcGD4fTT4bvv4KijYPnytCNyrsnI1cnfw8x6RluvjK2nmfWKce4pQF9JPaNKAMOBCTWOmUC4MwI4FHg66kOZAAyPRpn1BPoS+lWynlPSScA+wBFmVthp3EuWhF9g664LbdoU9KtdDKNGhXV6Vq6EefPSjsa5JiPOgmMNYmaVkkYAkwiVAMaY2UxJlwFTzWwCcCdwj6QKwp3L8OizMyU9AMwCKoHTzWwVQG3njL7yVsJyzi+FcQL808wuS+r61uCTLItb27bw6KOhP8abMJ0rGIUbhqapvLzcpk6N09pXh8mTYbfdYJdd4MUX1/58LnnffBMSj3Ou3iRNM7Pyuo6LVU3Z1cHvYErH8uVh3ZhBg0KTmXMuMbESTDTJsaukTau2pAMrKf37wx//CEcemXYkri7ffgsTJ8KUKXD++WlH41yjVmcTmaQzgIuB+YTyMQDWGGqR5a2JzJWWV14J85YqK+Gf/4SDDko7IudKSj6byM4EtjCzrcxsm2gr+eTimrCddqquGXfssfDOO+nG41wjFSfBfESYYe+yefLJ8JdwVV+MK35nnQWHHRZKyfzsZ2GouXMur+IkmNnAs1F147OrtqQDKykjR8Ihh8Cbb6YdiYtLgjFjYJtt4N13w0x/51xexZkH82G0tYw2V5OPIitN7drBww/DFVeEzTmXV3FWtPRVLOtSVTTTy8SUnt69q+uVOefyKlexy+ujx0clTai5FS7EIrd6NXzxRXjesWO6sbi1s2JFWD/mtdfSjsS5RiHXHcw90eOfChFIyVqyJCSZ9u2hRWEX0XR5dt11cMcd8PjjYShz9+51f8Y5l1XWBGNm06LH/xYunBLkdy+Nx1lnhUmY//0v7L9/KAHUvn3aUTlXsrxUzNr68svw6Amm9LVqFYabb7FFGBE4bFioku2caxBPMGtr4MBQ32rSpLQjcfnQoUNoIuvUKfxMTzstrPfjnKu3upZMLpPka83WpXXr8IvJNQ69e8OECeHnescdcOedaUfkXEnKmWCiNVh2VLTAinNNxi67wD/+AYceCkcfnXY0zpWkOE1k04FHJP1c0sFVW9KBlYwbbwyFE8ePTzsSl29DhoQk07p1eO1NZc7VS5wE0wH4AtgDOCDahiQZVEmZNQteeKF6NJlrnJYvD1WXH3ww7UicKxlxZvIfX4hASpYPU24axo+HRx4JAwDKyrzEv3Mx1HkHI6mbpIclLZA0X9JDkroVIriS4AmmaTj2WDjvvLCGzLBhYTizcy6nOE1kY4EJQFdgE+DRaJ+D6gTTqVO6cbhkSXDVVdVJ5vDD4aGH0o7KuaIWJ8F0NrOxZlYZbXcBXtWxyuefh0e/g2n8qpLM+edXJxnvk3EuqzgJ5nNJR0dzYsokHU3o9Hdm3kTW1Ejwxz+GJLNqVRhl5pyrVZz1YE4A/gJcBxjwYrTPrVoFJ54YCl62bZt2NK5QqpJMv35wxBFpR+Nc0cqZYCSVAYeY2YEFiqe0NG8ON9+cdhQuDRIcd1z162XL4O674dRTw3vOuVgz+YcWKBbnSpNZuJM57TQ46SQvkOlcJE4fzAuS/iJpN0k7VG2JR1YKPv8cXn0VPv447UhcmiQ45ZQw43/MGNh3X1i8OO2onEtdnATzQ2Ar4DLgz9EWaxEySYMlvSupQtIFtbzfStL46P1XJPXIeO/CaP+7kvap65ySRkT7TFJhxgz/5z+w005hHRHXtO2/Pzz7LHTpAk89FWqZzZ6ddlTOpaquasrNgNFm9pMa2x51nTjqv7kZ2BfoDxwhqX+Nw04EFplZH8Iggqujz/YHhhMS22DglqpRbDnO+QKwF/BBnAvPCx9B5jLttFNYCXOrreCdd8LrF15IOyrnUlNXH8xqYEQDzz0QqDCz2Wa2EhjH9/tzhgJ/i54/COwZVW4eCowzsxVmNgeoiM6X9ZxmNt3M5jYw1obxBONq6tEDXnwRBg8OTajjxqUdkXOpidNE9h9J50jqLqlD1Rbjc5sAH2W8nhftq/UYM6sElgAdc3w2zjkLx2fxu9q0bw+PPgrXXgt//nPa0TiXmjgJ5gTgdOA5YFq0TY3xudrGatasd57tmPruj03SyZKmSpq6cOHC+nz0+3wWv8umefPQN9eyZXi9eDEcdhjMmZNuXM4VUJ0Jxsx61rL1inHueUD3jNfdgE+yHSOpObAe8GWOz8YM3INGAAAQ2ElEQVQ5Z13Xc5uZlZtZeefOa1nxxpvIXFwXXRTKyuy4oxfKdE1G1gQj6byM54fVeO/KGOeeAvSV1FNSS0Kn/YQax0wAjo2eHwo8bWYW7R8ejTLrCfQFXo15zsLxBOPiGjkSDjgAFi2CQw6Bk08OkzOda8Ry3cEMz3h+YY33Btd14qhPZQQwCXgbeMDMZkq6TFJVZYA7gY6SKoCzgQuiz84EHgBmAROB081sVbZzAkj6laR5hLuaNyTdUVeMa+2hh8Kooa23TvyrXInbYIOwnsyNN0KrVnD77VBeDjNmpB2Zc4mRZVkGVtJ0M9u+5vPaXpeq8vJymzo1TneSc3n0xhth5v+sWaGPpqICunev+3POFQlJ08ysvK7jctUisyzPa3vtnItr221hyhQ491xYvdqTi2u0ciWY7SR9RRi51SZ6TvS6deKRFbvFi+GXvwy/HEaNSjsaV2ratg2FUlevrt733HPw9NNw4YWhGc25Epe1D8bMysysvZmta2bNo+dVr1sUMsiitHBhmETnqxq6tdEs+idYWRkKZV56KWyzDTzxRLpxOZcHcebBuNosWhQeN9gg3Thc49C8eej433JLeO892GcfGDbMC6m6kuYJpqE8wbh82333MKrsqqtCE9o//gFbbAGXXQbffJN2dM7VmyeYhvIE45LQsmVYjvntt+Hgg8NcmdGj1+yrca5ExFky2dXGE4xL0qabhv69//43DChp1y7s//prmDw5NKH5ypmuyPkdTEN5gnGFsPvuMDSjCPmf/hQWNNt1V5g4Maym6VyR8gTTUBtvHP6R9+2bdiSuKenSJZQmeumlkGh22gkee8wTjStKWWfyNwU+k9+VpK+/hltvhWuugQULwr7ttw/LAwwalGpormmIO5Pf72CcKzXt2sE554TS/9deCxttBNOnw8qVaUfm3Bo8wTTUokWwYkXaUbimrG3bsObM7Nlw772w997V7x17LJx2Grz5ZnrxuSbPE0xD7bkntG4N06alHYlr6tq0gaOOqh5V9tln8Pe/h+HN224LO+8Md9wRmtacKyBPMA1VNYps/fXTjcO5mjbaCF5/PdTKa98+LCnxi1+EgSknnwzz5qUdoWsiPME0lA9TdsVsq63gllvg00/hrrvgRz8KdzBjxkCLjFKCCxb4CDSXGJ9o2RCrVsGSJeH5euulG4tzubRtG/pjjj02VAeYPBk23DC8t3p1WMK5ffuwyubQobDDDj6B0+WNJ5iGWLw4PK63HpSVpRuLc3FtuWXYqsydC8uXhyazWbPg8suhWzc48MCQbAYNCqVrnGsgbyJrCG8ec41Br16hCW3iRDj1VOjaNSSbW24JpWimTKk+dulSb0pz9eYJpiGqEkyHDunG4dzaatEiJJPRo+Gjj+DVV+G3vw0lanbeufq4Aw+Enj3DYIEHHqie4OlcDj6TvyEz+b/8Ep55JgwP3W+//AfmXDGprAzFNz/9dM39ffqEckknngi77ZZObC4VcWfyex9MQ3ToEDpFnWsKmjcPTWfTp8N//gNPPhlqoVVUhG3QoOoE89RT8OyzYfDAjjuGPh0fNNBkeYJxztWtWbPqpHHBBfDdd2GuzQsvwF57VR/3yCNw003Vrzt3rv7czjvDkCGFj92lxpvIGtJE9swz4S+4PfZYs53auabumWdg0qRQ4WLatOr+SoBddoEXXwzPKyvhmGOgX78wZ2fLLcOgg9at04nb1Ys3kSVp4kQYNQquvNITjHOZfvKTsEEYdTZ3bnWy2WST6uPefx/uv3/Nz0qhSa1377DuzY47hv2ffx6GS7dvX5BLcPnjCaYhfJiyc3WTwsiznj3h0EPXfK9z51BhYOZMeOst+N//QjL66KOwNc/41XTxxWHodPv20L37mtvWW8NBBxXyqlw9eIJpCB+m7Nza6dAhVBfI9N138OGHYeDA5ptX71+5MozY/OqrkJBmzqx+b9Cg6gSzdGmYy7PhhmFhtg03XHPbbz/o0SMcu2JFSGI+UTpRiSYYSYOBG4Ay4A4zu6rG+62Au4EdgS+Aw81sbvTehcCJwCrgV2Y2Kdc5JfUExgEdgNeAn5tZMgtk+B2Mc/nXokVoHuvde839t98Ot90WpgdU3eFUbZttVn3cggWh3trXX4cmuJoee6w6wYwcGSoXrLde+HfcoUN43GCDcMw111R/7pFHQmzrrrvm1q5dSHw+Si6rxBKMpDLgZmBvYB4wRdIEM5uVcdiJwCIz6yNpOHA1cLik/sBwYCugK/CkpKo/abKd82rgOjMbJ+nW6NyjE7k4TzDOFZYUloru2BEGDKj9mF69wr/N+fNDspk/f80t867o22/DOZcsCdvcudXvbbHFmgnmqKNg2bLav3PUKDj33PB84sTQnNeuHayzTkg+mdsf/witWoVjH3kkJMyax7RpE+62uncPx1VWhmtq2bJ6K6G7riTvYAYCFWY2G0DSOGAokJlghgKXRM8fBP4iSdH+cWa2ApgjqSI6H7WdU9LbwB7AkdExf4vO6wnGuaZCCstnrL9+SBK5jBoFV10VksuiReGX/ZdfhueZ1abNYP/9Q/Pc0qXf39q1qz523rxQCSGbq65a8/nLL9d+3HHHwdix4fnbb4c1fWpeZ1WyefJJGBj9arz8chg3rvq9Fi3CY/PmYVLsLbdUn+Ppp8Mo2IQlmWA2AT7KeD0P2CnbMWZWKWkJ0DHa/3KNz1YNQantnB2BxWZWWcvx+deuXdg8wThXupo1q24W69Wr9mMkGD8++zkyp3n87Gdh0MHSpfDNN2Fbvrx6y0xcBxwQkmDm+1VbZrOfWbhrW7ky9FGtWBH2rVgRtmYZ1b6qipbWZuHCNV+//37JJ5jaGiZrTrrJdky2/bXVTst1/PeDkk4GTgbYdNNNazukbm+80bDPOecal8z+l06dwhbHRRfFO27bbcMw7UyrVoWEs3JlaIqrcuml8KtfVb9Xta1aFZZtqGJWPZQ8YUkmmHlA94zX3YBPshwzT1JzYD3gyzo+W9v+z4H1JTWP7mJq+y4AzOw24DYIEy3rf1nOOZeisrLq/ppMG20UtrpIocmsAJKspjwF6Cupp6SWhE77CTWOmQBUjVU8FHjaQmmBCcBwSa2i0WF9gVeznTP6zDPROYjO+UiC1+acc64Oid3BRH0qI4BJhCHFY8xspqTLgKlmNgG4E7gn6sT/kpAwiI57gDAgoBI43cxWAdR2zugrzwfGSboCmB6d2znnXEq8FllDapE551wTFrcWmS845pxzLhGeYJxzziXCE4xzzrlEeIJxzjmXCE8wzjnnEtGkR5FJWgh80MCPdyJM8GwMGsu1NJbrAL+WYtVYrmVtr2MzM+tc10FNOsGsDUlT4wzTKwWN5Voay3WAX0uxaizXUqjr8CYy55xzifAE45xzLhGeYBrutrQDyKPGci2N5TrAr6VYNZZrKch1eB+Mc865RPgdjHPOuUR4gmkASYMlvSupQtIFacfTUJLmSnpT0gxJJVX1U9IYSQskvZWxr4Ok/0h6L3osiSVHs1zLJZI+jn42MyTtl2aMcUjqLukZSW9LminpzGh/yf1cclxLKf5cWkt6VdLr0bVcGu3vKemV6OcyPloCJb/f7U1k9SOpDPgfsDdhYbQpwBFmlmWt0uIlaS5QbmYlN65f0o+Br4G7zWzraN8o4EszuypK/BuY2flpxhlHlmu5BPjazP6UZmz1IWljYGMze03SusA04GfAcZTYzyXHtQyj9H4uAtYxs68ltQAmA2cCZwP/NLNxkm4FXjez0fn8br+Dqb+BQIWZzTazlcA4YGjKMTU5ZvYcYQ2hTEOBv0XP/0b4hVD0slxLyTGzT83stej5UuBtYBNK8OeS41pKjgVfRy9bRJsBewAPRvsT+bl4gqm/TYCPMl7Po0T/xyP8T/aEpGmSTk47mDzY0Mw+hfALAuiScjxra4SkN6ImtKJvVsokqQewPfAKJf5zqXEtUII/F0llkmYAC4D/AO8Di6Ml5iGh32OeYOpPtewr1XbGXc1sB2Bf4PSoqcYVh9FAb2AA8Cnw53TDiU9SO+Ah4Ndm9lXa8ayNWq6lJH8uZrbKzAYA3QitMFvWdli+v9cTTP3NA7pnvO4GfJJSLGvFzD6JHhcADxP+xytl86O286o29AUpx9NgZjY/+qWwGridEvnZRG38DwF/N7N/RrtL8udS27WU6s+lipktBp4FdgbWl9Q8eiuR32OeYOpvCtA3GoHREhgOTEg5pnqTtE7UeYmkdYCfAm/l/lTRmwAcGz0/FngkxVjWStUv5MhBlMDPJupMvhN428yuzXir5H4u2a6lRH8unSWtHz1vA+xF6FN6Bjg0OiyRn4uPImuAaGji9UAZMMbMRqYcUr1J6kW4awFoDtxXStch6X5gEKEq7HzgYuBfwAPApsCHwGFmVvSd51muZRChGcaAucApVf0YxUrSj4DngTeB1dHuiwh9FyX1c8lxLUdQej+XbQmd+GWEm4oHzOyy6HfAOKADMB042sxW5PW7PcE455xLgjeROeecS4QnGOecc4nwBOOccy4RnmCcc84lwhOMc865RHiCcU2OpI4Z1XA/q1Ed98UCxtEjs4JyIUn6Osv+NpL+GxV1zfbZJ0ulRIpLV/O6D3GucTGzLwhzGUqyanHCTiBU2F2V45h7gNOAkpk35dLhdzDOZaj6y17SoOgv+Qck/U/SVZKOitbVeFNS7+i4zpIekjQl2nat5Zxlkq6J3n9D0im1HNND0vOSXou2H2bE8ZykhyXNknSrpGbROe+S9FYUz1nR8b0lTYwKmD4vqV+0v6ekl6IYLs/xn+AoohndkjaOvntG9D27RcdMIEw4dC4nv4NxLrvtCEUBvwRmA3eY2UCFxafOAH4N3ABcZ2aTJW0KTOL7hQRPBJaY2Q8ktQJekPQEaxYXXADsbWbfSuoL3A+UR+8NBPoDHwATgYOBOcAmGevHrB8dextwqpm9J2kn4BZCWfYbgNFmdrek02u72Kj0US8zmxvtOhKYZGYjoyaztgBmtkhSK0kdo7tB52rlCca57KZUlQGR9D7wRLT/TeAn0fO9gP6hdBUA7SWtG60hUuWnwLaSquo+rQf0JSxcV6UF8BdJA4BVwOYZ771qZrOjOO4HfgQ8BfSSdBPwOGHZhXbAD4F/ZMTTKnrcFTgken4PcHUt19sJWJx5/cCYqOjjv8xsRsZ7C4CugCcYl5UnGOeyy6zLtDrj9Wqq/+00A3Yxs+U5ziPgDDObtMbOsM5IlbMIdci2i875bcZ7Nes5WXQXsR2wD3A6YaXFXxPW+BiQJY666kItB1pnfMlz0RIO+wP3SLrGzO6O3m4dHe9cVt4H49zaeQIYUfUiugOpaRLwy+hOAEmbRxWsM60HfBqVgf85oTBhlYFRH0oz4HBgsqROQDMzewj4PbBDtF7JHEmHRd+jKAkBvECo/A2hn+V7zGwRUCapdfT5zYAFZnY7obLwDlXnBTYiFHt0LitPMM6tnV8B5VHn/Szg1FqOuQOYBbwWDUv+K99vPbgFOFbSy4TmsWUZ770EXEUoDT+HUAV7E+BZhVUK7wIujI49CjhR0uvATKqX8z6TsKjcFEIyy+YJQhMchIrOMyRNJzSv3RDt3xF4OWM1ROdq5dWUnStikgYB55jZkAJ93/bA2Wb28xzH3ABMMLOnChGTK11+B+Oc+z9mNh14JtdES+AtTy4uDr+Dcc45lwi/g3HOOZcITzDOOecS4QnGOedcIjzBOOecS4QnGOecc4nwBOOccy4R/x8FSXSY8RnkLQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x205a52ff390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write mass error calculation and plotting code here...\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "mass_error = np.empty(301)\n",
    "\n",
    "for i in range(301):\n",
    "    mass_error[i] = abs(massEuler[i]-true_masses[i])\n",
    "  \n",
    "plt.ylabel('Error in Mass of Sample (kg)')\n",
    "plt.xlabel('Time elapsed (s)')\n",
    "plt.plot(time, mass_error, lw=2, c='r',linestyle='--')\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b67b488d737e4c00",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## 4) Solving ODEs using a Python ODE solver <a class=\"anchor\" id=\"ode-solver\"></a>\n",
    "\n",
    "One of the main advantages of using python is it's modularity. This coupled with it's wide use in the scientific community means that for many problems rather than having to create your own solution it is easy to import a module which has a function built in to do it for you. \n",
    "\n",
    "So in our case we are going to use the ${\\tt scipy.integrate}$ module as this has functions that will allow us to solve our radioactive sample problem. We will be using the odeint function this function solves a system of ordinary differential equations using lsoda from the FORTRAN library odepack. In order to run it requires us to input 3 things, firstly the function of the ODE, which we need to define, so in our case:  \n",
    "\n",
    "$$\\frac{{\\rm d}M(t)}{{\\rm d}t}=-0.3M(t)$$\n",
    "\n",
    "The next two inputs are the initial value of $M$ (so at $t=0$) and then the time array for which to solve for $M$. Additionally we also give it a value for rtol, this is an optional argument (unlike the previous three) which allows us to define the tolerance or error we are permitting for the solution from the ODE. If you plot the error and play around with this value you can see how it effects the directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f2046e67581e1a89",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Again we can check the accuracy of this by comparing it to the true values we have calculated as well as plotting it on a graph.\n",
    "\n",
    "You'll often find that when using Python there is already a function that does whatever you're trying to do, typically it is just a matter of adjusting your problem to be in the same format. It's rare you will have a problem in which you will need to undertake entirely using your own functionality.  \n",
    "\n",
    "### Task 2 <a class=\"anchor\" id=\"Task-2\"></a>\n",
    "\n",
    "Write a python program that utilises scipy.integrate in order to solve an ODE. Compare the errors of this solution to that of the Euler method. You should find that the errors in this method are far lower, that is because all of the functions that are used in numpy and scipy have been well written to optimise both speed and accuracy (although obviously there is a trade-off between the two). How does the error compare as you change the value of rtol? What happens if instead of using rtol (relative tolerance) you pass an input of atol (absolute tolerance)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b2c7478caaa185ed",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x205a58ad080>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEKCAYAAAC2bZqoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xuc1VW9//HX2xFF8ZKCFwQRNCzxBjqiZiGaAVmJFiZ0kchz6CLnmJa3X6dAPRXqSU9GapaiVopmqeQlEZVQ0wQE7xIcJR0hccC8hJfAz++P7xrcDHvv2czee2YP834+Hvsx3+/6ru/6ru9snQ9rfdd3LUUEZmZmtWST9q6AmZlZcw5OZmZWcxyczMys5jg4mZlZzXFwMjOzmuPgZGZmNaddg5OkEZIWSlos6aw8xzeXdEM6/hdJfXOOnZ3SF0oanpN+laTlkp5sVtb2ku6WtCj93C6lS9IlqazHJR1QvTs2M7NStFtwklQH/Az4JDAAGCNpQLNsJwGvRsQHgYuB89O5A4DRwN7ACODSVB7A1SmtubOAeyKiP3BP2iddv3/6jAcuq8T9mZlZ67Vny2kwsDginouId4FpwMhmeUYC16Ttm4CPS1JKnxYR70TE88DiVB4RMRtYmed6uWVdAxybk35tZB4GPiCpZ0Xu0MzMWmXTdrx2L+DFnP0G4OBCeSJitaTXgO4p/eFm5/Zq4Xo7RcSyVNYySTsWqUcvYFnzAiSNJ2td0a1btwM//OEPt3BJMzNrMm/evMaI2KGUvO0ZnJQnrflcSoXylHJuJeuRJUZcAVwBUF9fH3Pnzm3lJc3MOh9Jfys1b3t26zUAu+bs9waWFsojaVNgW7Iuu1LObe7lpu669HP5BtTDzMzaUHsGpzlAf0n9JG1GNsBherM804GxaXsUcG9kM9VOB0an0Xz9yAYzPNLC9XLLGgvcmpN+Yhq1dwjwWlP3n5mZtY9269ZLz5AmAHcBdcBVEfGUpHOBuRExHbgS+JWkxWQtptHp3Kck3Qg8DawGTo6INQCSrgeGAj0kNQATI+JKYDJwo6STgBeA41NV7gCOJhtUsQoYV/27NzOzYuQlM1rHz5zMasO//vUvGhoaePvtt9u7KpZ07dqV3r1706VLl3XSJc2LiPpSymjPARFmZmVraGhg6623pm/fvmRvmlh7ighWrFhBQ0MD/fr1a3U5nr7IzDq0t99+m+7duzsw1QhJdO/eveyWrIOTmXV4Dky1pRLfh4OTmZnVHAcnM7MyNTQ0MHLkSPr3788ee+zBKaecwrvvvgvArFmz2HbbbRk0aBAf+tCHGDJkCLfddtvacydNmkSvXr0YOHDg2s8//vGPsurTt29fGhsbyyqjvTk4mZmVISL47Gc/y7HHHsuiRYv461//yptvvsl3v/vdtXk+9rGPMX/+fBYuXMgll1zChAkTuOeee9YeP/XUU1mwYMHazwc+8IH2uJV1RATvvfdeu13fwcnMrAz33nsvXbt2Zdy47BXJuro6Lr74Yq666ipWrVq1Xv6BAwfy/e9/nylTppR8jWXLljFkyBAGDhzIPvvsw/333w/A9ddfz7777ss+++zDmWeeud55Z555Jpdeeuna/UmTJvHjH/8YgAsvvJCDDjqI/fbbj4kTJwKwZMkS9tprL775zW9ywAEH8OKLL65XZltxcDIzK8NTTz3FgQceuE7aNttsQ58+fVi8eHHecw444ACeffbZtfsXX3zx2i69I444Yr381113HcOHD2fBggU89thjDBw4kKVLl3LmmWdy7733smDBAubMmcMtt9yyznmjR4/mhhtuWLt/4403cvzxxzNjxgwWLVrEI488woIFC5g3bx6zZ88GYOHChZx44onMnz+f3XbbrdW/l3I5OJlZ59PYCBdemP0sU0TkHZ1WKL3pWK7cbr377rtvvfwHHXQQU6dOZdKkSTzxxBNsvfXWzJkzh6FDh7LDDjuw6aab8sUvfnFtgGkyaNAgli9fztKlS3nsscfYbrvt6NOnDzNmzGDGjBkMGjRobaBctGgRALvtthuHHHJIa38dFeOXcM2s85k6Fc44I9s+/fSyitp777353e9+t07a66+/zosvvsgee+zBihUr1jtn/vz57LXXXiVfY8iQIcyePZvbb7+dL3/5y5x++ulss802JZ07atQobrrpJv7+978zevRoIAuOZ599Nl/72tfWybtkyRK6detWcr2qyS0nM+t8xo2DCy7Ifpbp4x//OKtWreLaa68FYM2aNXz729/mK1/5CltuueV6+R9//HHOO+88Tj755JKv8be//Y0dd9yRf//3f+ekk07i0Ucf5eCDD+ZPf/oTjY2NrFmzhuuvv57DDz98vXNHjx7NtGnTuOmmmxg1ahQAw4cP56qrruLNN98E4KWXXmL58uXrndue3HIys86nR4+yW0xNJHHzzTfzzW9+k/POO4/33nuPo48+mh/+8Idr89x///0MGjSIVatWseOOO3LJJZfw8Y9/fO3xiy++mF//+tdr92+55Rb69u27dn/WrFlceOGFdOnSha222oprr72Wnj178qMf/YgjjjiCiODoo49m5Mjmi4lnLbs33niDXr160bNntsj3sGHDeOaZZzj00EMB2Gqrrfj1r39NXV1dRX4nleCJX1vJE7+a1YZnnnlmg7rIrG3k+142ZOJXd+uZmVnNcXAyM7Oa4+BkZmY1x8HJzMxqjoOTmZnVnHYNTpJGSFooabGks/Ic31zSDen4XyT1zTl2dkpfKGl4S2VKul/SgvRZKumWlD5U0ms5x75f3bs2M7OWtFtwklQH/Az4JDAAGCNpQLNsJwGvRsQHgYuB89O5A4DRwN7ACOBSSXXFyoyIj0XEwIgYCDwE/D7nOvc3HYuIc6t0y2a2kaqrq1tnyYvJkycXzX/11VczYcKEVl9vyZIlSOJ73/ve2rTGxka6dOnSYrmzZs3iz3/+89r9r3zlK9x0002trku55xfSni/hDgYWR8RzAJKmASOBp3PyjAQmpe2bgCnKJqsaCUyLiHeA5yUtTuXRUpmStgaOBMp/NdzMDNhiiy1YsGBB1cpfvXo1m2667p/r3Xffndtuu43zzjsPgN/+9rfsvffeLZY1a9YsttpqKz7ykY9Upa6V0p7der2A3PnYG1Ja3jwRsRp4Dehe5NxSyjwOuCciXs9JO1TSY5LulNTyt2tmVoLcRf/mzp3L0KFD18vzyiuv8LnPfY6DDjqIgw46iAcffBDIlrcYP348w4YN48QTT1zvvC222IK99tqLpskAbrjhBj7/+c8XLXfJkiVcfvnla2dBb1p6Y/bs2XzkIx9h9913X9sKighOP/109tlnH/bdd9+1s5tHBBMmTGDAgAF86lOfqtq0R+3Zcso3XW/z6SoK5SmUni/YNi9zDPDLnP1Hgd0i4k1JRwO3AP3zVlgaD4wH6NOnT74sZtYJvfXWWwwcOHDt/tlnn80JJ5xQ0rmnnHIKp556Kh/96Ed54YUXGD58OM888wwA8+bN44EHHmCLLbbIe27TvHk777wzdXV17LLLLixdurRouV//+tfZaqut+M53vgPAlVdeybJly3jggQd49tlnOeaYYxg1ahS///3v1y7R0djYyEEHHcSQIUN46KGHWLhwIU888QQvv/wyAwYM4Ktf/Wo5v7682jM4NQC75uz3BpYWyNMgaVNgW2BlC+cWLFNSd7Luv+Oa0nJbUBFxh6RLJfWIiPXm0o+IK4ArIJu+qLTbNLNa09iYTUw+blw2zV65yunWmzlzJk8//f7TjNdff5033ngDgGOOOaZgYAIYMWIE3/ve99hpp53WC4bFym3u2GOPZZNNNmHAgAG8/PLLADzwwAOMGTOGuro6dtppJw4//HDmzJnD7Nmz16bvsssuHHnkka2675a0Z3CaA/SX1A94iWyAwxea5ZkOjCUbwDAKuDciQtJ04DpJFwG7kLV0HiFrURUr83jgtoh4uylB0s7Ay6ncwWStr/XnuDezjUYFV8woatNNN1271Pnbb7+dN897773HQw89lDcItbR8xWabbcaBBx7Ij3/8Y5566in+8Ic/lFRuc5tvvvna7ab5VovNu1ponapKardnTukZ0gTgLuAZ4MaIeErSuZKOSdmuBLqnAQ+nAWelc58CbiQb6PBH4OSIWFOozJzLjgaub1aVUcCTkh4DLgFGh2fDNduoVXDFjKL69u3LvHnzANZb86nJsGHD1lmyfUNbYN/+9rc5//zz6d69e0nlbr311gVbULmGDBnCDTfcwJo1a3jllVeYPXs2gwcPZsiQIUybNo01a9awbNmyvIsjVkK7LpkREXcAdzRL+37O9ttkrZ185/4A+EEpZeYcG5onbQowZf3cZraxquCKGcD6z5xGjBjB5MmTmThxIieddBI//OEPOfjgg/Oee8kll3DyySez3377sXr1aoYMGcLll19e8rX33nvvvKP0CpX7mc98hlGjRnHrrbfy05/+tGC5xx13HA899BD7778/krjgggvYeeedOe6447j33nvZd9992XPPPfOuIVUJXjKjlbxkhllt8JIZtclLZpiZ2UbHwcnMzGqOg5OZdXh+PFFbKvF9ODiZWYfWtWtXVqxY4QBVIyKCFStW0LVr17LKadfRemZm5erduzcNDQ288sor7V0VS7p27Urv3r3LKsPBycw6tC5dutCvX7/2roZVmLv1zMys5hRtOUnqCnwa+BjZNEFvAU8CtzebecHMzKxiCgYnSZOAzwCzgL8Ay4GuwJ7A5BS4vh0Rj1e/mmZm1pkUaznNiYhJBY5dJGlHwOtGmJlZxRUMThFxe7ETI2I5WWvKzMysolocrSfpD6y/YN9rwFzg57nLT5iZmVVCKaP1ngPeBH6RPq8DL5M9e/pF9apmZmadVSnvOQ2KiCE5+3+QNDsihkjyiD0zM6u4UlpOO0haO/Ahbe+Qdt+tSq3MzKxTK6Xl9G3gAUn/R7YMej/gm5K6AddUs3JmZtY5tRicIuIOSf2BD5MFp2ez5HgH+N8q18/MzDqhFrv1JF0VEe9ExGMRsQCoo8Ay6GZmZpVQyjOnlyRdBiBpO+Bu4NeVuLikEZIWSlos6aw8xzeXdEM6/hdJfXOOnZ3SF0oa3lKZkq6W9LykBekzMKVL0iUp/+OSDqjEvZmZWeu1GJwi4nvA65IuB2YAP46IqeVeWFId8DPgk8AAYIykAc2ynQS8GhEfBC4Gzk/nDgBGA3sDI4BLJdWVUObpETEwfRaktE8C/dNnPHBZufdmZmblKRicJH226QM8AhwCzAcipZVrMLA4Ip6LiHeBacDIZnlG8v6gi5uAj0tSSp+WuhufBxan8kops7mRwLWReRj4gKSeFbg/MzNrpWItp8/kfD5NFpi65OyXqxfwYs5+Q0rLmyciVpPNTNG9yLktlfmD1HV3saTNN6AeAEgaL2mupLle2MzMrHqKza03rsrXVr7LlpinUHq+YNtU5tnA34HNgCuAM4FzS6xHlhhxRTqX+vp6rwltZlYlxbr1/isNgCh0/EhJ5bSgGoBdc/Z7A0sL5ZG0KbAtsLLIuQXLjIhlqevuHWAqWRdgqfUwM7M2VOw9pyeA2yS9DTwKvEK2nlN/YCAwE/hhGdeeA/SX1A94iWyAwxea5ZkOjAUeAkYB90ZESJoOXCfpIrJFEPuTPRdToTIl9YyIZemZ1bFkiyY2XWOCpGnAwcBrEbGsjPsyM7MyFevWuxW4Nb2AexjQk2zS118D4yPirXIuHBGrJU0A7iJ7d+qqiHhK0rnA3IiYDlwJ/ErSYrIW0+h07lOSbgSeBlYDJ0fEGoB8ZaZL/kbSDmQBbAHw9ZR+B3A02aCKVUC1uzPNzKwFivCjk9aor6+PuXPntnc1zMw6DEnzIqK+lLylvIRrZmbWphyczMys5jg4mZlZzSll4tc9Jd0j6cm0v5+k/6p+1czMrLMqpeX0C7IXWP8FEBGPk0bNmZmZVUMpwWnLiHikWdrqalTGzMwMSgtOjZL2IE3pI2kU4JdUzcysakpZpv1ksvnkPizpJeB54EtVrZWZmXVqpSzT/hxwlKRuwCYR8Ub1q2VmZp1ZweAk6bQC6QBExEVVqpOZmXVyxVpOW7dZLczMzHIUm/j1nLasiJmZWZNSXsLdXdIfJL0iabmkWyXt3haVMzOzzqmUoeTXATeSLZmxC/Bb4PpqVsrMzDq3UoKTIuJXEbE6fX5NgWXMzczMKqGU95zuk3QWMI0sKJ0A3C5pe4CIWFnF+pmZWSdUSnA6If38WrP0r5IFKz9/MjOziirlJdx+bVERMzOzJqWM1quTdIyk/5R0WtOnEheXNELSQkmLU9dh8+ObS7ohHf+LpL45x85O6QslDW+pTEm/SelPSrpKUpeUPlTSa5IWpM/3K3FvZmbWeqUMiPgD8BWgO9mLuU2fskiqA34GfBIYAIyRNKBZtpOAVyPig8DFwPnp3AFky3bsDYwALk1BtFiZvwE+DOwLbAH8W8517o+Igelzbrn3ZmZm5SnlmVPviNivCtceDCxOc/chaRowEng6J89IYFLavgmYomz+pJHAtIh4B3he0uJUHoXKjIg7mgqV9AjQuwr3ZGZmFVBKy+lOScOqcO1ewIs5+w0pLW+eiFgNvEbWgit0botlpu68LwN/zEk+VNJjku6UtHehCksaL2mupLmvvPJKy3doZmatUkpwehi4WdJbkl6X9Iak1ytwbeVJa/7+VKE8G5qe61JgdkTcn/YfBXaLiP2BnwK3FKpwRFwREfURUb/DDjsUymZmZmUqJTj9GDiUbEXcbSJi64jYpgLXbgB2zdnvDSwtlEfSpsC2wMoi5xYtU9JEYAdg7YCOiHg9It5M23cAXST1KOfGzMysPKUEp0XAkxFR6Vkh5gD9JfWTtBnZAIfpzfJMB8am7VHAvake04HRaTRfP6A/8EixMiX9GzAcGBMR7zVdQNLO6TkWkgaT/U5WVPhezcxsA5QyIGIZMEvSncA7TYnlrucUEaslTQDuAuqAqyLiKUnnAnMjYjpwJfCrNOBhJVmwIeW7kWzwxGrg5IhYA5CvzHTJy4G/AQ+lWPT7NDJvFPANSauBt4DRVQjEZma2AdTS3+HUFbaezr6kRn19fcydO7e9q2Fm1mFImhcR9aXkLWWGiE4dhMzMrO21GJwk7QCcQfbCa9em9Ig4sor1MjOzTqyUARG/AZ4F+gHnAEvIBh6YmZlVRSnBqXtEXAn8KyL+FBFfBQ6pcr3MzKwTK2W03r/Sz2WSPkX23pCn/jEzs6opJTj9t6RtgW+TzaCwDXBqVWtlZmadWimj9W5Lm68BR1S3OmZmZkWeOUn6d0n907YkTU3rHj0uaVDbVdHMzDqbYgMiTiEbmQcwBtiPbEn204BLqlstMzPrzIoFp9UR0TQY4tPAtRGxIiJmAt2qXzUzM+usigWn9yT1lNQV+DgwM+fYFtWtlpmZdWbFBkR8H5hLNoHq9KYJVCUdDjzXBnUzM7NOqmBwiojbJO0GbB0Rr+YcmgucUPWamZlZp1V0KHlaGv3VZmn/rGqNzMys0ytl+iIzM7M2Vew9p8PSz83brjpmZmbFW05N7zI91BYVMTMza1LsmdO/JE0Fekla76XbiPjP6lXLzMw6s2Itp08DdwFvA/PyfMomaYSkhZIWSzorz/HNJd2Qjv9FUt+cY2en9IWShrdUpqR+qYxFqczNWrqGmZm1j2JDyRuBaZKeiYjHKn1hSXXAz4BPAA3AHEnTI+LpnGwnAa9GxAcljQbOB06QNAAYTbY67y7ATEl7pnMKlXk+cHFETJN0eSr7skLXqPT9mplZ6UpZMmOFpJuBw4AAHgBOiYiGMq89GFgcEc8BSJoGjARyg9NIYFLavgmYIkkpfVpEvAM8L2lxKo98ZUp6BjgS+ELKc00q97JC14iIKPP+8mpcuIKpE+Yx7kN/pseWq8orbMstYcwYuP56WFVmWVaYf89tw7/nghpXbcmU+R9h1btdyi5ry25izFl9uP5Hf2vVr3nLbmLCz/elx4e6l12XYkoJTlOB64Dj0/6XUtonyrx2L+DFnP0G4OBCeSJitaTXgO4p/eFm5/ZK2/nK7A78I7231Tx/oWs0Nq+wpPHAeIA+ffqUep/rmHraE5wxcxjMvJvT+Z9WlbGOOXPgjjvKL8eK8++5bVT599xId6ZwMqs62PSg8xnITIZVrLw5Tz/CHa+0fgWkbqfN4vTbh1asPvmUEpx2jIipOftXS/pWBa6tPGnNWyuF8hRKz/cMrVj+UuuRJUZcAVwBUF9f36qW1biL9oV3ZzDuQ91gy9NbU8T7mv6ledBB/pdmNfn33Daa/Z4r2VpoMn95L2a+sGfLGWvQUbv9lUE7vFR2OU0tp4N+dF+rW07jLtq37Hq0pJTg9IqkLwHXp/0xwIoKXLsB2DVnvzfZEvD58jRI2hTYFljZwrn50huBD0jaNLWecvMXukZV9PhQd06/exhU8F9BTJpUubKsMP+eN1hjI0yZUlpM35Lsj8v1TGIVMH8hzPxz5et01FEwqIOtSLflljBhwp706FG5wDppxO4VK6saSglOXwWmABeTtSj+nNLKNQfoL6kf8BLZAIcvNMszHRhL9q7VKODeiAhJ04HrJF1ENiCiP/AIWStovTLTOfelMqalMm8tdo0K3J/ZRmlDAs78+TBzZsv5mjTv1at0IMn+yEOPHpUr06qjlGXaXwCOqfSF0/OdCWTD1euAqyLiKUnnAnMjYjpwJfCrNOBhJVmwIeW7kWzwxGrg5IhYA5CvzHTJM8lGH/43MD+VTaFrVFNjI0ydCuPG+X8Sqz0tBZ8NDTilBpjmvacOJJ2b3Ehonfr6+pg7d26rzr3wQjjjDLjgAji9zMdOZq1RLACVEnw2JOA4wFgTSfMior6UvKV061mFjRu37k+zaikUhFoKQMWCjwOOtQUHp3bQo4dbTFZ5+QJRsSBUKAA5+FgtaDE4STqF7L2mN4BfAoOAsyJiRpXrtlHzcycrx4YEonxByAHIal1Jo/Ui4idp/rodgHFkwcrBqQxTp2bPncCtKCuunEDkIGQdVSnBqekl1aOBqRHxWJpCyMrg505WSPNg5EBknVEpwWmepBlAP+BsSVsD71W3Whs/P3eyXLkBKV8wciCyzqaU4HQSMBB4LiJWSdqerGvPzFqppdZRbjByILLOqJTgdCiwICL+maYxOgD4SXWr1Tl4UETnUmrryMHIrLTgdBmwv6T9gTPIZlS4Fji8mhXrDDwoYuNXLCC5dWRWWCnBaXWam24k8JOIuFLS2GpXrDPwoIiNUykBycHIrLhSgtMbks4mW8dpSFrBtnJz2HdiHhSx8XBAMqusUoLTCWSzhZ8UEX+X1Ae4sLrV6jz83KljawpKDz7ogGRWSaXMSv534KKc/RfInjlZBfi5U8dTqJXkgGRWOaVMX3QI8FNgL2AzsqUo3oyIbatct07Bz506jmKtpMMOc0Ayq6RSuvWmkK1x9FugHjiRbHE/qwA/d6ptTd2uxxwDp532/kJ4biWZVVdJs5JHxGJJdWlBv6mSqrB4cufl5061JV+33axZWWByK8msbZQSnFZJ2gxYIOkCYBnQrbrV6lz83Kk2FOq2O/pouOgiGDrU/4AwayulBKcvkz1nmgCcCuwKfK6aleps/NypfeULSvm67fwPB7O2U8povb+lzbeAcypx0TQ/3w1AX2AJ8PmIeDVPvrHAf6Xd/46Ia1L6gcDVwBbAHcAp6UXhvOVK+iJwZirnTeAbEfFYKmsJ2VpVa8heOC5pCeFK6tEjC0zu2mtbhYKSu+3M2l/B4CTp8WInRsR+ZVz3LOCeiJgs6ay0f2ZuhhRoJpINwgiy2dGnpyB2GTAeeJgsOI0A7ixS7vPA4SlQfRK4Ajg453JHRERjGfdTNnfttR0HJbPaV6zl9B5ZULgO+ANZy6lSRgJD0/Y1wCyaBSdgOHB3RKwEkHQ3MELSLGCbiHgopV8LHEsWnPKWGxG5AzgeBnpX8F4qwl171eegZNZxFAxOETFQ0oeBMWQB6un0c0ZErC7zujtFxLJ0nWWSdsyTpxfwYs5+Q0rrlbabp5da7klkgaxJADMkBfDziLiiUKUljSdrsdGnT58it7fh3LVXXY2NMHbsukPBHZTMalfRZ04R8SxZ19pESSeQzQxxPiVMXyRpJrBznkPfLbFu+VbbjSLpLRcoHUEWnD6ak3xYRCxNgexuSc9GxOx856fAdQVAfX19SdfcEO7aq7zmrSUHJbOOoWhwktSL7AXc44BXyUbr3VxKwRFxVJFyX5bUM7VuegLL82Rr4P0uOsi64mal9N7N0pem7YLlStoP+CXwyYhYkVPPpennckk3A4OBvMGp2ty1Vzn5uvCOPhquucZByawj2KTQAUl/InvW1AX4CjAWuB3YLA1WKMf0VB7p56158twFDJO0naTtgGHAXanb7g1Jh0gS2YwVTefnLTdNVvt74MsR8dece+yWlp1HUrd0jSfLvLdWy+3aa2zX4RkdV2MjTJoEY8bAOee831qaONGByawjKdZy2o2su+xrpOcsiVL67mVcdzJwo6STgBeA4wEk1QNfj4h/i4iVks4D5qRzzm0aHAF8g/eHkt/J+8+Q8pYLfB/oDlyaxbO1Q8Z3Am5OaZsC10XEH8u4r7K5a6/1/FzJbOOhiIo/OukU6uvrY+7cuRUv11MZbTg/VzLrGCTNK/Vd0pLm1rO241F7pfNzJbONl4NTDXLXXsvchWe2cXNwqkEetVeYu/DMOoeSgpOkOrLBA2vzpxVxrQrctZdf89aSu/DMNl6lrIT7H2Qv4r5MNqURZKP1yplbz1rgrr33ubVk1vmU0nI6BfhQ7ourVn3jxsE//5l9Ghs77x9ht5bMOqdSgtOLwGvVroitq0cP6NYtaz1169b5Wk9uLZl1bqUEp+eAWZJuB95pSoyIi6pWKwM678AIt5bMrJTg9EL6bJY+1kY628AIt5bMrEkpK+FWZPVba53OMjDCrSUzy1VsJdz/jYhvSfoDeZakiIhjqlozAzb+gRFuLZlZPsVaTr9KP/+nLSpi+W3MAyPcWjKzQoqthDsv/fxT21XH8tnYWk9uLZlZSzx9UQewMbWe3Foys1I4OHUQHb2OX8fAAAAQAklEQVT15NaSmW2IlpZprwMmR0QH/rf6xiG39TRnTsdqbbi1ZGYbqmhwiog1kg6UpPCqhO1u3DiYNSv7Iz91au1377m1ZGatVUq33nzgVkm/Bf7ZlBgRv69arSyvHj2yFseUKbXfvefWkpmVY5MS8mwPrACOBD6TPp8u56KStpd0t6RF6ed2BfKNTXkWSRqbk36gpCckLZZ0iSQVK1fSUEmvSVqQPt/PKWuEpIWprLPKua+20NS9d8452R//xsb2rtG6Ghth0iQYMyYLTEcdBRMnOjCZ2YZRe/TWSboAWBkRk1NA2C4izmyWZ3tgLlBP9hLwPODAiHhV0iNks6U/DNwBXBIRdxYqV9JQ4DsR8elm16gD/gp8AmgA5gBjIuLplu6hvr4+5s6dW86vodVyWyUTJ2bBqr2nN/KS6WbWEknzIqK+lLwttpwk9ZZ0s6Tlkl6W9DtJvcus40jgmrR9DXBsnjzDgbsjYmVEvArcDYyQ1BPYJiIeSs/Brs05v5Rycw0GFkfEcxHxLjAtlVHTmrr3Lrgg2z/jjOwZVHtpCpbnnPP+syW3lsysHKU8c5oKXAccn/a/lNI+UcZ1d4qIZQARsUzSjnny9CJbrqNJQ0rrlbabp7dU7qGSHgOWkrWinipwjYMLVVrSeGA8QJ8+fVq8yWrq0SMbENHUrdcez6A84MHMqqWU4LRDROT+u/xqSd9q6SRJM4Gd8xz6bol1U560KJJezKPAbhHxpqSjgVuA/htaVkRcAVwBWbdeC9dsE+0xxNxdeGZWbaUEp0ZJXwKuT/tjyAZIFBURRxU6lroHe6bWTU9geZ5sDcDQnP3ewKyU3rtZ+tK0nbfciHg9p153SLpUUo9U1q4FyuowcoeYjxlTvdZLvqDk1pKZVUMpwemrwBTgYrJWxZ9TWjmmA2OByennrXny3AX8MGck3zDg7IhYKekNSYcAfwFOBH5arFxJOwMvR0RIGkz2rG0F8A+gv6R+wEvAaOALZd5bm2t6BtU0SGLmzMq2ohyUzKytlTJDxOeqsDzGZOBGSSeRLWR4fLpePfD1iPi3FITOIxtBB3BuRKxM298Arga2AO5Mn4LlAqOAb0haDbwFjE6DKVZLmkAWCOuAq9KzqA4n9x2oBx8svxXVFJBWrYL58x2UzKxttTiUXNKsiBjaNtXpONpzKHlLmr8Ae9RRMGgQbLll4aCSG4xg3YDUVIaDkpmVY0OGkpfSrfegpCnADaw7Q8SjrayfVVnzVtTMme8HmgcfzAJVc82DEZQW1MzMqqGUltN9eZIjIo6sTpU6hlpuOeUq1D2XT1MwAgckM6u8irWcJG0CXBYRN1akZtbmevTIphOC9bvucjkYmVktaWlW8vfSgAEHp41AbqAyM6tlpUz8erek70jaNU2sun2a987MzKwqSn3PCeDknLQAdq98dczMzEoIThHRry0qYmZm1qRgt56kM3K2j2927IfVrJSZmXVuxZ45jc7ZPrvZsRFVqIuZmRlQPDipwHa+fTMzs4opFpyiwHa+fTMzs4opNiBif0mvk7WStkjbpP2uVa+ZmZl1WgWDU0TUtWVFzMzMmpTyEq6ZmVmbcnAyM7Oa4+BkZmY1x8HJzMxqjoOTmZnVnHYJTmlm87slLUo/tyuQb2zKs0jS2Jz0AyU9IWmxpEskqVi5kk6XtCB9npS0pmlmdUlLUlkLJNX+6oFmZp1Ae7WczgLuiYj+wD1pfx0peEwEDgYGAxNzgthlwHigf/o0TaeUt9yIuDAiBkbEQLKpmP4UEStzLndEOl7SCo1mZlZd7RWcRgLXpO1rgGPz5BkO3B0RKyPiVeBuYISknsA2EfFQZGvMX5tzfinljgGur8xtmJlZNbRXcNopIpYBpJ875snTC3gxZ78hpfVK283TWyxX0pZkrazf5SQHMEPSPEnjW31HZmZWMaUsNtgqkmYCO+c59N1Si8iTFkXSS/EZ4MFmXXqHRcRSSTuSrfr7bETMzluhLHiNB+jTp0+JlzQzsw1VteAUEUcVOibpZUk9I2JZ6qZbnidbAzA0Z783MCul926WvjRtt1TuaJp16UXE0vRzuaSbyZ5v5Q1OEXEFcAVAfX29J781M6uS9urWmw40jb4bC9yaJ89dwDBJ26WBEMOAu1J33RuSDkmj9E7MOb9guZK2BQ5vltZN0tZN2+kaT1bmFs3MrLXaKzhNBj4haRHwibSPpHpJvwRIXW/nAXPS59yc7rhvAL8EFgP/B9xZrNzkOGBGRPwzJ20n4AFJjwGPALdHxB8rfbNmZrZhlA14sw1VX18fc+f6tSgzs1JJmlfqKzueIcLMzGqOg5OZmdUcByczM6s5Dk5mZlZzHJzMzKzmODiZmVnNcXAyM7Oa4+BkZmY1x8HJzMxqjoOTmZnVHAcnMzOrOQ5OZmZWcxyczMys5jg4mZlZzXFwMjOzmuPgZGZmNcfByczMao6Dk5mZ1RwHJzMzqzntEpwkbS/pbkmL0s/tCuQbm/IskjQ2J/1ASU9IWizpEklK6cdLekrSe5Lqm5V1dsq/UNLwnPQRKW2xpLOqdc9mZla69mo5nQXcExH9gXvS/jokbQ9MBA4GBgMTc4LYZcB4oH/6jEjpTwKfBWY3K2sAMBrYO+W9VFKdpDrgZ8AngQHAmJTXzMzaUXsFp5HANWn7GuDYPHmGA3dHxMqIeBW4GxghqSewTUQ8FBEBXNt0fkQ8ExELC1xvWkS8ExHPA4vJAt5gYHFEPBcR7wLTUl4zM2tHm7bTdXeKiGUAEbFM0o558vQCXszZb0hpvdJ28/RiegEPFzin+TUOLlSIpPFkLTaANyXlC4Sl6AE0tvLcWrOx3MvGch/ge6lFG8t9QHn3slupGasWnCTNBHbOc+i7pRaRJy2KpLemrHwtx4JlRcQVwBUtXKtFkuZGRH3LOWvfxnIvG8t9gO+lFm0s9wFtdy9VC04RcVShY5JeltQztZp6AsvzZGsAhubs9wZmpfTezdKXtlCdBmDXAucUSjczs3bSXs+cpgNNo+/GArfmyXMXMEzSdmkgxDDgrtQd+IakQ9IovRMLnN/8eqMlbS6pH9kgikeAOUB/Sf0kbUY2aGJ6uTdnZmblaa/gNBn4hKRFwCfSPpLqJf0SICJWAueRBZA5wLkpDeAbwC/JBjb8H3BnOv84SQ3AocDtku5KZT0F3Ag8DfwRODki1kTEamACWSB8Brgx5a22srsGa8jGci8by32A76UWbSz3AW10L8oGvJmZmdUOzxBhZmY1x8HJzMxqjoNTG9tYpkuStCRNIbVA0tz2rs+GkHSVpOWSnsxJK2lKrVpT4F4mSXopfTcLJB3dnnUshaRdJd0n6Zk0BdkpKb1DfS9F7qMjfiddJT0i6bF0L+ek9H6S/pK+kxvSYLLKX9/PnNpOmi7pr2SDQBrIBnqMiYin27VirSBpCVAfER3uxUJJQ4A3gWsjYp+UdgGwMiImp380bBcRZ7ZnPUtR4F4mAW9GxP+0Z902RHqlpGdEPCppa2Ae2cwvX6EDfS9F7uPzdLzvREC3iHhTUhfgAeAU4DTg9xExTdLlwGMRcVmlr++WU9vydEk1ICJmAyubJZcypVbNKXAvHU5ELIuIR9P2G2SjZ3vRwb6XIvfR4UTmzbTbJX0COBK4KaVX7TtxcGpbhaZk6ogCmCFpXprWqaNbZ0otIN+UWh3JBEmPp26/mu4Ka05SX2AQ8Bc68PfS7D6gA34naYLsBWQTJdxN9urOP9JrOFDFv2EOTm2rNVMv1arDIuIAshndT07dS1YbLgP2AAYCy4Aft291SidpK+B3wLci4vX2rk9r5bmPDvmdpPdBB5LNnjMY2Ctftmpc28GpbRWbRqlDiYil6edy4Gay/3A7spfT84Km5wb5ptTqECLi5fRH5T3gF3SQ7yY91/gd8JuI+H1K7nDfS7776KjfSZOI+AfZ9HGHAB+Q1DT1XdX+hjk4ta2NYrokSd3Sw14kdSObWurJ4mfVvFKm1OoQmv6YJ8fRAb6b9PD9SuCZiLgo51CH+l4K3UcH/U52kPSBtL0FcBTZM7T7gFEpW9W+E4/Wa2NpCOn/AnXAVRHxg3au0gaTtDtZawmyyYOv60j3Iel6skmFewAvky1qeQvZFFd9gBeA43Omy6pZBe5lKFn3UQBLgK81PbepVZI+CtwPPAG8l5L/H9nzmg7zvRS5jzF0vO9kP7IBD3VkDZkbI+Lc9P//NGB7YD7wpYh4p+LXd3AyM7Na4249MzOrOQ5OZmZWcxyczMys5jg4mZlZzXFwMjOzmuPgZFYiSd1zZpX+e7NZpv/chvXomzsLeVuS9GaB9C0k/SlNblzo3JkdZdoea3+btpzFzAAiYgXZuyodcubvKvsq2UzVa4rk+RXwTaDDvBNn7cctJ7MKaGpRSBqaWhA3SvqrpMmSvpjWxXlC0h4p3w6SfidpTvoclqfMOkkXpuOPS/panjx9Jd0v6dH0+UhOPWZLulnS05Iul7RJKvNqSU+m+pya8u8h6Y9pIt/7JX04pfeT9FCqw3lFfgVfJM0UIKlnuvaCdJ2PpTzTyV5GNWuRW05mlbc/2QSZK4HngF9GxGBlC8/9B/At4CfAxRHxgKQ+wF2sP6nmScBrEXGQpM2BByXNYN2JNpcDn4iItyX1B64H6tOxwcAA4G/AH4HPAs8DvXLWfvpAynsF8PWIWCTpYOBSsqURfgJcFhHXSjo5382mqbh2j4glKekLwF0R8YPUzbclQES8KmlzSd1TK9SsIAcns8qb0zQ1jaT/A2ak9CeAI9L2UcCAbCo2ALaRtHVaA6jJMGA/SU3zmG0L9CdbsLJJF2CKpIHAGmDPnGOPRMRzqR7XAx8F7gF2l/RT4HayZU+2Aj4C/DanPpunn4cBn0vbvwLOz3O/PYB/5N4/cFWaAPWWiFiQc2w5sAvg4GRFOTiZVV7uPGPv5ey/x/v/z20CHBoRbxUpR8B/RMRd6yRm6wQ1OZVsTr39U5lv5xxrPjdZpNbL/sBw4GSyFVq/RbZGz8AC9WhpjrO3gK45F5mdllD5FPArSRdGxLXpcNeU36woP3Myax8zgAlNO6nl09xdwDdSCwRJe6ZZ4HNtCyxLSzF8mWySziaD0zOjTYATgAck9QA2iYjfAd8DDkjrDT0v6fh0HaUABvAg2ez5kD1XWk9EvArUSeqazt8NWB4RvyCbofuApnKBnckmPjUrysHJrH38J1CfBjo8DXw9T55fAk8Dj6ah4z9n/d6OS4Gxkh4m69L7Z86xh4DJZMszPE82k3wvYJay1U2vBs5Oeb8InCTpMeApsuXRAU4hW0xyDlkgLGQGWbchZLOiL5A0n6xL8Ccp/UDg4ZxVVM0K8qzkZhshSUOB70TEp9voeoOA0yLiy0Xy/ASYHhH3tEWdrGNzy8nMyhYR84H7ir2ECzzpwGSlcsvJzMxqjltOZmZWcxyczMys5jg4mZlZzXFwMjOzmuPgZGZmNef/A66tDVX5EZ7GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x205a5305b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.integrate \n",
    "\n",
    "def ODE(m,t):\n",
    "    dmdt = -0.3*m\n",
    "    return dmdt\n",
    "\n",
    "massODEINT = scipy.integrate.odeint(ODE, mass_initial, time, rtol=1e-5)\n",
    "\n",
    "# Put the code to plot the errors of the two solving methods here...\n",
    "### BEGIN SOLUTION\n",
    "plt.ylabel('Error in Mass of Sample (kg)')\n",
    "plt.xlabel('Time elapsed (s)')\n",
    "\n",
    "plt.scatter(time, massODEINT[:,0]-true_masses, s=1, c='r', label='ODE solver')\n",
    "plt.scatter(time, massEuler-true_masses, s=1, c='b', label='Euler Method')\n",
    "plt.ylim([-1e-3,1e-3])\n",
    "plt.legend()\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f7449647f1fa4f70",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## 5) $F = ma$ as an example of a second-order differential equation <a class=\"anchor\" id=\"2nd-order-ode\"></a>\n",
    "\n",
    "Here we will study the second-order ODE we get from a mass $m$ attached to a spring\n",
    "of spring constant $k$\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}^2x}{{\\rm d}t^2}=-\\frac{k}{m}x\n",
    "$$\n",
    "\n",
    "where we want to solve for the position as a function of time, $x(t)$. \n",
    "This is done by splitting it into two first-order ODEs, using\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}x}{{\\rm d}t}=v\n",
    "$$\n",
    "\n",
    "for $v(t)$ the velocity of the particle. Then we have two coupled\n",
    "first-order ODEs\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}x}{{\\rm d}t}=v~~~~~~~\\frac{{\\rm d}v}{{\\rm d}t}=-\\frac{k}{m}x\n",
    "$$\n",
    "\n",
    "which we solve simultaneously for $x(t)$ and $v(t)$.\n",
    "To do this we need the values of $k$ and $m$, as examples\n",
    "we take $k=5$N/m and $m=1$kg. Then the two ODEs are\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}x}{{\\rm d}t}=v~~~~~~~\\frac{{\\rm d}v}{{\\rm d}t}=-5x\n",
    "$$\n",
    "\n",
    "As we have two first-order ODEs we also need two BCs. Here these\n",
    "are the initial positions and velocities, e.g.,\n",
    "$v(t=0)=0$m/s and $x(t=0)=0.1$m.\n",
    "\n",
    "We start with the Euler method for integrating two coupled ODEs simultaneously. The\n",
    "code for one Euler step is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-cf3e9d626dee4065",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# initialising arrays\n",
    "t=np.empty(2)\n",
    "x=np.empty(2)\n",
    "v=np.empty(2)\n",
    "\n",
    "t[0] = 0\n",
    "x[0] = 0.1\n",
    "v[0] = 0\n",
    "h = 0.01\n",
    "i = 0\n",
    "\n",
    "dxdt = v[i]\n",
    "dvdt = -5*x[i]\n",
    "t[i+1] = t[i] + h\n",
    "x[i+1] = x[i] + h*dxdt\n",
    "v[i+1] = v[i] + h*dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ad5a7aadcb2c2f9c",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Iterating this integrates the particle's position and velocity in time. Note here\n",
    "${\\tt x}$ and ${\\tt v}$ are arrays but ${\\tt dxdt}$ and ${\\tt dvdt}$ are just variables.\n",
    "\n",
    "Like the previous example this too can be done using the scipy.integrate package. However there is an added complexity due to it being second order so the function we need to define has to contain both parts of the solution we are interested in. This means you will need the function to return a list of values rather than just the one value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c8e68ca4225cc14c",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### Task 3 <a class=\"anchor\" id=\"Task-3\"></a>\n",
    "\n",
    "Solve the two coupled ODEs\n",
    "\n",
    "$$\n",
    "\\frac{{\\rm d}x}{{\\rm d}t}=v~~~~~~~\\frac{{\\rm d}v}{{\\rm d}t}=-5x\n",
    "$$\n",
    "\n",
    "with the two BCs that\n",
    "are the initial positions and velocities:\n",
    "$v(t=0)=0$m/s and $x(t=0)=0.1$m. Solve from $t=0$\n",
    "to $t=10$s. Use 1000 steps.\n",
    "\n",
    "First use the Euler method, then after you have checked that works,\n",
    "do it using the scipy.integrate package. Make sure you understand the ODE function as it is more complicated now it is second order. \n",
    "\n",
    "If you compare the two solutions at the same step size\n",
    "you will see that the solution using the scipy.optimise method is a lot more accurate.\n",
    "The exact solution of these ODEs with these BCs is known, it is\n",
    "\n",
    "$$\n",
    "x(t)=0.1\\cos(\\sqrt{5}t)\n",
    "$$\n",
    "\n",
    "Show that the solution is $x(t)=0.1\\cos(\\sqrt{5}t)$,\n",
    "by plotting the array you obtain for $x$ by numerically solving\n",
    "these ODEs, as well as the $x$ values given by this function.\n",
    "\n",
    "Plot both your calculated $x$  values and the exact $x$ values\n",
    "on the same plot. If your program is working correctly and you have\n",
    "chosen a small enough step length the two curves will be almost on top\n",
    "of each other. You can quantify this by in addition plotting the\n",
    "difference between the exact solution and your numerical solution.\n",
    "How big is the largest difference between your numerical solution, and the\n",
    "exact solution, for 1000 steps and $h=0.01$s? How big is the largest\n",
    "difference when you reduce the number of steps to 100 (and so $h=0.1$s)?\n",
    "\n",
    "What this task should demonstrate to you is building a model of an oscillator. Since there is no driving force it should show constant amplitude oscillations. In the assignment you will see oscillation with a driving force so variable amplitudes can be observed for instance if you make the driving force frequency equal to the natural frequency you should see that it diverges. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2db184c2c51643cc",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnWd4VNXWgN896TPpE0ihJKFLb1IUATs2ENRrwV6wXfWq116ujWv3w94LXlREBEXsIlEElN4CoaQRCGmTPpM2M/v7MZkQIQlTzplJOe/zzEOY2WedNWfO2WvvtddeS0gp0dDQ0NDQcAWdvxXQ0NDQ0Og4aEZDQ0NDQ8NlNKOhoaGhoeEymtHQ0NDQ0HAZzWhoaGhoaLiMZjQ0NDQ0NFxGMxoaGhoaGi6jGQ0NDQ0NDZfRjIaGhoaGhssE+lsBpYmLi5MpKSkeH282mzEYDMop1AHoat+5q31f0L5zV8Gb77xx48YSKWW3Y7XrdEYjJSWFDRs2eHx8WloaU6dOVU6hDkBX+85d7fuC9p27Ct58ZyFErivtNPeUhoaGhobLaEZDQ0NDQ8NlNKOhoaGhoeEymtHQ0NDQ0HAZzWhoaGhoaLiMZjQ0NDQ0NFxGMxoaGhoaGi7T6fZpaHjGvqIqftlVRIw+iPNGJKEP1m4Nb6m32vkhvYC8Ugsn9DUyqneMv1XqFGSXmPl5ZwHhIUGcOyKRyNAgf6vUpdB6Bg1+3d/AJz+twmZ31It/9dd9zL92HH27hftZs45LuaWeKz9Yx7YDFU3v3XBSKg+efRxCCD9q1rFZuvkA9y7eRoPNca++smIvC64fR7/uEX7WrOuguae6OL/sLOR/O+uZ3D+OdQ+dysI5E6ipt3H9/A1U1jb4W70Oic0umfPxRjIOVfHaZaPY+ugZXD6hN++uymb+mhx/q9dhWb2vhLsXbWVscix/PnAqX958Ala75PL31lFmrve3el0GzWh0YSprG3hg6XZ6Reh48/IxdI8IZUIfI29ePob9pRae/2G3v1XskMxfk8O6nFKeuWAY5w5PIkofxJMzhnLacd2Z+90uckrM/laxw2Gpt3Lv4m2kGA28d9VYEqJCGZMcw0fXHE9JdR1PfrvT3yp2GTSj0YV5My2Tkuo6rhkaTGhQQNP741JjmT2+N5+u28/ewio/atjxqKhp4P9+3sPUgd2YOapH0/tCCP47cxhBATqe+zHDjxp2TD5ak8PB8hqeuWA4hpDDXvWhPaK4cUoflmw6yI6DFW1I0FAKzWh0USosDfxvbS5nD0skNSrgqM//ddoAggN0vPVblh+067jMX5NDVZ2Ve84ceNTaRffIUK6flMp32wvYV6QZY1ex1Ft59/csTh7YjXGpsUd9PmdyXyJDA3llxV4/aNf10IxGF2XRhjyq66zcMrVvi5/HGoL5x9ieLNt6kKLKWh9r1zGx2uz8789cpg7sxpCkqBbbXHlCCsGBOj5cneNb5Tow32zNp8zSwC0n92vx86iwIK6cmMLPuwo5WF7jY+26HprR6IJIKVm0IY9RvaNb7dwArj4xlQabZMnmgz7UruOycncxxVV1zB6f3GqbuPAQZoxIYunmg1jqrT7UruPy6bo8+ncPZ2xy6yHLFx/fC4AvNuT5Sq0ui1+NhhBimhBitxBinxDi/hY+v0sIsVMIsU0IsUII0frTqOEyWw9UsLeomn+M7dVmu9Q4A6N7R7N000GklD7SruPy+fr9dIsI4eSBbdexuWBMTyz1Nn7eWegjzTouuw5VsjWvnEvG9W4zVLlXrJ5J/eL4YsOBptBxDXXwm9EQQgQArwNnAYOBS4UQg49othkYK6UcDiwGnvOtlp2TJZsOEBqk45zhicdsO3NUD3YXVrHrkOaDb4sycz0rdxcza1QPAgPafqzGpcSSFBXK11vyfaRdx+WbrfkE6MTfggpa46KxvThYXsPG3DIfaNZ18edMYxywT0qZJaWsBxYCM5o3kFKulFJaGv/7J9DTxzp2Oux2yU/phUwZ0M2lnbTnDk8iUCf4ZpvWwbXFiowibHbJ2cOObYh1OsF5I5P4bU+xtr/gGPyQXsCEPrHEGoKP2faUQd0JDtTxY3qBDzTruvjTaPQAmjsgDzS+1xrXAd+rqlEXYPvBCgoqazljcIJL7WMMwYxLjWXFLs2V0hY/pReQEBnK8J6trxE156yhidjskt/2FKusWcdlX1EVWcVmpg1x7V4NDwlkUr84fthRoLlTVcSfaURaclC2+EsLIS4HxgJTWvl8DjAHID4+nrS0NI+Vqq6u9ur49s6Xe+rRCQgx7SUtbR9w7O+cHNTAmsJ6vvjuV7rpO37shNK/cZ1NkpZh4aSegfz2228uHWOXkshgWPj7dqIr1A8V7Yj39bJMxywsvCKLtLQcl45JCWzg1/J6Pv7mV4y6mg73nb3FF7+zP43GAaD5SmxP4CgfiBDiNOAhYIqUsq4lQVLKd4B3AMaOHSu9KSbf2YvRP7Pld45PieDcMyY2vXes75xSYuazjDQs0X2YekKK+kqqjNK/8cqMIurt67nqtNFMHtD2InhzzijZyk/pBUw6afIx10G8pSPe129krGVoDyszp53k8jFDq+v4YMcvVEckkywOdLjv7C2++J39OWxcD/QXQqQKIYKBS4BlzRsIIUYBbwPTpZRFftCxU1FSXUdGQZVbHRtASpyBPt0M/Jqh/QQtsXpfCcGBuhY3nrXFKYO6U1lrZdP+cpU067iY66xs2l/GpH7u3atx4SEMToxk1V7N7acWfjMaUkor8E/gR2AXsEhKmS6EeEIIMb2x2fNAOPCFEGKLEGJZK+I0XGBtpgmAE/oa3T52Ur84NuSU0mCzK61Wh2d1pokxvWP+lorFFU7sF4cQsCazRCXNOi7rckqx2iWT+sW5feyk/nFsyi2nzqqta6iBXx3UUsrvpJQDpJR9pZRzG997VEq5rPHv06SU8VLKkY2v6W1L1GiLNZkmIkICGdbDtcXa5kzoY8Rcb2O7lt/nb5iq69h1qJJJ/d3v3KLCghicGMmfWSYVNOvYrN7rmL2NTXG/BsmkfnHU2+zsLrOpoJlGx1/V1HCZNZkljO8T65H/3Ol60Tq4v7M2y/PZGziM8ab95dQ2aB1cc/7YV8LYZPdnb+C4V4MDdew0abNiNdCMRhfhQJmFXJOFiX3dHxGDw1c8ID6cP7NKFdasY7N6n+ezN4DxqbHUW+1szdPWNZyYGtfeTvTANQUQGhTAiJ5R7NVmGqqgGY0uwl+Nnf3EPp6NiMExKtbWNf7O+pxSjk/1bPYGjlGxEGjGuBnOwAB3AwuaMyY5lpxKuzaDUwHNaHQRNu0vIzwkkIEJnpfFnNDHiEVb12iiwtLAvqJqRveO9lhGtD6Y4xK0dY3mbMwtIyhAeDx7AxibHINNos3gVEAzGl2EzfvLGdkrmgCd5/WpxzRmGd2shYgCsDnPkeNodG/3F2ubc3xKDFvyyrFqMzjAMcAZkhTl0XqGE+e9ukHLQ6U4mtHoApjrrGQUVDLKixExQHxkKElRoWzRRm+Aw42iEzCil3fXdWTvaGoabOwtqlZIs45Lg83OtgPlXhviGEMwiQahJS9UAc1odAG2HajALr0fEYOjg9uSpz2IAJv3lzEwIfJv5Uc9YWQvx++iGWNHKvTaBnvTTMEb+scEsDG3DLuWKl1RNKPRBdi039HJj/RyROyUkVdaQ0l1ixldugw2u2Tz/nKv1jOcpBj1RIUFsUVz+zXNDEYne39d+0brqKhpINtk9lqWxmE0o9EF2Ly/nD5xBmJcSC99LJpGxV28g9tXVE11nVWR2ZsQghG9orWZBo57NTEqlMSoMK9l9YlyrIns0AI3FEUzGl2ArQfKGanAiBhgWI8oAnSiy3dw2w44vv+IXp5H+DRnZK9o9hRVUV3XtUvA7siv8CpqqjlJBkFIoI5tBzSjoSSa0ejkFFXWUlxVx9A2aoG7Q1hwAIMSIrq80UjPr0QfHEBqXLgi8kb1ikZK2N6FO7jqOivZJWbFjEaATjA4KbJLX1M10IxGJyc9vxKAoQo9iOCYbaTnV3TpQjfp+RUMToz0KoS5Oc7iTV3ZlbLrUCVSKnuvDm+8V7W64cqhGY1OjrMTOi7R8019RzI4KZIySwMFlbWKyexI2O2S9PxKhiRFKibTGB5CQmQoOw9VKiazo+G8V4f0UO66Du0RhbneRnaJFs6sFJrR6OSk51eSGmcgwoV64K7i7CzTD3bNDi7bZMZSb2OIgiNicBjjnfld85oC7DhYSfeIELpHhComc3hPx1qetq6hHJrR6OTsyK9gsIIjYoCBCZEIQZcdFTe5/BRaJ3IyODGSfcXVXTZfUnp+haKuKYC+3QyEBQVoRkNBNKPRiamwNHCgrEZRNwpAeEggKUYD6fld80FMP1hBcICO/vHKLII7GZwUic0u2VvY9VwptY074ocqfK8GBugYmBBBRkHXHOCogWY0OjHOTl3pETE0ulK66ExjR34FAxMiCFK4rvfgREeHufNQ1zPGGQVV2OxScZcfONbzMgqqunTghpJoRqMT43SjKD3TAEcHl1daQ0VNg+Ky2zNSSnYcrGSogou1TnrH6gkPCeyS6xpNi+Aq3KuDEiIptzRQWNm1sxgohWY0OjHp+RUkRIZiDA9RXLbz4d7VxWYbhypqqahpaJoVKIlOJzguMaLJ2HclMgoqiQwNpEe09zvBj2RQYzmAXZqLShE0o9GJySioUjTUtjnOxfWuNireXVgFOIIB1GBwYiS7DlV2uSR7ewqqGZgQgRDK7HtpzqDG3yrjUJXisrsimtHopFhtdrKKzQyIV8dodI8IJS48pMuta+wpcHQ8AxReBHcyOCkSc72N/aUWVeS3R6SU7C6sUu1ejdIHkRQVym5tpqEImtHopOSYLNTb7Ko9iODoOPcWdq3R2+7CKuIjQ4jWe5/8sSWcM5jdXei6FlXVUVHT4FVVyWMxKDGSjIKuc03VRDManZQ9TW4UNY1GBHuLqruUK2WPiiNigP7dHTOYrmSMdzfN3lQ0GgkR7Cuqpt6qVUf0Fs1odFL2FFYhBPTtpo4bBaB/fDiWehsHy2tUO0d7wrmHYqCKnZshxLEYvKcL7dVwDnDUNBrHJUZitUsyi7vOdVULzWh0UvYUVpEcqycs2PM6y8fC+ZDvLeoao+L9pRbqrHYGqDh7A4fbb08Xm2l0iwghVoF6L63hDAjpatF+aqAZjU7K7gJ13SgAA7o3Go0uMip2ulHUnGmAwxhnFZux2rqGK2VPYZXq1zTFaCAoQGh12BVAMxqdkDqrjRyTRXWjEaUPontESJdxpThH/0qnDzmSft3DqbfZu0QEld0u2VNYrfq9GhigI8Vo6DIDHDXRjEYnJKvYjM0uVXejgHMxvGu4UvYUVtE7Vo8+OFDV8zg70K5gjPPKLNQ02BiYoK4hBoex19Y0vEczGp2QpsgplUdv4HgQ9xZ2jQgqtSOnnPTrQhFUvoicctKvWzi5JjN11q6ZRVgpNKPRCdlTWEWgTpAaZ1D9XAPiI6hp6PwRVPVWx2ZJX4yIDSGB9IwJY08X8L8fdvn5wGjER2CXkF1iVv1cnRnNaHRCdhdUkxpnIDhQ/Z/XuTO6s0f7ZJeYsdqlT0bE0Oj26+TXFGBvUTU9osMID1HX5QeOmQbAvi5gjNVEMxqdkH1FVaov1jrp171r+N+dvnA19700p398eJeIoMosrqZvd99c0z7dDAjRdaL91EIzGp2MOqsjb1E/H3VuUWFBJESGdvpRcWbj6LRPN/VdfuAIZ6632ckxdd4IKiklWcVm+vromoYGBdA7Vs8+bTHcKzSj0cnYb7Jgl9DHR0YDHKPizv4gZhY73ChqR045ado42YmNcUFlLZZ6m89mb+BwUe3TZhpeoRmNTkZmsWORz1cjYoA+cQayis2dujJaVonZp9c0tfFcWZ140TazyPf3ar/4cMf6VCd3+6mJX42GEGKaEGK3EGKfEOL+Fj6fLITYJISwCiEu9IeOHQ2n792XM40+3cKprrNSXNU5K6NJKcksqvbpiDg8JJD4yBCyijux0Wi8V33lSnWeq6tsnFQLvxkNIUQA8DpwFjAYuFQIMfiIZvuBq4FPfatdxyWr2Ex8ZIhPolGcOEeKmZ20gyuqqsNcb/PpiBigT1w4WSWd15WSVVxNREgg3SKUryzZGs49MFoElef4c6YxDtgnpcySUtYDC4EZzRtIKXOklNsAbS7pIpnFvh0Rw+FZTWft4JyL4L6/rp3b7ZdZbKZP93BVqvW1RtPGSc1oeIw/jUYPIK/Z/w80vucXqqvhIuOvRJ99L2mTH6U4vchfqniMIxql2ucj4sTIUEKDdE0+aoDcXzP5bfSd7NKPITPkOP7scQG/vppOR+z/Mkt873t3nC+cipoGSs31ANitdtbc8Tmzzyxh4EC48bjfSZv4AIc2HPSpXkqRWVxNXx9sQG1ORGjQUW6/g2v3kzbuXnYajicrZBB/Jc5g7d2LsWu1N1rEdz6Mo2lpeOFRlyKEmAPMAYiPjyctLc1tGRUVgfTpV4ZtewAnrZqLeeg8Fl30ON1vGeWJSn6hok5SWWtFlheQlmZy+bjq6mqPrllzuoXCht37WRlexPrXTdzx5WUkAFsiJlKrTyS1YAMTb48m6eNSHn54J1FRVq/O5w3uft/fd9UREgAZm/5ktw9HxeZixzX68uc/SCqspuddczmhZj0fRS8gccTpGLN3MSnjeWqPf5VF5zxM939PaFWWEr+xktRaJYcqahHVRarp1dp3jglsYGtWPmlpZSxfnkjWvF0stM1jS/hEiiIT6Ve0mfiXvuOyb37hH88EERtbr4p+auCT31lK6ZcXMBH4sdn/HwAeaKXtR8CFrsgdM2aM9IaVK1fKzO8y5LbwCVKCTJv1slfyfMmfmSUy+b7lMm13kVvHrVy50utz3/LJRnnSs7/KO+6QUodVft7/IVmwOb/p84Z6u3z9dSnDgq3yo+g7ZP66PK/P6Snuft8r3v9LnvPK7+oo0wa5JWaZfN9y+eoH62SBLkFWECFX3TBf2hpsh9uszJSboqdKCXLlaU+1KkuJ31hJth8ol8n3LZffb88/dmMPae073//lNjny8R/lC3N2S5DyzDPsMm/dYT2sdVb55aNbZFiYlL16SZmzp041HZXGm98Z2CBd6GP96Z5aD/QXQqQKIYKBS4BlftSniT5nDWTAwTTW9LiIT5aE8skn/tbINZwL0b7aLNWcvnEG8kxmFr12gNv/FcBFu58ifmRi0+eBQYJbboHfP8xkZvkHVE4+h8oDHaMgjq8jp5z0iAkjSCfI/98n6KSNwqVrmfTOleiapYfpPbUPQw/+xB99rmDqLw/z4y1f+1xPT/BHlJ+TPnEGyiwNXP7xROaev55vvxP0PP7wvRoQHMCsx0fwxx9wpulTrENGUJZZ6nM92yt+MxpSSivwT+BHYBewSEqZLoR4QggxHUAIcbwQ4gBwEfC2ECLdV/qFRIYwZt/n7JkyhxtugPQd7d8Zn1VcTWiQjqSoMN+ffNU+pBDMHfUGL70ErXlxxl42gL3/XUzf2nR2nHADsp1nx62pt5FfUUOfON93bkhBSKWOnJgESj7+nv7nD2mxWZA+iPHb3+e9gc9zwbvT2LjRx3p6QGaxGZ2AZKPe9yf/KweAHwf+g/s+H01AK8UtR4+G219KoVdDJntOvLrd36u+wq/7NKSU30kpB0gp+0op5za+96iUclnj3+ullD2llAYppVFK2fJToxIhoYKFC+Hi4KXUjT+JBkuDL0/vNpnF1aTGhaPT+c7vDlCeXcbo958EIPbVG1s1GE7GPHAGf5zxJCfkLWLNre17GpddYkZK6Nvd97O3l16CooNx7B48geMuH9Nm2yB9ELPW/JuYhBCuudhCTVmtj7T0jMzianrH6gkJVK8ccUtU5Vcx5Ln7ANC9cBcBxyiHPOzGE1g76wXGF37Dqsvf9oWK7R5tR/gxSEiAm27RMdqymtX/eNnf6rSJr3ctO9k+7R5GmXYAcKDGNcN60jf3si3iBFLeeYCyovZrjJvcKD6eaRxcux/zQ3PpawylwmahwYUdzLGxsOD1Cr7OHMJfM5/xgZae4y+X38bpjzO6dCs6oNDuWnTU5C9uY1PMKYz47H6KdxSqq2AHQDMaLjD+vzP4K/48xn77GPl/5R37AD9QZ7WRV2rx+YOY/uE6TtrzPptH3OLWDuaA4ACCF3zISXIV/3kqSGUtPSer2IwQ+KQ2SXPyLryTexvmcs1pVqx2yYEy1+qVTJkeRX7yRCb89gy5v2aqrKVn2O2SbD8McLZvh7RNkaw+7iZ6G/Uup2gROkH0J28QRD2Lb/xZZS3bP5rRcJEei18hECuZVz7mb1VaJLcxUaEvF8GlhIfeTeE1/T2M+eoRt3cwD5o+gLNvSeGNN2DXlvaZgiSrpJqkqDDCjuHGUJJtb/zBhPwlrDv9YcZMTHbo4UZCyL5LXqCBIPKvPCozT7ugoLKWOqudVB/P3u6/H16JfpSRv71CapyBbDcyGPQ5ayDz7sjhljWXs3mzikp6we97itlerH4ou2Y0XKTnpBT+GnUzE/fMJ/u3/f5W5yicI3xfjoh/+AG+Xtsdnn2OiKQIj3YwP/6YZAmzKJt+lYqaek6OyeLzWYbtoUco1CUw7tN/NQ0C3MlBlTA6iY0n3cnEg4vZvWirWmp6TI7J8V1SfLgIvvXDTYjvlnPfvRKjEVLjwhvXq1y/V295rDvR0fDmXXtV1NRz3kzL5OtM9V29mtFwg8Ef38+MkB95+J3e/lblKHIbH8Rko286OGmXVF5zB/9IXMWcOY73jtzB7ArGOEHkuEFMyFvEvmU7VdLWc3JKzKTE+a5z2/zir4wqTyNj5oPo4/RE64OJNQS7naJl5Pw7qSCKXQ/+TyVNPSe3sUZIsg+Nse3ue/hQXMc/r3O4+VK7GahpsFFY6foMNzoa3j9rMe+kDWD7O2vVUtVjck1m4vXqd+ma0XCDbkPjGXL7qSxcCNnZ/tbm7+SWWog1BBMV5pv1gW2vr+Liwle4dcoOgoMd7/Vp7AScI0lXGfbhXVjQU3j7U0qr6RXllnoqahpI8ZEhBnh3YQRfhs1m/Hs3NL3XJ87gdjLI6NQYPr59AzMzn2f7dqW19I4ck5ngQB2JkaE+Od+O9/9idNmvpE+/H0M3xwAg1ehMPe+eMT5z3lmUiljMjz2nuJ7eUNtgI7+ilu569SMn2zQaQohQIcSFQoiXhRBfCCE+FkLcK4Twaehre+L22yQvyLvZe8nD/lblb+SazPSO9d2IuO6p5ygW3Tj+9aub3ktpNBrZJe6lnTYOjGPDmJsYn7uoXQUaOKvm+Wr2tmULvLnheDIfW0Bo9OEONSXO0DSTdIfZ/+lHWJjg1ZfaV3RabomF3rF6n4WGVz7+f1QQyZg3rm96z1mvJNvNeiWG7ga2nXAz4w59Tc7P7cdNldeY6r27P2caQojHgNU40n38BbwNLAKswDNCiJ+FEMNV17Cd0bOXYFxyIRPXvUxFbrm/1Wkip8TiMx/x3qU7GFf0Lekn30ZY7OGNhD2iwwjQCY86uH7z/olAsvvOt5RU1StySnzre197x0IGhuZyww1/fz/FqKewsg5LvXuLnLGxMG/SYp76qAfm7PaT1TXHZPbZNT2wOpdxeYvZcvwcIpIimt5PjAwlJFDn1mK4k8Fv/JMGgsi9a56SqnqFc4AT7+eZxnop5Rgp5d1Syk+llL9IKZdLKV+SUp4HzAaCVdewHWKcexcRVLP5lnf9rQrgCLfNr6jx2Yj40EOvUUMow9+65W/vBwfq6BEd5lFd656TUnhuwlKu2XUvlnZSHyfH5Ai37eWDGVzhlkNc9/uVvN3/eWJi/v6Z83f1pHDQ6XcOpTvF1L+yQgk1vUZKSa7J4rN7ddmrueSQQv9Xbvvb+zqdcERQeVAZsfvwBNb3u4xBOxZjym8fUX/OgZpfZxpSym/bOlBKWSSl3KC8Su2fQZeNZkvUFFJ/eqtdpE/OK61BSnyyYFtVBd9kDuaXYXcR29941OfJRn3TCN1dTnx2OrnlUe0m11euyUJSVBihQeqH2+66612CaSD5xTuO+swZveXJdU09axDru53NCVs/axcZDYqr6qhpsPlkplFXp+OhHyfz8D/2kjTh6OCV1DgD2R7MigFi3/ovA9nNx5/7roBUW+SYzESFBREe7Oc1DQAhxFghxNLGsqvbhBDbhRDbVNesnWOePYdkaxabX1rpb1V8Gjm1cCG8UH87cW/PbfHz1DgDOSbPCgeddBLcmbyEhAeu8VZNRcguMfskN5Kt3ka/399nY+xppJze/6jPezfq4MkMDkDedDPxspBNT7Y5DvQJvlwn2vpNPZbyOm68qeWONDXOwH6TxaN64YNPTWTQ+GjefZd2kZMq1+Q797Qrc5lPgA+BC4DzgHMb/+3SjJk7i3dD/sn/Vvb0typND6IvonzSn13O6CF1TGildEOy0UBVrZUyD0a1QsCMsQc4z/QRe77w//6CXJO5aXFfTba8uIKetv3UXn59i59HhgZhNAR7tFYEMPrBaRSIRMQH73ujpiIc3qOh/nWdOv8J/gqZzJQpLX+eGmdwa7f9kfz7vN38b9cYtr+9xgstlSHHZPaZy88Vo1EspVwmpcyWUuY6X6pr1s4JjQ4l/aZXeWPFQIqL/atLrslMRGggMXp1w213LdjIvMzzeHHIB60mJUxpGhV71sENe3o2dQRz6L8feqqmIlRYGiizNPhk9LZ90S4OiSTGPHF+q20cbj/PZhqBoYF8etLj3Fr8GAf9XOQv12QmUCdIilY33Dbnl32Mq15F+ZTz0bXSyznTmHhapnja9T3pz14qX/Tv2ma91c7Bspp2NdP4jxDiPSHEpUKIWc6X6pp1AG64AUY3/MlvD/3kVz1yTBZSjAbVay0XP/0eFsIY9dylrbZJ8cL/DhDb38imnjMYunUBdW5svFKa3FLfuPxKSmBO+h28cEsWoVGt+8dTGt1+nhI/Zzgb5Bjmz/dYhCLkmCz0itUTGKDugm3OI+9jJYBBT7eeacD52+Z66PYLjzew5bjLGJ25yK+1YQ6W12CXvgsNd+WXuwYYCUzD4ZZyuqi6PENkV+/kAAAgAElEQVSGwPuG2xm64D6/6pFrUt/3Xl9dz9Bdi9icfD5RydGttusZE4ZOeO5/Bwi+6VqM0sSmx7/xWIa3OKNq1E4hsvADCw0NcPWNbS+ophgNHKqopbbB5tF5evSo4Zox24h7/l6/+uB9ca/a6m0ct+4j1sacTsLopFbbGQ3BGIIDPDYaALF3XY2eGrY9tsRjGd7S5PLzUeYCV4zGCCnlWCnlVVLKaxpf16quWQeh5MzLGVSzhX1f+6w+1N9osNk5UFajuo9489M/ECtLCb56dpvtQgIDSIoO89j/DjDyntP5IfR8vlsVcezGKuHsSNTeMHny41NYaryOYcPabufsaD0Ju3Vy7egtzCl/nh3v/emxDG+QUpJbYlH9Xt3yyu/E2ws4dPoZbbYTQpBs9GzjpJMh144nN7Av+iULPJbhLbklvk0h5IrR+FMIMVh1TToogx+/GCsBHHjWP3GiB8tqsNllU4SNWpQu+oUSEcfIe9t+EKExgspD9xQ40qan3bGUZzafSUmJx2K8IsdkJjEqVNVw2+yf9jLEsoHoicd+vJwdrTfXdfij51NDKKWvfeqxDG8oNddTVWdVfabxevpUTtOvIWb20GO2TYnTezXTEDrBxnP/w2tls/22XpRjshAeEojR4Jttc64YjUnAFiHEbi3k9mi6DY1nc9zp9Fv3qV/2bPgiGqWyEmblvcy8KzYR5MJie7JR75V7CuCSSyDSamLFa7u8kuMpOT4It815+jPsCAY+evEx26Z46X8HiOwZydYe5zB4xyKsteqn0D4SX0T51dbCl0sEvf4xkaDIY9+rvWMN5JVZsHnhshv23BV8yDUsXOixCK9wphBSe03TiStGYxrQHzgDLeS2ReovnI3eVsmmL32fxdDprlAzcmLpUqitE5xzUy+X2qcYDVTUNFBucT3b7ZGMGAFrgk9mwMu3HLuxCuSqnBJd2iXJaz5lW9RkEo8/dth2lD6IGH2Qx5vRmrjsMrrJIrbO8/3+osP7idS7VzfO/YGnK2/hivNcS/GTYtTTYJPkl3sWdgvQvz+cOaqIslf8k1E412TxaSbmtnJPhQM0D7M9MuTW2aarM2LuP+gbdoj30/r6/Nw5JRbCggLoFqHeztTk+y5mXvRjre7NOKq90bNkcM0RAg6dcAEjyn/j0AbfzvsraxswmetV9RHv/nwLfep3U3lO65FoR+Kt/x1g5INnkyX6svF735ctzTFZ0AnoGaNeB6f78H0uEl8y+WzXuiZvI6icPJC6kKf2X+nz9P5Wm528Mt+lZYG2ZxpfCyFeFEJMFkI0aSSE6COEuE4I8SOOWUiXJzw2mDOnh/DlYonN6tvIFGc0ilpT05JdxZxUuJgRQ6yt7s04ktTGUY+3D2Lvey5Gh2TP3C+8kuMu+03qz94+WZPKzbq3GfrYhS4fk+LFXg0nodGhPHHFXu7dejl1Po5ozjWZSYoOIzhQnXDbygOVjDy4nJ3DLiYwNNClY5yzHmeItacMeuQi7AgOzFvslRx3OVRRS4NN+rSgVVu5p04FVgA3AulCiAohhAlYACQAV0kpfXuF2jHXjk/nj5KBbH/tN5+e15ExVL1Rxq5nlxGAncTbXO/cesboEcLzDX5O+pw9iIywkRh//swrOe7inCGptRvcbof5X0eTd9acFvN3tUZKnIH8ihrqrJ6F3Tq59DJBRYXkpy8qvJLjLs79RGqx7amvCaOW2Ftdn70lRIYSHKjzeoATPzKR7ZGT6LHGtwOcHB+mEHLSpsmXUn4npZwtpUyRUkZJKY1SyhOklHOllAW+UrIjcOLsFHpygIoPfGdHbXZJXmkNySr6M0O+X8qBgGQGXDTC5WNCgwJIigrz+kEEKJh6KUPN68hdfcBrWa7idAGpFW67Y9FOzs57i8vOq3LruBSjASkP107wlFNPhY2B44l92LfrRWrv0QhctpRDuh4Mvd5FPyqObLe9Y/Veu/0Ayk+7kP51O8j+PsNrWa7iyxRCTrTKfQph6G5ga4+zGZT+JbYG30RRHaqood5mV+2GqcqvYkTRz2QOn4lws2BOSpzeqzUNJ6lPXksqWSz+03c5vnJMFuIjQ9AHu+bicBfTvI95lds4+0z3ZgzODtdbF1VgIJj7DGdo7nKf7bovt9RTblGvCqLFAhuKerN59HXu36tG78JunfS//wJs6Eh/Z7XXslxlv8lMaJCO7iquaR6JZjQURM66kHh7ATve9s1N01RrWaXR28pvLbzNjUTOucTtY5VYtAVIHhNH9MhUlvhww21OiXouP2mXpG76kq3GU4hOaX1nfUs07dVQ4LqGXjaLKCrZNu9Xr2W5gtr36k8/wW22eYQ887jbxzruVYtHmZmbk3R8D6YfX8AjOdd5JccdckwWkmMNPquCCJrRUJRh959DLSGUv+cbF5XaezQWroznybhXGH7DeLePTTHqKbM0UKFADYcbT9jOnWsupHDLIa9luYKavve9X6WT0rAPy5nup2+LaawBr8SoePidp1JBJLWffOm1LFfIdYaGq7RO9PuC/cRESyZPdv/YZKOemgYbxVXez7pOu7QbW7bAvn1ei3IJX6RlORKXjIYQIkAIkSSE6O18qa1YRyQiKYLPhzzJ+wfOxO4DD1WuyUJwoI6ESOUzhtZX11Py9WpmnGcnwINN0UqOik89TXAhX7L7ua+9lnUsquuslFTXqbZOdPDVJdgRDLpvhkfHpxj1ilzTkMgQdiSfy3F7v/bJRj9nqgs11omstVYeXjKKz+NvJ8iDRM/JTfeq98b4gvPqWcZ55N7xkteyjoXdLhv3aPhuPQNcK8J0G1AI/Ax82/harrJeHZbAB+7hf6az+dMH6X1yShw7QdWYmm57eSU/WSYxp/cPHh3flO1WgQ6u34whZAf1x/CT+j6qXJVnb5Vbs9gWOYnuwxM8Oj7Z6F222+bU3XQH18n3WO0Db2qOyaJaWpbtb6wiVpYSMf1kj45PjnWGiHt/XXv3CybVUEzSSvXTChVW1VJntbfLmcYdwEAp5RAp5bDG13C1FeuonHsuDA7cw9aX01Q/l5rVuiyfLKWKcIb/6xSPjneOKL1dtAVHfp+c0bMYblpJeXaZ1/LawqmvGkZj3z44v+wjfnvoZ49lpBj1HCyroV6BlDXjbxvHT6EzWPyVOgv+zVHTjVI5fyk1hDLs32d6dHyPmDACdEIRtx9AyeRZHFeziYNr9ysirzXUvFfbwhWjkQf4NqC7AxMVBR9H3ca0pXNUTUEtpSS3VJ1qXXarnYG7v2ZHz7MI9bBYTmhQAIlRoV5vmnLS7YaZBGFlx7PqTnJzVEx1sXSJ436YfpHnkS7JRgN2CXll3ndwBgNcdVIWCR8/p3reNLXWiaRd0i/9K7YlnIGhu2fygwJ09IwJU2wGl3y7w/W476VlishrDbVDw1vDFaORBaQJIR4QQtzlfKmtWEfGfPr5pDbsJfNb9eK1i6rqqG2wqzLTSH//T+LtBcjzZ3olJ1mhUEaAwVcdz6qQU1m3Ud1Rca7JTLeIEAwhyp/nlLmnMr/7v0lN9VyGM8fQfoWu6+V91vBQ5X3snL9eEXkt0bROpILRyPh0Ez1seTSc49292jtW71Xa+eakThtIZvAgIlZ8pYi81sgxWQgKECRFh6l6niNxxWjsx7GeEQxENHtptEL/u6cDcOA19W6aHBVz6Bd/8A31BDHknrO9kpOiUNgtgC5Qx+fX/8LD6ZdiUebZbpGcEgupKlzTgk35jKlcSe9h7oXZHkmyggEGAMMeOJcGAil5W70oqsPrRMoPcBbtGMxM8RWD75/ulZwUo4HsErPXYbdOtk+9ncXlp1Gmojc112SmV6yeAB+G24ILRkNK+biU8nHgJeDFZv/XaIXEsT3YqR9L3Gr1on1yVdoJKiXcWvI4d09YQ1TvKK9kJRsNlFTXU1XrfdgtwKxZYKupY+XnRYrIa4kck1mV2iS7n3MMIHre7l2lZKMhmPCQQMVmcFHJ0Ww1nkLKlqWquVMP79FQ3hh/sTyM8ikziO0X65WcZKOeqlor5QqEiAMkPH4zT8v7+e47RcS1iNppWVrDleipoUKIzcAOHDmoNgohhqivWsem6ITzGWDexKH0UlXk55jMBOoESR6uObRGejpkZAUz9OqxXstyjiyV6uAmnyTZrTuOqCfvVkTekVjqrRRV1amSEj3856VkBw2g77nHeSXHUW1OmbBbJ5bTzyelYR+Z36iToVWtdaLsX7OZlf4El53qfcbepmy3Crmoxo2D1HgLOz5Ux+0npfTLHg1wzT31DnCXlDJZSpkM3A28q65aHZ+Ex28mnkK+XuXdCKg1ck0WesXqCQxQdn9mzkPv8hQPM2O696NOpdJOOwkMEhxInczQnOU0KDQibI7Tp630wmJFbjnDS9PIHXW+2ykuWiKlcQezUgy4ZwYWwtj2yXbFZDYnt8SiyjpR7ktf8gT/4aypntfCcHJ4gKOQO1UHH8bezYMrTlElVUtxdR2Welv7nGkABillU8UWKWUaoIimQohpjRUB9wkh7m/h8xAhxOeNn/8lhEhR4ry+YODEWIx9Y/haJQ9VjkqjjL4/v8l5EStJSPS+c2vKlaTgqDjoHzOJluWqZBNWK4Tx5x/tPMGjxNzievbVtkg26skrtWC1KRPxlDA6ibPHmZib6X66GFdwZGJW/l6N/X0pGWGj6DkpxWtZvWKVnRUD6C+ZTgTVqqRq2a9yWpa2cCl6SgjxiBAipfH1MOB1iTohRADwOnAWMBi4tIVa5NcBZVLKfsD/Ac96e15fIQTcNfZ37vnxVCoPupfN9Fg4pqbK+zPz/sjluJrNlE72LhLFiSEkkG4RIYqN3gCG3XU6FsKoWqB8kMH+xvBgpdc0Fv0Sy7vxjzDsipGKyEsxGrDaJfnltYrIAzj7gjA2bYL9ucqva+SalC8SVLStgKFVaymYcL4i8pwh4koOcIbdcQpVhGNZqPzI0R/ZbZ24YjSuBboBS4CljX9fo8C5xwH7pJRZUsp6YCFwZG6FGcD8xr8XA6cKXxXCVYATJwlOkb+y4wXPdlW3hslcT3WdVfFRRtZLjo449U5lHkRwpr1QcPQWp2db4jQGpH+l+N6CHJOFGH0QUWEe5KJohbrKOgK++YoLz7agU8iTqFThoObMPK2KvxhH1t2vKSYToKbeRkFlreIzjYznlqFD0uOfygxwQNkQcYDQqBC29zyLQbu/VvxezTWZCdAJesT4NtwWXIueKpNS3i6lHC2lHCWlvENKqUQgWQ8cGwedHGh8r8U2Ukorjk2GrletcQNTdR3//HQTO0q8K3DTnKFzTqBExCGXKjsqVqvWctSvS9kbMpTkU/spJlOpbLfNKbnxIWbZv2DjJmXHD46FRWVHbtteXslntTO5sneaYjIPp2hRroPrPzqC2BAzUSuUTdXiXCdS+rrm/nWIPUFD6Hf+UMVkJscqu1YEYD93BvH2AtLnb1BUbo7JQs+YMIIUXtN0hVZXpoQQ86SU/xJCfAMcNWeVUnoXGA0tPfFHnseVNggh5gBzAOLj40lLS3NbmTqbZPk2CyRLj45vjfqkUxmX+y0rflhBQKgyeXdWH3QsAhfuSyetYJfX8qqrq/lqySpqKxLRjRrPQQW/v72ynsLKBn78ZSUhgcp08vZRgazTnci8l/dzww3ue0qrq6tb/I13H7TQL1qn6O9f8/6nVBFO5dggxeRKKQnWwR9bMuhV69r3b+07N6ei/ymcu+MNvv/ft4T1UqaT31joSIZYkrOLtLI9isisrg7g2pxHuHDmldz4e+trW6585+bYKuspqW7g+19WEqbQvVozNY7Rb2/iuBUx3NDXdV2OxY7sGiKCxVHfz93v7BFSyhZfwJjGf6e09GrtOFdfwETgx2b/fwB44Ig2PwITG/8OBEoA0ZbcMWPGSE8ZP/cXOfuVHzw+viX+fPBrKUFufPZnxWS++GOGTL1/uaxtsCoib+XKlfL996UEKTduVERkE8u2HJTJ9y2Xuw5VKCr3xrEb5Ktxj3l07MqVK496r67BJlPvXy5f/Gm3l5odxtZgk4d0iXJNjwsVk+nkzP/7TV730XqX27f0nY9k5/82SAly1XUfeq7YEbz92z6ZfN9yWW6uV0zmwgUNEqT844+227nynZvz7bZ8mXzfcpl+UNl79ZRTpBw0SDl5drtdDv3PD/KRr7Yf9Zm737k5wAbpQt/dVo3wjY1/jpRS/tb8BSixorce6C+ESBVCBAOXAEcma1kGXNX494XAr41fThWSjXqKLMr6HofdeRq/6E5n1RrlsnvmmCwkRYcREqiczF8XFpGcDKNGKSYSaJYiXYHEhc25NGUt/yx5TLHSmgfKLNilsruWd368gQT7IezneZYGvS2UKlHanEGXjeZgQC+Cvl2qmMymdSK9cutEPR+4gh+Dz2OC61VdXaK3gtlum3PlpCzuzrierB/3KiKv3NJAVa3V5zmnnLjiELuqhfeu9vbE0rFG8U8cs4ldwCIpZboQ4gkhhNP19T5gFELsA+4CjgrLVZIUo4FCi7I2SR+n5/XpP/HCxpNRytzlmpStLFdfWse7PyfzWo+nUTrMoLfC8e9O+v/b0RHnvqLMepEaleUOfvwLVgIYfM85isl0khJnILfUgl3BXdxCJ1g1+SHeK5mJWaGfS+l1orrKOoblfUton0SP6ry0xeEQcWUHOKeeEcD1vM/+l5UxxmoXXzsWrRoNIcSljesZqUKIZc1eKwGTEieXUn4npRwgpewrpZzb+N6jUspljX/XSikvklL2k1KOk1JmKXHe1kiO01NRJzHXKVuUZsYMMB8oZeuvJYrIyy21KNq5lS/cSRi19PzHCYrJdBIVFkSsIVjxBzFpfC9HqpZVyhiNw7uWlXsQ7yx8gGsn7CKmT4xiMp0kG/XUW+0UVCoXdgsQ/8iNvGe9mh8UCvjLKVE2ff+2/1tBJFWEXapc1JSTiNAgjIbgptBrpeh5YjK7wkYR+4cyobdNKYRUKhR2LNqaaawBXgQyGv91vu4Gpqmvmu9JUXgHs5PzplaRTxLlj7/itaxySz3llgZFRxnG3//AJIwMvfFExWQ2xxHKqOyDCFB04kyGmv/i0IaDXsvKNVkwBAdgNAQroBns2QO7MgRjL+2viLwjUbIyYnNOOgkGxRSy5y3vN6TVWW3kV9QoaohrFn5FJREe13k5FslGveKuVIDC8TMYWrWWom0FXsvKMZkRAnrGtDOjIaXMlVKmSSknHrGmsanRtdTpSFbJlWJMiSAjajxJ670fFSvtRqmvrmdM4Qp29ptOYKg6aceVTnvhpOdtMykmjvULdnsty+lGUWobUPaD7/IONzDjHHUelWSF83o5CQyEd4z3c/OKC7xO1ZJXWoOUyo2IbfU2jtv9FTt6nU1IpOc1SdpCyczMzUm6dSY6JBnPeV9jI9dkISkqTJUqiK7Qlnvqj8Z/q4QQlc1eVUKISt+p6DuUrBV8JBUnz2RA7XZyV3hXcb7Jn6lQUr1tL68kmgpCLlZ+uu8k2agnv6KG2gbl9sAA9D1nEJP7F/DaTu9HnY5ay8qN3BJ/ms9JYRtI7quOIU6MCiM4QKf4TAMguDFVy7ZX0rySk6uwy2/tKisPyyepu/ZmReS1RLLRQH5FreL3av9Zw1gXOpmtm7yXq1YKIVdpa6YxqfHfCCllZLNXhJQy0ncq+o7wkEAig4UqI42+/3Z0ytkverd5yjmyVCpyYsHe8dwQ+D7D7z5dEXktkWI0IKUjQklJhE4wfWYAK3+VlJs8fxhtdklemYXescp0bsU7ChlatYZD45TbWX8kATpBr9gwclVwpQy/+3TM6Kn+xLuZsdKpLr5cHsL8kBsZe/cUReS1hHPgkKdQtlsnQidYfNtv3LX3ZsrLvZO1X4W0LO7gSmr0vkKIkMa/pwohbhdCeFdJph0TrxeqjN56npjsWLj93btiNzkmMwmRoYpMTW02+Oz7aPZNOtfjsq6uoJYrBeDiE/LYY+tD+iMLPZaRX15Dg00qtmCb8cJydEiSblY+1LY5KUaDKvdqWGwY25OmMWCnd6lack1mIkIDiVEg3FbaJeLj+Vw4uYgIFUvAqeltmDULbFY7Py30vFxCZW0DJnO9KgkgXcWVkNsvAZsQoh+OENhU4FNVtfIj3fU6VTo3gPVXvsr55k/Iyzt229ZQ0o2ydf4WLi56hVPH5SgirzVSVHwQR57TgzBdHbplno+Km1KiK/QgBv/wNQcCkhlw0QhF5LVGstHA/lKLYtXmmmObPpNEez47Fu7wWEaOyRHlp8Q6UcZnm3ip9Gpu7Pmt17LaIrVpX5HyxnjcONgUNJ5eT97gsQznzLJduqeaYW9c+J4JzJNS3gkkqquW/4g3CA5V1FJTr6xPE2DinRPIpB9LvQjXzilRbo9G1WvzeY57Of54FWtSAtH6ICJDA1Vx++kCdewZNINhB7+nptSzugpKxr1XVsKakgFsPf46RWpntEVKnB5LvY3iauXrNQx+8HxSA/bz6Y7hHstQcj9R4ZtLHXte7jtPEXmtEaUPIlofpMoMTqeDioHjGZH/PRYP3YpKr2l6gitGo0EIcSmOTX7LG99TbntnO6O73nFJlCoy35wBA+CG5B+JfOkxj453Tk2V8GdKu6TvtiVsjT+DkDh1IlGcCCFIiTOoMtMAMFw+k3DMbJ/3i0fH55osBAfqSIj03kX37bdwl+0Fop5/xGtZx0LpIlfNiekVTr+Te/GVhxO4BpudA2U1ihmNXuuXsD16MsaBcYrIa4tklaL9ACKumoWeGrY9/6NHxzcFFyi0/uYJrhiNa3DkiZorpcwWQqQCC9RVy3/E6x2jQzVGGgCX9V7NFblPUpzufp3r/U0Li95PTXct2EhP234azvWuZrWrqJHt1smw26ZSQRR1Cz3r4XJNZnrH6tEpMDNY/eEekhLsnKD8PsmjSG4MhlDDlQJwzcQMnts9nczl7ifFPFBWg80uFRkRZ32XQd/6XVSe5pt7NdWoJ1ulazrs1smUilisX3gWEJNjspAQGUpYsH/CbcG11Og7gX8D24UQQ4EDUspnVNfMTzhnGmp1cIn/vIAA7Ox61v14bSWnpkVvLfHJdN9JilHPgbIaGhSqNtec4PBgvjz+aV4vvBCbB15FR0Er7w1xdUE1z/88go+SHlSsdkZb9IgJI0AnVBsVnzw9gul8Q94r7vtTD7v8FBjgvL8aO4IB96oXjdYcR9htDXVW5V3UQWGBpPedzrDsb6ivrnf7+JwS/4bbgmvRU1OBvTiq7L0B7BFCTFZZL79hCBLE6INUc6UMuHA4+wP7EPa9+1FUSm3skxJK0/PZEnsqsf1VKU9yFMlGAza75GCZ9/WcWyLinpv5vPIsVq927zjZWAVRCZff9ue+J4xaul15lteyXCEoQEfPmDDVZsWJY3uwwzCObqvdn8Hllii3R+Ox/ddxzsh8Eo/v6bUsV0iJ0yOlY3OiGthvvZ1L+Iy0390fWeSoULHTXVzR+kXgDCnlFCnlZOBMHKVXOy1qulKETpA16gJGlKygIte9gO3sEjPdI0LQB3u3YWznTrig8iM2PK5uJEpzUpqSwalzXadNgxFBO9n58s9uHVdcVUdNg02R0ZtcvJhi0Z1hN0/yWparqOl/ByiZNJMhlvXkrzvg1nE5jWlZ4sK9S8uSkwMbNsCUixO8kuMOh9eK1LlXx980ilWGs1j8lXvPcXWdlZLqOpL9lHPKiStGI0hK2ZSnQUq5h068EA6NJUpV2DTlpNuNs8gmlVX/y3HrOKWiUb74tAGdDs6/UJ3dyi2h5qItQEQEvBNzH+d+fb1bewtymmZv3l1XS4mFEXnL2XXcLAJ86G92lNM1qxJ2C9D7Dsem1L3PuzfbyDGZSYnzPi1L5m3zWMr5XDzTfVeOpzjDbtVa1wgNhWunZtF/wX/cStWSq2CUnze4YjQ2CCHeb9zYN1UI8S6w8ZhHdWDU9GkCHHf1eE7vmcFbf7pXliRHgT0a0i65/PkRfNz7YRJ8N3gjLjwYfXCAajMNgIZZF9PTtp8d7651+ZjD0SjeXddtz/2AAQvhV1/olRx3STYaqKq1UuZlnqjW6HPWQH6KmMUfO9zbz5urkBulxy/z6WcoIHWgMokkXeFwiLh6A8fZw7dzT80TbHlxhcvHOAeyHcFo3AykA7cDdwA7gZvUVMrfqO3T1AUILrlUsPKHOkryXDuHuc5KcVWd1yPiXQs20q9hF72npHolx12EEKq7UoY/MoMaQil/w/W9p7kmCwE6QY+YMK/O/WrWOVwa+S3Db1MvxUVLqO32A1hz15c8svtyDrqYTNhqs5OnQPr+zG8zGFS7BdMZl3olx10Oh4ird01H3j+NchFN3UefuXzM4fT97dw9JaWsA14DHgceBV5vfK/TorZPE+Cq0/M5YEtg533zXWqv1Aa0olc+o54ghv3HN+GLzXG6UtQiIimCzb2mc1z6Fy5P+3NLLfSMCSMowPNwp8pKWPJtCMYrzlYtU3Br+OJenT0bQmQNP7ziWo3vg+U1WBUIt817fiF2BAMfvsgrOZ6gVooWJyGRIWzvN4sRWUtc3pSaazLTLSIEQ4hv77EjcSV66hwgE3gZh/HYJ4TwTXiIn1Az7YWTwacmYgpJIurbT1xqr0TklN1qZ+Dmz9mccBbRqcoXBjoWyUYDeaUWbApWmzuSgMsvI0JWsv6D7S61d+7R8IYNj3zNQ7UPc/mFyhZEcoVesWEIoXw53eb07w+rIs7hxFcvdqm9EokKpV2SvOYztkZPJWF0ksdyPCXFqOdgWQ31XuTeOhaG6y8lgmq2zHUtICVHodBwb3E1eupkKeVUKeUU4GQ6efRUjD6ICJXSXjgROsGBybMZUfkHB1bnHrO9Ens0tr2+ikT7QWwX+Xa67yTFqKfBJskvV8ftBzDygbMYGFnAW+tGu9ReibQsUQte54rAzxg/Wd2d9S0REhhAUlSYKhkMmt6bLE0AACAASURBVGM+/XwG1Wxh37Kdx2ybq8Aeja0bGviwYTYlF9/qsQxvSDYasKuQmbk5I/51MjkBfdj1s2vJ6JRMIeQNrhiNIill8yIQWYD725k7EEKIxumpug9i34cdnfe+p46doTW3xEJceAjhXkxNP90wgEeDnmbEw77Z0HckakdQAYREBHP6RdEsXQo1x6j3Xm6pp7LW6tXsrWBTPiNLV5A9cbbquaZaIyVOXbcfwODHL8ZKAAeePfbMOLvEjD44gG4RnhvRBYuCeTrwUUbPvcBjGd7gDDhR814NCA5g3s17uGXPnVRUtN3WUm+lqKrOrzmnnLhiNNKFEN8JIa4WQlwFfAOsF0LMEkL43jHuI9QqUdqcXpNT2RZxAklpnyCP4bLJNpm9GrnV1MA73ySSffH9GLr758ZzPohqd3BXn5HPL9Xj2XTf5222UyLcNuPxhQRgp/cDsz2W4S1qBxgAdBsazxbjafRd9+kx71XnZklPw20bLA2Y3l3CrHPrMfpm7+lRpKgcdutk9pUB1NXBko/armundMVOb3DFaIQChcAUYCpQDMQC5wHnqqaZn0kxGlRLe9GcnOvnckXtu2w8RhCzsxypp/zxxK+cUbGI665WJ4zYFeIjQgkJ1KlujE+YlUBSYBFhn7zXZjtv3SjSLkn8aT7phuPpc9ZAj2QoQXKsnlJzPRU16oTdOqm76HJ6WXPY/Pa6NtvleDnA2fTkt3xYeQH/HuHeRk0liTUEExGirosaYOxYWBZ1BeMeOq3Ndu1ljwa4Fj11TRuva32hpD9INupVTXvhZMp/prJDP55332t9VGapt1JYWUeqF3s04t56kucDH2TKyT5IitQKOp0g2ahX3e2nC9SROeU6RpetYH9aVqvtnIvHvTxcCN/4u5lttQMonuXfCHTnYGK/ytd11JOzOMXwFy/9Ma7VNofDbb3o3D78gAJdIqPuP9NzGV4ihCA5Tv17VQiIOHksQ8zr2bN4W6vtsttBHQ0nbvUgQohNainS3nD6DtV2pURFwW3T9jL2g1uoPlTVYhvnIqenD2LOikxGlaeRffK1fvO7O1EzRUtzBj59NTZ0ZD30fqttckxmekSHeVwF8a0F4Vyt/4LRr/l37OQrt58+Ts+Qa8ax+EuBydRym0MVtTTYpMcDnMIthxhT+B0Z467yefjykSSrHHbrZNizl1NHMIeeav1ezTWZiQsPJiLU/8k43B12+rfH8SFqlig9ktlnlnCD9U0239+yD96Z+trTqWn2ox9gQ8fA/17lsY5KkWLUk2uyYFcx7BYg8fiebOx+FoP+/BBrrbXFNlklZlI9XFisOlTNxk93c8klEBnpjabe4wwZ9oUxnnNNA6/W3cDmG99q8fPDG9A8u667HviYQGykPOF/J0aqj1zUxgFGNvY6n2HbFlBX2fIWuBwv3dNK4q7R8F2GOz/TLTxE9bQXToZeP4G9IUOI/fKdFj9vWrD1YPRWX1XH4D8/YGP3s0gc28MrPZUg2WigzmqnsEr9PQ11t93DI/bH+f7box96KSXZxdUeG43N93zK5ppB3H6K5+VQlUIfHEh8ZIjqrhSAYaODmBCeTp9v5rW4IO7NAMdqBZn2G1uippByen+vdfUWX7moAUJuuZ5YWcqGB1rOfp1r8n6HvVK0ajSEEL2OfE9K+XDjZyepqVR7wBdpL5rOpRPknzOHIeb17Jy//qjPc01mjIZgIj2Ymn73Zi7l9kh0t9+mhKpe07RxUsXNaE4m3DeFH3rewLw3js5bVGp2hNt6YjSkXZK0+GUywkYy/NIhSqjqNb5y+wFUXDyHPvW72frKb0d9lmOyEBqkIz7S/XDbZcvglNpvOfjyYiXU9BpfuagBRv37VO5LmM99a8/nyNyTNfU2DlXUNiVS9DdtzTR+E0LcK4RociwKIeKFEAuAl9RXzf+onfaiOaNfvYYKoqh45IWjPssp8WyUISU8sXAAswbtYswDZyihptccdvupf12DguCumywM+fUVMr74+w5xZyhlajf3H8SN//2RfnU7KbnyLr+vETlJ8UGAgZMxz11MiYij/ukXj/rMmYnZk3DbN1+qISVFMO1y9Uu6usLhAY7696ouUEfKo1eyerP+qJowTWua7WCPBrRtNMYAfYHNQohThBB3AOuAtcB4Xyjnb3yR9sJJRFIEf514F2vzepKd9ffzeZoS/c8vD5Kx2cIdd+raTeeWFB1GUIDwWQd37eX1/JcHMd3//N/ez3IaDQ+uq/i/lzikS2LcC66l1fAFyUYDxVV1mOtaXr9RkrDYMHZMvY1xRcvZ93X63z7L9rCyXMZnm1m0OolnzvqNAP9VMv0bceHBGIIDfHavXnkl3KT/mNJr7v7b+0pWQVSCVo2GlLJMSnkj8B7wC3APcKKU8nUppborQ+0EX6S9aM6Qzx/lvsAXmffy4Q6+pt5GfkWtR4tgAbfexHbdSK6Y3X5+rgCdoFes+hsnnUQlR7NxxHWMy/qMQxsOp2nNLjETqBP0dDO77c6VhQwqXc3u028j2MsCQ0qS4oPd9s0Z/vatvBV4K29+cjgKwGaX5JXWeLRrufS+ZwjEyrT7RiipplccdlH75l41GODiERmcs28euSsOJ+HIUbAKohK0taYRLYR4G7gGmAYsBr4XQpziK+X8jS99mgA9esBll0Hm279QtK3gb+fu46YbZef89YwrWk7eqVcTZvDf3oyW8EWKlub0eeVfAOy55umm97KLzfQ26gl0M7vtI6/FMyR8P8Pf9k9OpNbwpdsPILa/kZ03v8arX/UitzF1Wn55DfU2u9uz4n1fpzMh7ws2nnA7Ucnu1e1QmxQf7NVozuC3bqeOEPLmPNH0XlaxY00zKsz/4bbQtntqE47a4GOllD9JKf8FXAE8JYRwPQl8B8a5SJpV7JsHEeCx6/L4qm4auy5/6m/ndtdomO/5DyZhZMxH7WMBvDnOFC1qVZs7kl6TU1l73HVM3PEOeX84erjsEjN93BwRb/6tgiVLJFffbSQ2OUINVT0muamuhu86uHvugfH8xdoLHWsbTpefu9e16LYnMWP4//bOPD6q6mz832eykp0kZIdAQiDsCAgC1qXWorVqW7VqW7W+trWt1b7aVm1/ffv6dtXaRbtbtXWttlotFhWtCKIICELYIRtkIXvINlknM+f3x9zBCCFMJnPvzGTO9/OZT2bu3Nz7nDvnnGc55zyH+Y/d4XcZx8pUI0Q9aPK0Ww8Z87N4b8ktLK98mspXDgJQ2WIfdfs3k5GUxjlKqV8opY4HSZVSJUqpFcCb5osWeDIS3QkCK5vtlt1z2jmTebf4Jpbv+TM1Gw8fv/doZvns+s1bnNn8KntWfYfEnODq3MDdEHsGnDTbrduWpejx77NFVvDQ/Z24lOJI6+jXaAxeeQ1roy7ljuDr20iMjSI9IYbDLdbV1cmT4Ufzn+Oq7XdS9uJeKprc9y6YlOD1NXa9WMkyw8tILQpQoqkRmJoWz6BLcdSiEDXA3CfupJcJNH7l+4DHwPH+mZrNSGMap9xJXin1sDniBBciQsGk+OMWlFUUPfkDBonk6DV3UNnSTU5yLHHR3q2OdThgzw9foDYin6VPBp+XAUOsYgum3XrIPjOPNd/ewE9emsfW3Qn0D7qYNoqGuPV7qzmzZS0xF3804Iv5TkXBpHhLvWKA+c98FzuJdHz5W1Q220mMjSTdy7Eelwtuvq+Ay1PfYfFzd5ssqW94LHwrn2v6rEm8dekv+fHRG3npVQct9oGQ8TQ0QOGkhOMWlFVkL8nlvYvv4az6f7F7Z82oLLff/x6ua32A3Q9vJW6M+4mbRaFRHis9OIAf/ADmZrfS8ojbUfbW07A32Mn7+a2Uxs5j5TPBqYjBqKsWP9PUojRKrvwxS1pfZ8c7FRRMSvB6uu3fftPC1q1w9QPLg9Ijhg+8Jquf6wX/uJlDhZfw7R+O3nszG600TkNBejx1HX30DJg/lXEoK5+/na3x51M1IGTFe9e5lb98iEfvLuMTnxAu/mKmyRL6Tm7KBGIibZY3xIQEePi6jSyM2gh4P0604/xvkeusof+BPxEVFxyDkcNROCmeth4Hx7oHLL3vR575OnsSltPS10WOl/NlqzdU8snbp/PzGY/w+cBllT8tqfHRTIyLosJiDy42Fn7/O8XSrqeA0Y9pmolWGqfBo+HNzqt/IlFxUXS/9gqD0cJb/45n0DHyoHFnbSfqiiv5t+MiHn1oEB+3MrAEm02Ylh5veUMEOOu+T7Mh7zziBnqp+d3JK5pP5B9/bifn4Do2LL2TeTevsEBC3wmUB2eLtCFPPUxTYirvrk9l4DQ6q7upG/sln0WAzz32cWxB3gsFwoMDWHWRMLHYRaRzkIZfvnHa8197Dd55x/yFkQH5uUQkVUT+IyJlxt9hN6wWkbUi0i4ia6yW0UOhsWFRIDq4+Gz3PadseYvNxTfiHBh+L4yelh7K53+aqf0HOfazh8jKC2x2UG8ozAhMQwRoXD6H7PY2iv/fFex9dOspz1v7quIL30jh1hU7WLn+xxZK6BuBiL8fZ0keAOUliXzz8iM4eobf26OvvY99865hZs9ODn3/SXKXT7FSSp8IxFiRh/5LlpPZ2c6i+65l1283nvK8kj+8y62XV/PUU1NwmrxlTqB0/N3AOqVUEbDO+Dwc9+Oe5hsw3CkRrLfe4IMpjNcUN/ORysfZkXcpzfs+vNNu1ZsVHMk/h4Vt69l6819YdOfIm7kEC4Xp7qmM/YPWbwrV5IBpi2fTGpmJ68s38+Tjrg/l+3EOOHnrit/Qccm1zJmteOblpKAOS3nImxhHdIT1YT/4IOZ/5xdc/GDtcvbkXfShxZQAVUcUu6ddxtKmNWy69vcs/VFgth0eLYWTEmix95u+ydVwVLX3UTBvOnXRUym+7ULe+PLfcQ2Z/atcirdveIRZt5zPXyNu4qc/3WP6ivpAKY3LgceN948DnxruJKXUOmD4TSYsIjYqgtyUCQGxNCqb7cRE2rhizbfZeM0fmN+8jti5hbwz/Ys8eGs5V10FL1/4AHk9pWz//mrO/lNA9euoKMxIwKWsW8HsYWDQRXOPYlZROrHbN/GzhX/n+i/auHTuYTYs/y4bzvwORxLncu4L32TaJDvrX+klJbjWm52SCJswNT0uIF5xZXM3InDX3VmU3/hTZrZtIfnMIt4puI4N593Df13RQfEs4Z/dF/Pubc9yzt8Cu3HVaCgIUNjP6VIcbu1mVkEqE/e+w76JZ/PwIzB/Pvzqi7vZsOwuyuIX8JEnvsz+lJXMKnmW1FTzFZtYtcDqQzcVaVdKpQz53KaUOlWI6jzg20qpU24tKyJfAb4CkJmZufjZZ5/1WTa73U5CwodnKvxyex+dA4r/WzG6lBNj5dfv99Ha6+LHZ7tnQbVvbiT+t8+zqGE9l6nVlKYv4aJlZVx9WQUJM5J9vs9wZTabIx1O7tncxy0LYzgzy7pwWp3dxffe6eXL86JZmRuF0ym8+moWsU+9zv823oFC2Bu3mPKPf4b0WxZjiwzygPsJ/HZnH0e7XNx7zodnzpn9G/+xpI/KDhf3n+u+b8f7LcQ+8AJnHF1Hhmri0sTXkfPm8LnPVZGVZc36HH+VuaHbxd1v9/KledGcnWudx9nc4+I7G3u5cU40506OwjXo4s31mbz4rzzOP/AUv1K3cyB2AYc++hnSb1+GLdI2pjKff/757yullpz2RKWUKS/c+ar2DvO6HGg/4dy2Ea5zHrDG2/suXrxYjYX169efdOyel/aq4u+/qpxO15iuPVrOu3+9+tpT24f9bnDQf/cZrsxmY+9zqPy71qjfvVlm6X3X7q1X+XetUTur2076zulwKqfDaak8/ubnaw+owu++rAYGP1wOs3/jTzy4UV3/6NZhvwvUM/VXmQcGnarwuy+r+1494Jfrecv6g40q/641aktFy0nfuZyuYZ/rWMoMbFde9LGmmXhKqVMG10WkUUSylVL1IpINNJ3q3GCgcFICvQ4nDZ195KRY420MDLqoPtbDJfOyh/0+WDKB+kp8TCTZybGWr4EpN+5XOMwUxlDzKoajID2BQZei+ljP8dlUZuNyKSqbu1k6LXXY70P9uUZF2MhPi7N8rOiDFEIn/45ik4Blrg7Ur/kS4Nl79AZgdYDk8IpAzEqpNlKy+7qzXCgQiKmM5U12UmMlKPZaNoPCDGMxmoXKuKGzj16HM6gWoPmbgkkJlo9rVraMboW9VQRKadwLXCgiZcCFxmdEZImIPOI5SUTeBp4DLhCRWhFZFQhhj89/tzCvj2fQLZgW9fgbz1RGZeG4WnmTnez4IF7EMkaOGzgWrivydKaF49zAOdLabVniQjByTo1ihb1VBGRCv1KqFbhgmOPbgS8N+RwU28p6Ehdaab0dzxg6jq23wkkJdPUP0tzVT0ZSrOn3c7kU5U12zs4J7XDJSCTFRjEpMcbiumqE/DLGb10tmBSPw6mobfNtvxBfqGzuZnlB8CVxHL+tx48EInFhaWMXmUkxQZND3ww8Hly5RSGquo5eeh1OcoJsfxF/U5BubV2taLKTEBNJRuLo9wUPFQotzkHV2eegvqOP6ZnBp4jHd+vxIwXp1q4KLWu0MyMzOJO4+QurV9t7BsFzEsZ3tS/MSKC8yW5Z2K+00c70jOALo/iTQovHNcsa3XV1Rkbw9QHju/X4kcJJCRxt77VkD2ZPGGX6OHb3AbKSYomLjrAslBIuSqMoI4GOXodl+5WUNnYxIwgtYn+SEhdNWnw0ZU3WrDUua3TfJxgNx/HdevxIkfHjlVvQwR1td4dRgrHC+BMRsXQGVXmTnbT4aBKjx69FDB90NB5r1Uxa7f20dg+M+7oK7udaasEzBbf3NiEqYtR72FuBVhpeMjPL3SgONZpvaZQetzLGt/UG7ud6qMEa6628yT6uB2s9eDpwK56rpxMtCgOlMTMrkbLGLlwu88N+ZU1dTM9IwBagtRgjoZWGl0xJjSMm0kaphQ1xehDGM/3NzMxEmrr6aTN5DwilFGVhEPIDSE9w7wFhRSjFc49wMHBmZCbSPeC0ZOvX0sYuioK0rmql4SURNqEoM8EST6MsDGZOeZhhkQfXYh+go9cRtA3Rn4gIMzKt8eBKG7tIjIkky4Ip04FmZpa77pSaXFc7eh00dvYHrfemlcYocMc0LWiITV1hESMGt6cB5jdEz1hUOHga4AmlmD+DqrTRTlHm+J455aHoeF01d1yjLMjD01ppjIKZmYk0dvbT3mNeKMUzc6ooDEJTAJlJMSTFRppuFXuUUrg816LMRLr6B6nv6DPtHkopyhq7jo/3jXeSYqPITo413cApMwycYDUctdIYBZ5QipmWRm1bL30OV9BaGf5GRJiZZb4Hd7Chi4lxUWQmjd8FaEPxeHBmhv1a7AO09TjCRhEDloT9Shu7mGDs4xOMaKUxCqxoiMct4jBRGvDBDCozQykHGzopzkoKizAKfBDaKDOxrgbzWgKzmJmVSHmz3dQcVGVGyC8YZ06BVhqjIjs5lsSYSFNnUJUas1HCYeaUh5mZiXT2DdLYac5iNJdLcaihi+Ls8HmmKXHRZCTGcKjBPK84nKaGe5iRmcjAoIuqY+btOOmeORW8dVUrjVEgIszISjTV0zjU0EV2cmxYzJzyMMNkD66mrYeeASfFYRJ792B22K+0yU7yBHeCxHDh+MQNkwzHtu4Bmrr6g1oRa6UxSjwzqMwKpeyv62R2dpIp1w5WZpjcEA/Uu69bnBVez7UoI5GyJvMWox2o72RmVmLYhPwAI8eWeQbOgfpOAGbnBG9d1UpjlMzMTKC9x0FTl/9DKX0OJxXN9qCuMGYwMd4dSjloktI42NCJSHjF3gGKsxLpc5gTSnG6FAfru5gTZnV1QnQE+alxpg2G7zeUxqwgNhy10hglxcaP6bEI/Mmhhi5cirBriOAOpRxs8P8zBThY38W0tHgmRIf4HrmjxGN87Kvr8Pu1D7d00+twhp1XDO4Ofb8J7R9gX10nmUkxpCcEb8hPK41R8kFD9H+l8VTE2dnJfr92sDMnJ5nSxi4GBv0/K+VgQ2dYDYJ7mJGZSFSEsPeoeXV1Tk741dW5uclUtfbQ2efw+7X313UG/TPVSmOUJMVGMTUtjr1H/W+97a/rJDEmMigzW5rN3NwkHE7l94Hb7v5Bqo71hN14BkB0pI2ijERTPI19dR1ERUjYrLAfisdw3O9nw7HP4aS82R703ptWGj4wJyeZvSY0xP31nczKTgra+dlmMtewrvzdwR1s6EIpwm7mlIe5uUnsr+v0+8SN/XWdzMhMJDoy/LoQT131t+FY1mjH6VJBP6YZfr+4H5iTm0TNsV46evznnrpcigP1nUFfYcxiSmociTGRfg+leBr2vLzgdvnNYk5OMq3dA7T1+09pKKXCcpafh0mJMWQmxfjd09hf766rwf5ctdLwgeNWcb3/LI2qY+61BMFeYczCZhPm5Cb53YPbXdtBekJMWGRhHQ7PpIqqTv+NFTV1uTdeCscJGx7MiDbsr+skPjqCKalxfr2uv9FKwwc8jWWfH61iT1gmXD0NcCvjA/Wdfk3RsOdoO/PzksNqLcFQZmUnIQLVflQaHgt7Tm54em8Ac3OSKG+y0zvg9Ns199WFRnhaKw0fSEuIISc51q+Wxu7aDqIjbGGVc+pE5uYm0+dwUdnS7Zfr9QwMUt5kZ24Yd27xMZFMS4/niB+VhifkF67jROBWmC6F36aJDzpd7KvrDAnvTSsNH5mdk+zXgbCSmnZm5SQRExleawmGMjfX3WD21Prnue6v68SlYH4YKw1wh1L8GZ4qqWlnekYCibHhk+rmRDyd+14/jWuUNtrpdThZOCXFL9czE600fGRebjKVLd10+WGu9qDTxZ7aDs6YHPwVxkympScwISqCPX5Sxrtrw3sQ3MOCvGSO9SmaOse+t4ZSipKadhaGeV3NTZlAanw0u2va/XK9XbXu6yycPNEv1zMTrTR85IwpKSgFu2rG3sEdtzLCvCFG2IR5ecmU+Kkh7jnaQUZiDJlhOgjuYVG+uyPaUT3251pzrJfW7oGwr6siwhmTU9hR3eaX65VUt5MS514DFuxopeEjC6ekIIJfKo3HylgQ5g0RYHH+RPbVddDnGPsA447qtrDv3MAdSokU2OmHurqzxn2NM0IgjGI2i/InUtHc7ZedPEtq2lmQlxISEza00vCRpNgoijIS/KI0QsnKMJvFUybicKoxh6iau/qpau1hydTgd/fNJiYygvwkm3/qak07E6IijqcID2c8inPnGD1je/8gpU1dIWPgaKUxBhZNmciOqrYxp54OJSvDbDwN8f2qsXVw71cdA2BxfuqYZRoPFKbY2F3bgWOM05l3VrczLzeZyAjddSzIS8EmsHOMdXV3bTtKERKD4KCVxphYlD+Rzr5BKlt83x2to8dBaVMXi6Zoixjc05mnpcePWWlsP9JGdKTt+IyscGd6SgT9g64xZWfuczjZX9cZMp2b2cTHRFKclTTmsaL3j7jr+sK80HiuWmmMAU9Hv6PK90rz3pFjKAXLCrRF7GHRlInsrG4bU76k7VVtLMhLDuspzEMpTHE39bEo453V7Qw4XSybpuuqh0X5KZTUtOMcQ7Rh6+FjFGclMjE+2o+SmYdWGmOgID2eiXFRbD18zOdrbKlsJTrSFjLxTCtYnD+RFvsAh31c5NfncLKvrkOHpoaQNsFG3sQJbK5o9fkaWypbsQksmaqfq4czp6Zi7x/0OQ/VwKCL7VXHOKsgzc+SmYdWGmPAZhNWFKbzbkWLz1bx1sOtLJqSQmyUtog9LC90N6B3fezgdla343AqFufrkN9QVhams6Wy1WereHNlK3NyksNq//rT4amrmypafPr/3bXt9DlcnBVCkQatNMbIiulp1Hf0+WQVd/Q62FfXybJpoWNlWMHUtDhyUyawqdy3hripvIUIm+iQ3wmsmJ5GZ59vVnGfw0lJdfvxTlLjJiMxlqKMBJ/rqidKsTSE+gCtNMbIysJ0AJ8qzXZjPCOUXFMrEBFWTk/j3QrfrOK3y1tYkJdMUhinuRiO5QUeD270dXVHdRsDztCyiK1i5fR0th05Rv/g6NcWbalspTgrkdQQGc+AACkNEUkVkf+ISJnx96Q4gogsFJHNIrJPRHaLyNWBkPV05B+3ikcfSnm7rIWYSJteKDUMK6enG57Y6NZrdPQ42FPbztlFk0ySLHTJSIplekYCm3wI+22ucI9nnKnHM05iRWEafQ4XO0c5i6rP4WTbkdAaz4DAeRp3A+uUUkXAOuPzifQA1yul5gAXAQ+ISND1riLCisI0No8yVqyU4s2DTaycnq7HM4ZhheHBvTNKD25zZQsuBWdPTzdDrJBnZWEa2w4fG/WK+3UHmliSnxrWSQpPxbKCNGwy+mjD5opW+hwuzi/OMEkycwiU0rgceNx4/zjwqRNPUEqVKqXKjPd1QBMQlObjuTMn0dHrGNV0xormbqqP9YRchbGKSYkxzMpOYv3BplH93xsHmkiMjdTe2yk4vziDXodzVCGquvZe9td38tFZuq4OR/KEKBbnT+SNA6Orq28ebGJCVETITWGODNB9M5VS9QBKqXoRGbE2ishSIBqoOMX3XwG+ApCZmcmGDRt8Fsxut4/6/yMGFZE2eHTtNnpmxXj1P68edmfHjWurYMOGw6MV06/4UmYrmBk/wOpyB6tfW09yzOlXyztdirW7e5iXHsGmtzee8rxgLa+ZeMrscCliI+CxN0qwNXhXV9+sdtfVpK4qNmyoMVNMv2Ll71wY4+DZIwM898qbTIo7vS2ulOKVkl6KJ9rYsultv8lhSZmVUqa8gDeAvcO8LgfaTzi3bYTrZAOHgLO8ue/ixYvVWFi/fr1P/3fjX99TK+9dp1wul1fnX/GHTWrVr9/y6V7+xtcym82B+g6Vf9ca9fSWKq/O31TerPLvWqNe2V034nnBWl4zcI7ERAAADQZJREFUGVrmrz/9vlr8o/8op9O7unr9o1vVR+570+u6HSxY+TsfabGr/LvWqIc3Vnh1/p7adpV/1xr1zFbv6ra3jKXMwHblRR9rWnhKKfUxpdTcYV6rgUYRyQYw/g7r14lIEvAy8H2l1BazZPUHq+ZkUtvWyz4vpjPWHOthe1Ubly7IsUCy0GVmZiJT0+JYu6/Bq/Nf29tATKSNc2cGZRQzaPj47Exa7P28d+T0i1Kbu/p5p7yFT8zL1rnRRiA/LZ7irETW7vWurq4uOUpUhLBqTpbJkvmfQI1pvATcYLy/AVh94gkiEg28CDyhlHrOQtl8YtWcLKIjbDz/fu1pz31pVx0Al2mlMSIiwifmZbOpvIWGjpE3EOofdLJ6Vx0fm5VJXHSgoq6hwYWzM0mIifSqrq7ZXYfTpfjMolwLJAttLluYw/aqNiqbR85F53QpXtpVx7kzMkImdchQAqU07gUuFJEy4ELjMyKyREQeMc75LHAO8EURKTFeCwMj7ulJiYtm1dwsXtx5dMSZKUopXtx5lCX5E5mcqlOhn46rz5yM06V4bvvIsfTX9zXS3uPgmqWTLZIsdImLjuSSedm8sqcee//giOe+uPMoc3KSmKFToZ+WKxflEWETnjuNMn63ooXGzn4+fUZoKuKAKA2lVKtS6gKlVJHx95hxfLtS6kvG+6eUUlFKqYVDXiWBkNdbrl4ymY5ex4gu6qbyVsqb7Fx9pu7cvCE/LZ6V09N4dlvNiFOan91WTW7KhOOLLTUjc9WSPHoGnKwuOXrKc3ZWt7G7toOrFudZKFnokpEUy3kzJvH8+7UMDJ46Bf3j7x4hNT6aC0J0NppeEe5HVhSmUTApnj+9VXHKPTYefruS9IRoPZ4xCr6wLJ+j7b28vKd+2O931bSzqbyVzy2bgs2m4+7esDh/IvNyk/nzxkoGT7HHxsNvV5IUG8mVS7SB4y3Xr5hKc1c//9wxvLdR0Wxn3cEmvrBsSsiuz9JKw4/YbMI3zp/OwYYuXtl7cge3uaKVt0qbuensgpCtMIFg1ZwsZmQm8OAbpSdtIqSU4lf/KSV5QhTXL88PkIShh4hwy/mFVLX28MKOk72NXTXtvLKngeuW55MQo8eIvOWconTm5yXzuzfL6R04OUx9/9pDTIiK4LrlU60Xzk9opeFnLluQw+zsJH747/109DiOH+8dcPKD1XvJTo7lxpVTAydgCGKzCXeuKqaiuZs/bvjwUp2XdtXxVmkzt350ul6tPEo+PjuLxfkT+dmrB2ju6j9+vH/Qyf+s3kt6QgxfPbcwgBKGHiLC3RcXc7S9lwfWlX7ouzf2N7J2XwNfO7eQSYnerZEJRrTS8DORETbuu2I+bT0DfOmJbXT0OOgZGOT2v5dQ3mzn51fO116GD3xsdiaXLsjhgTdKecFw/TdXtHL3P/ewYHIKX1wxNbAChiA2m/DTT8+jz+HiS49vo9XeT5/DyV3P72Z3bQc//tRcrYh9YEVhOtecOZmH3qrkb1urUUqxs7qNO/5RwuzsJL58TkGgRRwT2u80gXl5yTxw9Rnc9uxOlt+7jggRuvoH+Z9PzuYjOpGez9x3xTyaOvu44x+7uPfVgzR19VOQHs8j1y/Re1b7yMysRH577Rl8/ekdrLzvTaIjbHT2DfKdVTO5aG7orSEIFv7v8jnUdfTxvRf38Os3Smmx95OTPIGHrlsc8kajVhomccn8bKamx/H3bTUMuhSfOSNX73g2RuKiI3nypmX8fXsNO6ramJ6RwA0rpuqY+xj52OxM1tx2Nn/bWk2fw8llC3JYoRM+jomYyAj+csMSnnu/lm2HjzE1PZ4blk8lOS70PTfd2kxkTk4yP7w8OdBijCuiI21cd1Y+152lB739yYzMRO65bE6gxRhXREbYuHbpFK5dOiXQovgV7dNrNBqNxmu00tBoNBqN12ilodFoNBqv0UpDo9FoNF6jlYZGo9FovEYrDY1Go9F4jVYaGo1Go/EarTQ0Go1G4zXi3hp2/CAizUDVGC6RDrT4SZxQIdzKHG7lBV3mcGEsZc5XSp02z9G4UxpjRUS2K6WWBFoOKwm3ModbeUGXOVywosw6PKXRaDQar9FKQ6PRaDReo5XGyfw50AIEgHArc7iVF3SZwwXTy6zHNDQajUbjNdrT0Gg0Go3XaKVhICIXicghESkXkbsDLY/ZiMhkEVkvIgdEZJ+IfDPQMlmFiESIyE4RWRNoWaxARFJE5HkROWj83ssDLZPZiMjtRr3eKyLPiEhsoGXyNyLyFxFpEpG9Q46lish/RKTM+DvR3/fVSgN3JwL8HrgYmA1cKyKzAyuV6QwC31JKzQLOAm4JgzJ7+CZwINBCWMiDwFqlVDGwgHFedhHJBW4Dliil5gIRwDWBlcoUHgMuOuHY3cA6pVQRsM747Fe00nCzFChXSlUqpQaAZ4HLAyyTqSil6pVSO4z3Xbg7ktzASmU+IpIHXAI8EmhZrEBEkoBzgEcBlFIDSqn2wEplCZHABBGJBOKAugDL43eUUhuBYyccvhx43Hj/OPApf99XKw03uUDNkM+1hEEH6kFEpgJnAFsDK4klPADcCbgCLYhFFADNwF+NkNwjIhIfaKHMRCl1FPgFUA3UAx1KqdcDK5VlZCql6sFtGAIZ/r6BVhpuZJhjYTGtTEQSgH8C/62U6gy0PGYiIp8EmpRS7wdaFguJBBYBf1RKnQF0Y0LIIpgw4viXA9OAHCBeRL4QWKnGD1ppuKkFJg/5nMc4dGdPRESicCuMp5VSLwRaHgtYCVwmIkdwhyA/KiJPBVYk06kFapVSHi/yedxKZDzzMeCwUqpZKeUAXgBWBFgmq2gUkWwA42+Tv2+glYabbUCRiEwTkWjcg2YvBVgmUxERwR3nPqCU+lWg5bECpdR3lVJ5SqmpuH/jN5VS49oCVUo1ADUiMtM4dAGwP4AiWUE1cJaIxBn1/ALG+eD/EF4CbjDe3wCs9vcNIv19wVBEKTUoIt8AXsM90+IvSql9ARbLbFYC1wF7RKTEOPY9pdQrAZRJYw63Ak8bBlElcGOA5TEVpdRWEXke2IF7luBOxuHqcBF5BjgPSBeRWuB/gXuBf4jITbiV51V+v69eEa7RaDQab9HhKY1Go9F4jVYaGo1Go/EarTQ0Go1G4zVaaWg0Go3Ga7TS0Gg0Go3XaKWhGReISJqIlBivBhE5OuTzuxbKMXVo1lErERH7KY5PEJG3jMScp/rfN8zIiKoZf+h1GppxgVKqFVgIICL3AHal1C8CKlTw8F/AC0op5wjnPAl8HfiJNSJpQhXtaWjGPR4LXETOMyzuf4hIqYjcKyKfF5H3RGSPiBQa500SkX+KyDbjtXKYa0aIyP3G97tF5OZhzpkqIm+LyA7jtWKIHBtF5EUR2S8ifxIRm3HNx4w9IPaIyO3G+YUislZE3jeuV2wcnyYimw0ZfjTCI/g8xspgEck27l1i3OcjxjkvAdf6/pQ14YL2NDThxgJgFu6U0pXAI0qppcYmVLcC/417/4lfK6XeEZEpuDMFzDrhOjfhzp56pojEAJtE5HU+nOiyCbhQKdUnIkXAM8AS47uluPduqQLWAp8BDgO5xh4QiEiKce6fga8qpcpEZBnwB+Cjhpx/VEo9ISK3DFdYYxV4gVLqiHHoc8BrSqmfGOGqOAClVJuIxIhImuG1aTTDopWGJtzY5kkdLSIVgCdl9h7gfOP9x4DZ7rRFACSJSKKx74iHjwPzReRK43MyUASUDjknCvidiCwEnMCMId+9p5SqNOR4Bjgb96Y5BSLyW+Bl4HUjC/EK4Lkh8sQYf1cCVxjvnwTuG6a86cDQ/TO2AX8xklX+SylVMuS7JtxZYbXS0JwSrTQ04Ub/kPeuIZ9dfNAebMBypVTvCNcR4Fal1GsfOujem8TD7UAjbu/GBvQN+e7E/D3KsPYXAKuAW4DP4vZ82pVSC08hx+nyAPUCx7c6VUptFJFzcG9E9aSI3K+UesL4OtY4X6M5JXpMQ6M5mdeBb3g+GJ7CibwGfM2w2BGRGcNsbpQM1CulXLiTQw6dvbTUGJOwAVcD74hIOmBTSv0T+B9gkbHHyWERucq4jxiKBWATH2xj+vnhCqKUagMixNgjW0Tyce8p8jDuLMeLPNcFsoAjIz8aTbijlYZGczK3AUuMAe79wFeHOecR3CnGdxhTbB/iZM/9D8ANIrIFd2iqe8h3m3FnJN2LeyzjRdy7RW4wsg4/BnzXOPfzwE0isgvYxwdbEX8T997u23ArqFPxOu7wF7izopaIyE7coa0HjeOLgS1KqcERrqPR6Cy3Go3ViMh5wLeVUp+06H5nAHcopa4b4ZwHgZeUUuuskEkTumhPQ6MZ5yildgLrR1rcB+zVCkPjDdrT0Gg0Go3XaE9Do9FoNF6jlYZGo9FovEYrDY1Go9F4jVYaGo1Go/EarTQ0Go1G4zVaaWg0Go3Ga/4/Fd3Bd4WlK+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x205a58df4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write Task 3 code here...\n",
    "\n",
    "### BEGIN SOLUTION\n",
    "\n",
    "def ODE(boundary_conds,t):\n",
    "    x = boundary_conds[0]\n",
    "    v = boundary_conds[1]\n",
    "    \n",
    "    # just needs to be a list with 2 items the empty lists are just used for None\n",
    "    x_dot = [[],[]]\n",
    "    \n",
    "    x_dot[0] = v\n",
    "    x_dot[1] = -5*x\n",
    "    return x_dot\n",
    "\n",
    "task3=scipy.integrate.odeint(ODE, [0.1,0], np.linspace(0,10,1000))\n",
    "true_values = [0.1*np.cos(np.sqrt(5)*t) for t in np.linspace(0,10,1000)]\n",
    "\n",
    "# plotting the x values wrt T \n",
    "plt.plot(np.linspace(0,10,1000),task3[:,0], c='b')\n",
    "plt.plot(np.linspace(0,10,1000),true_values, linestyle='--', c='r')\n",
    "plt.ylabel('X-position (m)')\n",
    "\n",
    "# velocities wrt T\n",
    "plt.plot(np.linspace(0,10,1000),task3[:,1])\n",
    "plt.xlabel('Time elapsed (s)')\n",
    "plt.grid()\n",
    "\n",
    "### END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}