{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cruise Control Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import necessary libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from collections import deque"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function allows us to ensure the correctness of symbolic manipulations. This is helpful to obtain the analytical solution of the cruise control example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def equals(a,b):\n",
    "    assert np.isclose(float(simplify(a - b)), 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following defines the symbolic equations for the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKkAAAArCAYAAAADraj8AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGfUlEQVR4Ae2djXXUOBSFMzkpgIUOQgewdAAdhN0KFjogJXCyHcBWsJvtIKGCQDpIqICfDsL9PJIiazyOZcsTO346x7EtPT09XV09SR7LWd3c3OxZKIfAarU6lLaXwvVjOa0PX5NwO1ItL4XbdVrb/TSi5L0KfqbjTMcPHR9K6i6tS/Y9imy96qMfHcr3wQjajp5w2uCFMPtfuU6URievhZUSaxFj3KhgCnntDBmjiGI6ZesXKfssW9/mKnV5qeeGN8jVtQT5lBe6p5N/En7P4/qP6kkpSAW/dAWexwWXvlY5RzreFdD7TDrOcvWo7DfKA7knT9CCWOXCFORlwwYvhN1PCZyn7XgQco138Uqqr50B45VSQHMTcBlqjyVLXS10Q2AbL94r+1cdf3s1o3tSFUSPGdWL+soUOG8DrlW1yF3No+bgRVsrstvERl44Z3YdOYy9op7UNRYeBbf9TQfkZPikd1RBMsw7GBqf6riSUaHHKI0Fy4ni7mtlXANO9rDifKyDORJ2bRvKkRulI04Nrz72KA+duJUXSo8DWL7WscZUwO+VOKSQBv6h49Dr0zVzOxZNj6I4GltRa3kf7+Ig6Wkc1/Va+SDKu67yqZzy0nkw7Ig0Bea3gEu9iHiT5vH3SqOeW9O9XJ+z9BbHSzp7Y5Vrj+Q78SLGRnlwYl98XJHh3vWuUyk+luLY23DNsy88K4soGv2CawV6yvfq6vYPvS3Of5sy/hVgEvzE/aOry2fFvdV1m3enXmld0DUoTA2vXHu68qIBJLDEaVThgL9O2ad1VOe/f6nhLp30ic4o/S/JTcPz/MuH78rj7/9QZJgGOAHI7Ens89TOspWyPKHiNMp/rPQ/40h3zcKNTtEWmI9SPp4rPH7SNXFtBEUnUwLkNsJAbAfhNQJWufZ05UWKG46Kjr8OaoTBw700Mcyfxbp0D2kYJvn1pVYGcS4tTAOQURgyXPcewlzZ1VQDG3QwfHe2RbLUf6Oeab373qNbBwAVwUt6hmLVyR6Vk8ULj4/ysY7R7Zo3B7oZFJyngJDeq3p9VISSmhYU21bRT3zmXZ5dHei5DOsM93j7K52ZqjTZn5pXG57SxAL3k8JL9bnTnp688FAxMoVp376PLXBOh2mG3Yq4MpgH7bfue+3Ka6RGRvL/FrCjjwqmHnGHgnQEOh/TIYattgCgADtWALsp4ZVjTw4vPH7oD9OnfR/b9yxPgzK8DYqroEZldYbL9r3hheT8NTJcx/KQAZlaQyC4o4BnCB7T1YmiY5vbTMHu2k95bcI90qaG15329OSFhwYsQ3sU+e1epIRk/+jwvYYCIC8eiLiNYVN5ThVPZZkL4sXuWpwgtjVIH56Yx1/huetW4SRBeZmD8mKIX9ThPZmbVs9ydWalT30ag2SZ2pAf+VGCyiiGl3T1xspXros9ksnmBfqVj/cnwsK8CEm94fd5LgH8EPtVPp3tlYja1fsOKW5Q3vvGqs142QaxeUYaOvxDIimVwyNv9Xht4AxNE7hMcZ6r/PD4aqjOsfI7ItwbVm31km2Mvt+EYxgRHwxJ2yq+qzQBzLSBJwST96a7wiSnHNd5dv+qXo6RD0CWHwwm/XL3xDFm3r3xo4t50sKtJm/AUwvbPpKJq3BjunTeNAoZSTPBNPHdIzD4OenuTbYSl4aAkXRpLT7D+q5kc/VmxwxtN5MXgoDNSRfS0HOupg33c269hdhuJF1IQ8+5mkbSObfeQmw3ki6koedcTSPpnFtvIbYbSRfS0HOuppF0zq23ENuNpAtp6DlX00g659ZbiO2DtzQvBKdQTb1Sxn4mNu4d6uCLKwT2DBFe6LjgrXLJsRGRHbN8E4t49kCFzWW6t9ARASNpR6AiMfYxHYuE1eY9xfNRjGqrg+LYwsJXrRFn82FFYt1DYl7o/Y0EC3kIGEkz8BLZ8I4XLgue9KcnqIvze++fKD72msRXe7CcnJ0yEDCSZoAlUb4pdek8JiRl2I8DJCak37iq7etfi9jfrgjYwqkrUpITQf1OVOal3MfekijI2PRVa+SZHljogYCRtAdoygIZm762AhnDByZQLa9LHEN9Fe+8MEkWOiJgJO0IVCLGt6NqXjSaAqQek92PLKL4xDZThMoLJ/rstgUBI2kLOE1Jjmh4xpSMFfkapgC/S9YTmj35NU/bVIbF1REwktbx6HKHN2Te6Ynn8xDf9D0rFlH8IzO27NqefI9WxvkX/4Nyhllp2tkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{- d v + k \\left(- v + vd\\right)}{m}$"
      ],
      "text/plain": [
       "-d⋅v + k⋅(-v + vd)\n",
       "──────────────────\n",
       "        m         "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t, v, m, k, vd, d, v0 = symbols('t, v m k vd d v0')\n",
    "dv_dt = (1/m)*(k*(vd-v) - d*v)\n",
    "dv_dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rewrite the previous equation into the form $av + b$, so that we can define the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAArCAYAAABB7ttrAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGOUlEQVR4Ae2cj3XcKBDGvX4uwC/pwOkgfzpIOnAuFSTpICnhntNBchXcSzqwr4Ik7sC+Cnx2B77vh4EFCa2EvCuxu/AeFhoGmG9mGJAsdnF/f39Q03o0sFgsTtTTa+n023p6nLcX4TmVBJfCc52S5NARxXisfK58q3zl6Ju6aoznwXhfNzXOVP0Ky7HG+jrEcebEjpyB3lfaWVh+CNOZ+JkU7SSGgzCL47cySojom7rXWIS+0031P1W/Vm8nOePNid3K22tn8TEpfqdw+cgjBpeeq3DubjZ5lUe/tv1fbHKcTfctHB80xi8pOBneU+M/Brvanip/SvWbQRtkZ2G6U58XqfGOwsEeAyjsJ6P8RrzXVsCMZsWxfpZEYMlJs2EfYec/Bexf5S8hwGbkmRoQkWfbo47ZD+REHWuAObFn2dlO7uvA6QyEKPKIEgESM7vtJ8ovlM+cgkRnHSRUP1O+Et17pOrYhMEbPXGIjpKZoYTBG2WchtCJV8+axuAJBEZHKydAgdgH2TnASBGMb+2Ve+1W7cZYRRzCb15VZk3F4AwE/UPAi3OIZOpuHd3ScJ7vDRp93Cr7DaXK7Kvo5DjknaMsGbLwhDJaHF43YR1lpbVjV5847KfmWEPu1W6wncP+1I5gEW2cj0RwCZAktzn6psZ3mjX/ifZRZRNJ7Cz6aTgfPJH6MBFdXjmCndXfdf9ZfYQbSsq8QyASzZZy8SQEZYI1dWDYCsU+yM4JnGDE8ZZJxjPRRxTetRAdVj6+qd5HCssfzQDREM4/eqtMv60II5pb3nz0c7JMeZUcWXiasqk9OuPFYAuH6I/CrvZnyrw6aWZ0x7hNOvdR1G/KpfpBdk60Y4sh8hKnB6wKs9zoynLFkhI5RdiIshJOknKKqJ14AHkette9C51JpYe8U5Ul0yA8TXksviQOW7d27Or3MctWlp0dXo3Zcp4jEQ9seCX8sjyxbPFm8UpXlpWuzWDXjv0pfZJsvzjKpSEs/2AoHLer7yXndKVePB2itMO5GEvEbmXKtbODzYNTuO04OLQ1f3ANjOnWcLPGaVDCZzMhROQU4mNG/N1k1L3bI7mqdyqYtrRRpq+5Uw6eUFYUimK7UknYx9jZ4UI/d+6G66G9Ydb5KCAnckyRp1led6HOG10OgKO9UlvvULYf+g352LUTAl3ftHFlkWdLvXg6JAMvrzKiVCj2MXZ2uMDofQTiQiAJsexx2CizXJkkGnsf8x5HV/Pk9VCz/CsenqJQOusokSt6twPNOtVfKroZiAA4J9EM2qqlUdXTpSF4mtKoDUswukNXUdoUdvVLhOe1h3+/Fg3ccaN2o+xMd2rLZvy9xvTBwThPx1iVPFADUiyT540UO0kEHes8A+G02DQeqwrveKIJUp2npap8gpTLUvxCyv2Y3zq/hTUmkd5tL/I7yWih8VglbjReFOmq82QocRWrFMySwNPqJNFnlSzrrLOO+o9wtfZ1bsO8zvH2tS/+78MLuF1L7GvB1ko18rRUMp6gWcpT5S59hspyfNEVTavzjPeVvW9Zl629d4HxCqjOM153e99yIQ2Y/3LOrQmtq8jSmdyuv5MhXRG91EqxqN8i8KdkK51W9zylW6hg+eqyVbBxShetOk/pFipYvuo8BRundNGq85RuoYLlq85TsHFKF606T+kWKli+6jwFG6d00arzlG6hguWrzlOwcUoXzRy96RNSr/D5TpePp0+UORFK4jtaEqdDf+rfC1/Ex4ftnIy4sXS+7Y0+mha9JquBrderjMvnjCuzsLqz3HwtR/YH+1Tm+1b+P8QH8/7wm8o4V3SOvW+cfauXfrZar72Rx0YTd/KByHMnI4ffsrozS09FD6MMdByrpoQGdkGvvc4j3Pz40qXA4gg4D8tXmFiqSM2fSonOCD2w1L+BBrZer70bZjmO+0KffQ/LWxhdIOEkqV/3gp8lrqaEBnZBr73OE+DGSfyBr4COk/jDgtDtRpBIZeg2alFVU1sDW6vXHOfhnHMUdYKlrBlh+Nre/H6veFjqTNRq661SpIGt1esg57EOQCRpOknXUvZSvM7ROMsURabqMg8a2Ha9DnIeQSV6sK9xDuHsD711Pl00Ns/8WDRHN3bxLJNgrSVttV7/B/9NG/o62/GKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{k vd}{m} - \\frac{v \\left(d + k\\right)}{m}$"
      ],
      "text/plain": [
       "k⋅vd   v⋅(d + k)\n",
       "──── - ─────────\n",
       " m         m    "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = -(1/m)*(k+d)\n",
    "b = (1/m)*(k*vd)\n",
    "dv_dt_norm = (a*v + b)\n",
    "equals(dv_dt, dv_dt_norm) # Proves correct rewrite.\n",
    "dv_dt_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define analytical solution, and prove that it is indeed a solution by differentiating it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAAyCAYAAAAaw83gAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAN6UlEQVR4Ae2djbXVNhLH8zgpIJtUsKQDAhUEOiCkAqADOFSwh3RAqCCBDshWAKEDSAWB1wH7/ykar6xr+/pDsmVfzTl+19bHzGhmNBp92O/qy5cvX1XIK4Grq6vnovBGsv4jL6WKfQ0JSJ/fiM5L6fOnNehdGo0bl9bgtdvrHdJX1SGtLfl89KTLa2F/Id3+mY/K5WIe7ZQYHXS90fVZ14ecIhP+WwGtFzlp5cStNjwR/psy4qc56fThFv2qsz7hLEz3g0x1TAvl2FX9SsLtSu9Nk6EzOrxTvce9hRJliBbM/SRarxOhXA2NeL8vYoT4/1qNaA+hqrMewSRIlmzfCM176XmTgSdBE4pDMTpSCji/pXsUkRWk7LuewO7WYcS7W3MQ/6WsOVSd5bNWdPwosNd8lC4E8ySnFAh+DUdxTzr4qBGI+fve4JUY/l28ryGnQdlUnQ2KZ3Gmt8+HQoTOKySQwCSnJHprOgoipc079VQZywkwbYP3UsL5qrOpSpxYXo6J5YVP0j27rLsEBi9dRNQnoHTWznrXkZU32G5f/2RTQOlPdDGraMFUp9RyFEJ4XxehK0zfBLN+WVyFGGks9Dag5w+6HjUJ/kZpN3VR/rku6iKcVaYcMS8JnlHQrwVFeGd1Rpsl80l6U/kj6SyB2t0ghO26fpACYWocXmcnEZ3n+ZZs9n0PTeoMBQgnjiXCw3LPuyiNHelflHbq0JRB5tlLlSHMwvN9ygvczpJ+MXoSHvn05/6X9M8hbj3jbV9Faa6c0tmlcnzonkaA8xtL28Ov+MXhwnfTli35Rn6en0GdwaNgtN5U9jA6S6kfyQX7fpESZ0pc4g37fBLjVFqrT3bk4zhc/47zeBYMtpl8Xc4G4/qkx3lnnZEh8ZWdo9A9Dsk5DH51mUNilLAOACMfrL5nHgac8ftn6n7W1Wqwnqn7Z1h3D/fiGaPkkORoueYsK16Q96DOoC8YrTeVPZTOUspfsrFBqbjBVLzhWJzThE9rt+7R54lToYwu+jkX9XoHWuWd1Df83r5cfZVjMANvq7yeW05xdOcBkS4cSAthRLxRhi/b8spKg6nGY4JLl+s0ER4a0TivMK/Ue/HLdJO23C2FR/FyVmfwKhitN4/zEDpLrSfkqAvZtOw+NZ25+MRXK0jwuqdPtvjVcxw8xDMeHBy2ZZc5O3tu+rin4QIPlWfwQ0aNU/T5rYF8ilOCMHNLPCfTq1ZDQkEpj4Z2GW6rjsrAbJuh/yu2mM4dtq3vXu1wDqAvf4t08TRaZ/AnOKs3lTmMznLoRPJhQbe4KF884QxO+FIakUvsRGiDi4z0y2DbimRiuSl/KFDBpsCH3+gsR36I84YSzoIWwmgQXg6kLE5xcJJFaQh2Qd+Oz3dW2OMEb7y45nCKztDCmqEp6feBmCmG5xk6Q5aDejugznLYz29CyhsJ9JeSgH7lFpvFG47I4KNuYl7pl598Aeq9VZ2+vu6L9f7g1H7zfuM2NqSLSCwEo+XSRjkllaTDEVVZpzMkMM/ODeFcCDSy5Ww8IygshrdRws96dnWpoysWWFR8+0fxiOCRRfZDpRNaO1VnoB6rt93rbIIcpxa1tw/mduKp9MaWtz5FxGL9mLo4qjvcBEDQ8Ux2bc7rez3P7YfUNZkwm3im61pXCK3nr8OcgXtG0KYhck7XYpjieNkuIL1RisrSYe+onjGHgwMHOJvGeiHQwc2hteoovVSwtjYyKoDRqTqD5UG9HUxnWVQkGX2UHdPJkP+vWYjMQApfqnbyapjXaQuj0rDjJLYsXA1N3Z/IQ7IiamoFK6PefVNFIgCmbo1TURoeFy/IusXJuRzls/6EIMjHCXUxhLN6qctGXgSBQom8SOOdoiTCEa5s4OVzW7xu/p6bNXKOzqh7Tm/KP4TOTE45fr3sWZOhfxQP4pfAgHUllmYmg+oPnXEaxKe6+JXGcVF4lFMaxFoz6cgs/vKSMqNjhQuXgOyBQZVNHTed2IM4xLOL9sXzakGAaBLYnAQ0Y6dve5DrJjxKsEQOXESFFaoEkICbHcg2ZkcQKcUoPtzW6hicKjumWMoybJi18F2sU5IgmMsSYs8KWQMp2pqYM8Qgvd5mkEBCvWXgrkHJ4jFA9GHroy5h7h+1G1xj3qfkUz8sgTSwp4gNpi/WKTUaW35jTqlGSstleRQMZgvJ1pT8tGq1qdWWirixJfGD0P7Wt8MM8SDNqs2YK4EgUjHbmIvqIuvVSGm52llPAloh8z9J9W8pEtD0Bz1x7sbtGMtxNNN25TH15rWmkx3ihfybbZxF4/njDM/fvjB8wtPFDXbVKZ01l7MF7JT6p7Mla4EtJfBMHfypX5vhuErjlDxTzbkilWFXCL3iIL6jni8z5YdBalSk5HniYOE9c0JKY1mArfI5tFUtHYgXnPkPujgaZG3iGZ7hk8OXnNpOsn5WnZKkuRBGj4YL6dTqMyXgO7idhePztfEAQsd3p5p9B2wckZ55q4B3v6g3Fc7ahnBTBif5MHBIpMETnb4EQF5Ekzh2nJEdg+G91ddqA0k/60rjlIRw9HYhlEsDCaW9nxgxqPZxZuRulMwjiv9W+QgzBj7DO9YI3cih8tcxkq5n0YPuf7vyBtIw2JbCq94m6e2T5GcHfx9Izv+JZI3uzGnhDJrDfNSTrF+ht7E69rjpyBZV+KTOHw4PA9giERoAP08n0nMVc/wxGQj3Q/AjC/2EMiXKxLGmAREc/aWAI5WV9DgS0PpqwZz2CQfKUNXLlOPa7V6iN9VlcOr9eoXy6Gzkc76o6Rd65nBs6036ML/rXuWJLFqf/OgpB73Bt/C76q2dRnuMpu7pO5zEdjKyPNItbclv3X2TJBcCI6KNHgtR1eqZJcCI3vXPKFg/AlgfAYhUQhgb9YR1uHe2ESd2PBe9mO0jo/cB38gpjIze+ygvLBMUn3Zb15SmyaurtBkwobrdd5WradtLgM7U6jjqTCcvhKpMlzMhipoC2MMYZ9Nbho6uiCNekJ/CQ5Ky4gG7bpYzYp703OSlIFgjpeVStC3cqUa7nHLFMFUCOACLhiy65UsULUelMvFaUPw8lm6vwwkQsJjdWvMkMtFFuq2DBcWPf1sjpeU6tuhoruEu56BiGCUBOR+OBLBozeYHaz6sibCwbWBOJB5geI4dl9Xp+6VOV8TVKi/6v4gfHCROiC80OlB6s9huaZfyO+iUJKhbEgQKvK2Lf664qqBEH8Uyd4U+q/3Jju0LXyo4a3ipCI3Fs6XeSteZbKh3qqE8vvGFY+oaYOx9trNq8DKg3PXZwioguptP0cbwuVaZQackYTE63JOQ2SHg4NSqgJGIIPQZQUYbxUgmU71HZKNrMy0YST9bsS31lllnyCyV3vrkT8SC43J0ZHusOb327eqrE6ebLdjUPs6vzwMSOLumJKXYfHeyMaBQXXb2YoCNs1lEbEmdIkY20dA6mRQOC+vdobLOQhskztVbyTpDjKn01qcS4Sdq4Z+m8kkNZgmcsu6NrnrwmFOyAaunWE3uksBgpOQr9G2jduFLnja3cyVnZBghxscUsyTYTG870Vmvrrxj6s0fkXHHl7EBa0SVWsQkcDZSUkEipclRkhFI8LtZ55rAO/IhmisJttTbHnSWU1cuUpJzq5HSDCm3IiWNcAiT3QjWcpgPW2eLj+UrazVodS7xyByfhUimS6W8Re12TcTb1FcRkgixQL3tQWdJZN+DpNX+njI1uUcCjVOSYSNIdrp+MA+vNFvH2SRSopOLHyIQ5xT1zPoUZzdwnraFmvpzE0I9GWwRHhnC32pQmt52pLMsOlL7sU3s1vpOFjpHRuqckjckHBIvAYYhJ/f8RxEipy2ATg78IR5xSO4j47pnG/6x+CrBIbH4yjF7ZMTLvas5pUL1tgudSU+5gEgeWM0O/iF3nL/OKak57DLg3X+PmoaBnRWuOgf1zRhDFOBc8iY+axN0dqZpzRkp7ySLcEhBY5Gd/fO+IDnr7Wy9VZ1l04utp4WDezZiR0Ts/sWSDPTkXwT5UZh0tkRnTd+Eg1Fj9sf5VZ9Tt+xgvIUPXXy/pciDZuIVp0zIPlteqjsJRDO53i5JZ5OEPbKw5MeZPmYcRdrpyGZsWuxrCZFohivevnSRz1yHtLRVni/m50zTmL4RsXF+hOnkLCe5lKeh+p5HRsfm4N1Q+aV5JeptbzpbqoO4vtpvZ/JKi+JjVot+vhFwRzQSAusjzlFJ2ByCxEGsCQ8gFjgg1pEAHCjvCjF1KQ3g6ZF4czyuxFxJetujzlKq6ZmQTT39nZL+IXDdUKe/VkuIPBqno07F2gi7XjYv5k1qu1fyKsB0rYmIPJ8QXpuP0Y0Vj4yQyBPjzAqF6m13OkulJPUZlioYjMIXfFOhvyg8Fikx5bgjwT7h0v07XRgYEQnPW2xvomC2/UNA4Y89T1uenQp5iu8fKgE5wn9uKE1ve9VZCj0xELE7XOygmaKRa+BwC925CKljLlrozsVXbrxqt3uBWAba7BjmppkK/6XqbIn8vMxeCse/fQS7BN3F183tlFy0cGmKkpEyFcYx/ai2xxsIRRudeL9Inc1VipfXX6rPv8tulhvm4qv1NDuTIKscMkjAj56cryrxG1AZWnyZKKVnDh3zra/dRcWlaszWlErlb7d8yUg5wsBRBoy2wgElIN2yIcT7jtUhJdRvdUoJhRmjMmOV8bJZUOFAEpBOOcfHGTq3IXSgpm3elOqUMqtARssO2ffeiDNTq+jXkIB0ybob0dGPa9C7NBr/A6KTYmS/JtiKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{k vd}{d + k} + \\left(- \\frac{k vd}{d + k} + v_{0}\\right) e^{- \\frac{t \\left(d + k\\right)}{m}}$"
      ],
      "text/plain": [
       "                        -t⋅(d + k) \n",
       "                        ───────────\n",
       " k⋅vd   ⎛   k⋅vd     ⎞       m     \n",
       "───── + ⎜- ───── + v₀⎟⋅ℯ           \n",
       "d + k   ⎝  d + k     ⎠             "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_t = exp(a*t)*(v0 + b/a) - b/a\n",
    "equals(dv_dt.subs(v, v_t), diff(v_t, t))\n",
    "equals(v_t.subs(t, 0), v0)\n",
    "v_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following represents the explicit euler method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Explicit_Euler:\n",
    "    def simulate(self, f, x0, t0, h, max_steps):\n",
    "        x = x0\n",
    "        xs = [x]\n",
    "        for i in range(1, max_steps):\n",
    "            x = x + f(x)*h\n",
    "            xs.append(x)\n",
    "        return xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following appplies the previous numerical solver to the cruise control example. \n",
    "The results are plotting against the analytical solution `v_t` obtained earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VPXZ//H3nclKIIR9S0ICouybIYqIdUGq1orWpbhbrWi1j9rWVm37a221fWytba1Pqw8ulbohohbriru4sYOgURBBEraw79nv3x8ZfCJkGUImJ8l8XteVK5mZc+Z8EuXcc77nu5i7IyIisSsu6AAiIhIsFQIRkRinQiAiEuNUCEREYpwKgYhIjFMhEBGJcSoEIiIxToVARCTGqRCIiMS4+KADRKJz586enZ0ddAwRkRZl/vz5m9y9S33btYhCkJ2dzbx584KOISLSopjZl5Fsp6YhEZEYp0IgIhLjVAhERGJc1O8RmFkImAescffTzSwHmAp0BBYAF7t76cG+b1lZGYWFhRQXFzdu4BiVnJxMRkYGCQkJQUcRkSbWFDeLrwfygbTw4z8Af3H3qWZ2H3AFcO/BvmlhYSHt2rUjOzsbM2u8tDHI3dm8eTOFhYXk5OQEHUdEmlhUm4bMLAP4FvBA+LEBJwLTw5tMAc5syHsXFxfTqVMnFYFGYGZ06tRJV1ciMSra9wj+CvwMqAw/7gRsc/fy8ONCoFdNO5rZJDObZ2bzNm7cWOObqwg0Hv0tRZqhgjkw666q71EUtUJgZqcDRe4+v/rTNWxa41qZ7j7Z3XPdPbdLl3rHQ4iINI1ITs6NsU3BHJhyBrzxu6rvUSwG0bxHMAY4w8xOA5KpukfwVyDdzOLDVwUZwNooZhARqVIwB1bNguyxkJnXsG32nZwrSiGUCJc+d8B2vno2/GsCVJTioQS+/NZUtnQcxu6SCvaUlrO7pII2RQsYN/f7hCrLqIhL4OHD/sayhAHsKauguLSCvWUVnLrtcS4oLyFEZdXxVs2qPfchilohcPdbgFsAzOx44EZ3v9DMngLOoarn0KXAjGhlEJFWIKAT+LJvPsb6tKHsLC5nZ3E5O4rLGLjiGcaET84V5SVMf+pxHk8qY1dxGbtLKthdUs4lFU/zo1AJ8VZJRZnz1PQn+EfFnq8d65rQ85wcX0acVVJZUcbe5W/zbmIXUhJDpCRUfX2ROoKKPVOJoxwLJVb9blESxBQTNwFTzex2YCHwYAAZDtlNN91E7969ueaaawC49dZbadeuHT/5yU8CTibSgjTCCXz/bfySGezsOpKtu0vZuqeMrbtL6bL4BQaWlxAXPoG/9Nw0preB7XvL2L63jB17y5hY8hQ3xP3fCfy5f0/jHxUlXzvUkXHdGJUYTwLllBPPwtBg0pLj6ZWeTGpiPKlJ8fQqOxn/+N9UVpZDKIHjTjyTo3qMIjUxRJvEeFKTQrTblEboqefwilLi4xO57tLvcd0Bv/8xUDCo/iLYCJqkELj7W8Bb4Z+/ABr1N/rNfz7mk7U7GvMtGdgzjV9/e1Ctr0+cOJEbbrjhq0Iwbdo0Xn755UbNINKiNcZJftUsvKIU8wq8opQ1C2eydHsmm3aVsnlXKZt2lZBb8CSnhz+ll5eX8NfJD/I/5Zu+9jYjLZ3HwifwMovnhR192VxZSvuUBHq2TyEtJYGu5Sfjn1adwC0+gZPGn81JGaNol5xAWnIC7ZLjaZN4GlY4BlbNIj57LHfUeHIeBHm9YdUs4rLHcnRN23QaW/X71neSz8yLagHYp0VMOtccjRgxgqKiItauXcvGjRvp0KEDWVlZQccSaRqN0RSzchaET/KVFaUsffd53u7agaKdJRTtLKZoZwldtqVwd2WIBJwyD3HdB6kseH/BV++R3iaBbcn9OIUEoIxKS6D7sJP5ZfcBpLdJpGNqAultEunQ5njKt+aRtO5DUnLGcm+NJ9chUND7q9/ryEM5OTfWNk2kVRSCuj65R9M555zD9OnTWb9+PRMnTgwkg0ija+yTfHkpi975DzM7prF++17Wbi9m/fZiuu+IZ0ro/07yt37UgQW+jPYpCXRtl0TXtCRS+x7DU/ydAaWL2dNzNDdkHkWntol0bptEx9REEkLhjo8FI2HVLELZY7motpNr5zHQb0zdv3szOjk3pVZRCIIyceJErrzySjZt2sTbb78ddByR+kWhuWb57Bd5e2VnCrfuoWDrXgq37qHT1gQesvBJnhC3L+3Ikrgv6JaWTM/2KQzPTKdH+njeqOxNv72LqOx9LHf3PYYu7ZJITgjtF2o4cF7dv1eMnsAbiwrBIRg0aBA7d+6kV69e9OjRI+g4InWL8CRPRSmET/LFy94i3/uxatPuqq/Ne0han8ZvPUSCV32Sv3l+exZ4Pu2S4sno2IbenVLJOOwkXrMsjihehGWP5b4jxtA5NYm4uP2HEg2ggZMLSCNSIThES5YsCTqCSJX6Pu1XO8lX75deXlHJ6i17WF60i52bcziDeOLCzTUXvpbAglffByDOoGd6Cjmdh/JI13sYXrGEyt5j+O1hY8js0Ia0lPj9RqgPQif5lkGFQKQlaIQmnfKsMcSFEqAcKiyeu5d15dV577By025KK/bNApPEK+1uZVzKMrZ1O4rTso7imk6pZHdOJbNjCknx+5ptjgJ0X6y1UCEQae4OskmHilKKP3+bj8oPI3/dDvLX7eCTdTv4bP1OBlXczNFx+XzoA9i8JYPDuqRwfP8uHNalLf26taNvl1TaJWsq8lijQiAStAY26eyzbU8pK+MGM9jiiQu3218wM54F/gEAHdokMLBnGhcf3ZsBPYZwRPd2XNe1bQ03ZSVWqRCIBCmST/vZY6teqyjFQ4l8ljyM995dyUeF21hcsI1Vm6umLxhpt3BquxVs73YUJ/UZzX/1SGNAjzS6pSVpdlmpkwqBSJDq+bS/fnsxc7ZksD7nryQWvscLO/oy9+kS4BO6pyUzLLM9543KZHhGOoMzxpOmZh1pABUCkWiqr9lnv0/7a9OP5N25q5m9cgtzV22hYMteANokpjMi60LyhqdzZUY6wzLT6ZaW3MS/jLRWKgQi0VJPs4+7syplEPlH/i+ly99hxvYc3nxsN7CEjqmJjMruwKWjs8nL6cjAHmnEh6K9jpTEKhWCZmbVqlUMGDCAI4444qvnfvzjH3PJJZfUus+tt95K27ZtufHGG5siouzTgJu8WzsO570Vm3h3+SZmLd/Emm17gXh6pX+Low7vyH/ndGRUdkf6dklVu740mdgqBJHMa94M9O3bl0WLFkXt/cvLy4mPj63/9I0uwpu8HkrEK0qpsHh+Pi+N6S++iju0S47nmL6duPr4vow9rDO9O7XRiV8CEztng0j+4R6Eg1mP4M4772TatGmUlJRw1lln8Zvf/Ia5c+dyxRVXMGfOHCoqKsjLy+PJJ5+kbdu2tR6zbdu27Nq1C4Dp06fz/PPP8/DDD39tmxUrVnDttdeyceNG2rRpw/3330///v257LLL6NixIwsXLmTkyJHcddddDf7dhTpv8m7fW8ZbnxXx6ifxbCv7BUPLlzCXgdB1EDec1IVj+3VmWEZ7NfVIsxG1QmBmycA7QFL4ONPd/ddm9jDwDWB7eNPL3D16H3/3qad3xsGKdD2CmTNnsnz5cubMmYO7c8YZZ/DOO+9w3HHHccYZZ/DLX/6SvXv3ctFFFzF48GBWrVrFihUrGD58+Ffvcc899zB2bGSrE02aNIn77ruPfv36MXv2bK655hreeOMNAJYtW8Zrr71GKKT+4/U6iJu8hBIp6jSKF99byav5G5j9xRbKK53ObRM5acg3GD7gXH7Qt5MGakmzFc0rghLgRHffZWYJwLtm9lL4tZ+6+/QoHvtA+/3DPdRl3yJdj2DmzJnMnDmTESNGALBr1y6WL1/Occcdx69+9StGjRpFcnIyf/vb377ap6FNQ7t27eL999/n3HPP/eq5kpL/W2Hp3HPPVRGIRCRXj5l5fHn6VFYvfIVntuTw7L92Ap/Qt0sq3x/bh5MHdmN4ZjqhAyZZE2l+orlmsQO7wg8Twl8erePVKzMvshWBDkIk6xG4O7fccgtXXXXVAa9t2bKFXbt2UVZWRnFxMampqXUer3obcnFx8QGvV1ZWkp6eXmsRqe/9JayOq8eCLXt4bvFa/rN4LZ+u302cHUtu7478fFRXxg3oRp8utTftiTRXUW2kNLOQmS0CioBX3X12+KXfmdlHZvYXM0uKZoavycyDsT9ptBvFEydOZOrUqUyfPp1zzjmnxm2++c1v8tBDD33Vtr9mzRqKioqAqmac2267jQsvvJCbbrqp3uN169aN/Px8KisrefbZZw94PS0tjZycHJ566imgqggtXry4ob9e61UwB2bdVfW9JvuuHi0EoUS2dj2Kh99byVn/eI+xf3yTO1/5jNSkeH47YRBzfjGOaVePZtJxfVUEpMWK6s1id68AhptZOvCsmQ0GbgHWA4nAZKoWs//t/vua2SRgEtBsl4CMZD2C8ePHk5+fz+jRo4GqG76PPvooL7/8MvHx8VxwwQVUVFRwzDHH8MYbb9CnT58D7hFcfvnlXHfdddxxxx2cfvrpZGZmMnjw4K+KS3WPPfYYP/jBD7j99tspKytj4sSJDBs2LDp/gJYowmafPec/y+dzX+KZzTn86+HtVPp2+ndvx02n9Ofbw3qQ0aFNMPlFosCqWnCa4EBmvwZ2u/ufqj13PHCju59e1765ubk+b968rz2Xn5/PgAEDohE1ZsXE33TWXfDG76qafSwEJ/6i6iqR8BVU4XamzlnNc4vXsqe0gqyObZgwvCdnDOtJv27tAg4vcnDMbL6759a3XTR7DXUBytx9m5mlAOOAP5hZD3dfZ1UN3mcCS6OVQeQANXQa2L63jBmL1vD47NV8un4nKQkhvj2sB98dlcXIrHT175dWL5pNQz2AKWYWoupexDR3f97M3ggXCQMWAVdHMUOTWbJkCRdffPHXnktKSmL27Nm17CFRUV+3z3CnAV85i0+Th3H/B4m8uOQ1issqGdwrjdvPHMyE4T3V1VNiSjR7DX0EjKjh+RMb8RjN5tPakCFDojoaONqaqokwqiJo/99bWsHTa7oxZV4uy4t20TZpA98ZmcH5o7IYktE+oOAiwWqxI4uTk5PZvHkznTp1ajbFoKVydzZv3kxycgufzbKObp8bd5bwyAereOTDL9m6p4whvdrzx7OH8q2hPUhNarH/DEQaRYv9F5CRkUFhYSEbN24MOkqrkJycTEZGRtAxDk0N7f/LN+zkwXdX8szCNZSWVzJuQFeuHNuHvJyO+gAhEtZiC0FCQgI5OTlBx5CmVtc9gGrt/0sShvCX1+DNz94hKT6Oc47M4Ipjc+irvv4iB2ixhUBiUATz+7+yPYt7FuXx8doddErdzo/GHc5FR2fRqW3TjVsUaWlUCKTlqOUegLvzxqdF/PnVZXy8dgc5nVO54ztDOHNELy3QLhIBFQJpOfa7B+C9j+WdZRv586vLWFywjayObfjTucM4c3hPTfEschBUCKTlqDZx4EfxQ/jtCxXM+3IOvdJTuOM7Qzj7yAwSVABEDpoKgTQfEawgN7fiMO7Kr+DDL7bQLW0Pt00YxHmjMkmKVxOQSEOpEEjzUM+N4IIte7j9hU945eMNdG6bxK9OH8gFR2XpHoBII1AhkOahlhvBu0vK+cdbn3P/rJWEzPjJyYfz/bF9SElUARBpLCoE0jzUcCP43wsLueOlT9mwo4Qzh/fk5lMH0L19Cx/9LNIMqRBI81DtRvDylOHc9Hw5C1YvZmhGe/5x4UiO7N0x6IQirZYKgTQbRelD+eP6JKbPL6RzW+eP5wzlnJEZxGndX5GoUiGQplFHj6DKSufxOau546VPKSmv4Kpv9OGHJxymqaBFmogKgURfHT2CVm/ew01Pf8QHX2xmzGGduP3MIeR0Tg04sEhsUSGQ6KuhR1Blr1E8/P4q7nzlM+LjjDu+M4TvjsrUjKAiAYjmUpXJwDtAUvg4093912aWA0wFOgILgIvdvTRaOaQZ2K9HUGH6kdzwvx8w78utHH9EF35/1hB6pqcEnVIkZkXziqAEONHdd5lZAvCumb0E/Bj4i7tPNbP7gCuAe6OYQ4IW7hFU8cU7/Gd7H256ci/JCaX8+bxhnDWil64CRAIWzaUqHdgVfpgQ/nLgROCC8PNTgFtRIWj1liUO4KdLSllcuJ3xA7tw+5mD6ZqmMQEizUFU7xGEF66fDxwG/B1YAWxz9/LwJoVAr2hmkGC5O//64Et+90I+bZPjuef8EZw+tIeuAkSakagWAnevAIabWTrwLDCgps1q2tfMJgGTALKysqKWURpBLV1Dt+8t4+anP+Klpes5sX9X7jxnqBaIEWmGmqTXkLtvM7O3gKOBdDOLD18VZABra9lnMjAZIDc3t8ZiIc1ALV1DFxds44dPLGDdtmJ+flp/vn9sHw0ME2mmojZ5u5l1CV8JYGYpwDggH3gTOCe82aXAjGhlkCawX9dQXzmLh95dyTn3vU9lJTx51WgmHddXRUCkGYvmFUEPYEr4PkEcMM3dnzezT4CpZnY7sBB4MIoZJNqqdQ31UAJ/+LQz933xCeMGdONP5w4lvU1i0AlFpB7R7DX0ETCihue/AGpedURannDX0LWLZnLrko68saoLv/xWf644Nkc3hEVaCI0slkPi7jy4qjN3fDCSbmnJPHX1CEZkdQg6logcBBUCabDS8kpueWYJTy8oZPzAbtx5zjDat9FEcSItjQqB1K2WrqFbd5dy1aPzmbNyCzeM68f1J/VTU5BIC6VCILWrpWvoio27uOLhuazdXszdE4czYbjGBIq0ZCoEUrsaZg19v7QPVz8yn4RQHE9ceZRWDhNpBVQIpHb7zRr62t5+XP3gHHI6p/LQZaPI7Ngm6IQi0ghUCKR24a6hlStn8ei6TH71Roix/Trx9wtHkqbVw0RaDRUCqdOebiO54c04Zn6ygYuOzuLWbw8iPhS1AekiEgAVAqnV5l0lXPbPuXy8dju/On0g3xuTrZ5BIq2QCoHUqGhHMRc+MJvVW/Zw/yW5nDSgW9CRRCRKVAjkAGu37eXCB2azYUcx//zeKI7p2znoSCISRSoEsayGwWKrN+/hggc+ZPueMh65Ik/dQ0VigApBrKphsNiK5IFceP9sissrePzKoxmS0T7olCLSBFQIYtV+g8U2Lnmd7y7YgbvzxJVHM6BHWtAJRaSJqBDEqmqDxSrjErhxbltCCfDY90dzWNe2QacTkSakQhCrwoPF1iyayU3z27MyeRDTrjyK3p1Sg04mIk0smktVZprZm2aWb2Yfm9n14edvNbM1ZrYo/HVatDJI3WaX9WX83FwKUgcz7erRKgIiMSqaVwTlwE/cfYGZtQPmm9mr4df+4u5/iuKxpR4LVm/lsn/OpVeHFB77/lF0S0sOOpKIBKTeQmBmXYExQE9gL7AUmOfulXXt5+7rgHXhn3eaWT6g+YqbgeUbdnL5w3PpmpbE41ceRdd2KgIisazWpiEzO8HMXgFeAE6lajH6gcAvgSVm9hszi6hriZllU7V+8ezwUz80s4/M7CEz07qGTWjNtr1c8tAcEkJxPHK5ioCI1H1FcBpwpbuv3v8FM4sHTgdOBp6u6wBm1ja8zQ3uvsPM7gVuAzz8/S7g8hr2mwRMAsjKyorol5FqahgstmV3KZc8OJtdJeVMu2o0WZ00jbSIgLl79N7cLAF4HnjF3f9cw+vZwPPuPriu98nNzfV58+ZFJWOrVMNgsd1dR3LBA7P5dN0OHrniKPJyNGJYpLUzs/nunlvfdvX2GjKz680szao8aGYLzGx8BPsZ8CCQX70ImFmPapudRdU9B2lM+w0WK//iHa5+dD5L12znfy4YqSIgIl8TSffRy919BzAe6AJ8D7gjgv3GABcDJ+7XVfSPZrbEzD4CTgB+1MDsUpt9g8UshIcSufvzbsxavon//s4QTh6oWURF5Osi6T66bwL604B/uvtii2BSend/t9q+1b14EPmkIcKDxXzlLB4o7MU9H7Xj5lP7c15uZtDJRKQZiqQQzDezmUAOcEt4TECdXUelGcjM455lHfjzR8uYdFwfrv5G36ATiUgzVWshMLN4dy8HrgCGA1+4+x4z60RV85A0Y49++CV/fnUZZ4/M4OZT+gcdR0SasbquCD40s0LgZeBld98G4O6bgc1NEU4a5t3lm/jVjKWc2L8rd5w9hLg4LS8pIrWrtRC4e66Z9aZqMNlfzawX8C7wEvC2u5c0UUY5CKs37+GHTyygX9d23HP+CBK00LyI1KPOs4S7f+nu97n7mcAxwH+AccAsM3uhKQJK5HaXlDPpkXlUVjqTLzmS1CRNLisi9Yv4TOHuZcAb4S/CVwjSTLg7P52+mGUbdvLw9/I0k6iIRCySAWWnm9lCM9tqZjvMbKeZ7XD3NU0RUCLzj7dW8OKS9dx8an+OO7xL0HFEpAWJ5Irgr8B3gCUezfkoJHL7zSP05qdF/GnmZ5wxrCdXju0TdDoRaWEiKQQFwFIVgWZiv3mE1kyYynVPlzCgexp/OHsoEYz1ExH5mkgKwc+AF83sbeCrnkI1TSInTaDaPEJeUcrLz08nITSByZccSUpiKOh0ItICRVIIfgfsApKBxOjGkXqF5xHyilLKiOelXYfx98tHktFBU0qLSMNEUgg6unu9s41KEwnPI/Tea8/y52VdOeP0CYzu2ynoVCLSgkUy2ui1SKadlqbzyo4sLvrsWPqMPJFLj8kOOo6ItHCRFIJrgZfNbG/17qPRDiY1K9iyhxunLWZYRntuP3Owbg6LyCGrt2nI3ds1RRCpX3lFJT96chEA/3PBSJITdHNYRA5dXYvXZ9e1Y3jFsozGDiS1u/etFcz7ciu3nTmYzI66OSwijaOupqE7zexpM7vEzAaZWVczyzKzE83sNuA9YEBtO5tZppm9aWb5ZvaxmV0ffr6jmb1qZsvD3zs08u/UKi0q2MZfX1/OhOE9OXOEZvcQkcZTayFw93OB/wccAfwdmAXMAL4PfAac6O6v1vHe5cBP3H0AcDRwrZkNBG4GXnf3fsDr4cdSh90l5Vw/dSHd05L57YTBQccRkVamznsE7v4J8IuGvLG7rwPWhX/eaWb5QC9gAnB8eLMpwFvATQ05Rqz4zX8+pmDLHqZOGk37lISg44hIK9Mkk9WH7zeMAGYD3cJFYl+x6NoUGVqql5asY9q8Qn5wfF/ycjoGHUdEWqGoT1hvZm2Bp4Eb3H1HpN0dzWwSMAkgKysregGbm2oTyq1PG8rNzyxhaEZ7bhh3eNDJRKSVimohMLMEqorAY+7+TPjpDWbWw93XmVkPoKimfd19MjAZIDc3NzYmvKs2oZyHErm3w+8pLe/J3RO10piIRE9di9ePrGtHd19Q1+tW9dH/QSB/vwnqngMuBe4If58RcdrWrvqEcuWlpK79kF9P+CU5nbXIjIhET11XBHeFvycDucBiwIChVLX1H1vPe48BLgaWmNmi8HM/p6oATDOzK4DVwLkNi94KVZtQrqQyhGeP4bujMoNOJSKtXF2L158AYGZTgUnuviT8eDBwY31v7O7vUlU4anLSwUeNAZl5lFzwLI8++Rjvlw/gTxecrykkRCTqIrlH0H9fEQBw96VmNjyKmWLa75e0Y8r2U3nkijw6pGrWbxGJvkgKQb6ZPQA8CjhwEZAf1VQxav6XW5nywZdcdkw2Y/tp3WERaRqRFILvAT8Arg8/fge4N2qJYlRZRSU/f2YJPdsn89NvHhF0HBGJIZHMPlpsZvcBL7r7Z02QKSY9MGsln23YyQOX5JKaFPXhHSIiX6m3c7qZnQEsAl4OPx5uZs9FO1gsWb15D3e/voxTBnVn3MBuQccRkRgTySilXwN5wDYAd18EZEcxU0xxd345YynxcXHcesagoOOISAyKpBCUu/v2qCeJUc9/tI53lm3kxvGH0719ctBxRCQGRdIYvdTMLgBCZtYPuA54P7qxYsP2vWX85j+fMDSjPRePzg46jojEqEiuCP4LGASUAI8D24EbohkqVvzx5U/ZsruE3581hFCcBo6JSDAi6TW0B/iFmf3e3Xc3QabWq9rMovMr+/HY7NV8/9gcBvdqH3QyEYlh9RYCMzsGeABoC2SZ2TDgKne/JtrhWpX9ZhZ9JOm39Gzflx+drOmlRSRYkTQN/QX4JrAZwN0XA8dFM1SrtN/Moj23zeO3EwZrzICIBC6iSe7dvWC/pyqikKV12zezqIUo8RCefazGDIhIsxDJx9GCcPOQm1kiVb2GNNfQwcrMwy+ZwVNPT2XGtj7cdd55QScSEQEiuyK4GriWqoXn1wDDw4/lID2/NZOfbRjHyeNP15gBEWk2Iuk1tAm4sAmytGp7Syu4/YVPGNJLYwZEpHmJZK6hPmb2HzPbaGZFZjbDzPpEsN9D4e2XVnvuVjNbY2aLwl+nHeov0FI8MOsLNuwo4dffHqgxAyLSrETSNPQ4MA3oAfQEngKeiGC/h4FTanj+L+4+PPz1YqRBW7KNO0u47+0VfHNQN3KzOwYdR0TkayIpBObuj7h7efhr3wI1dXL3d4Ath5ywFbj79WUUl1fys1P6Bx1FROQAkRSCN83sZjPLNrPeZvYz4AUz62hmDfl4+0Mz+yjcdNShAfu3KJ8X7eKJOQVckJdF3y5tg44jInKASLqPfjf8/ar9nr+cqiuDeu8XVHMvcFt4v9uAu8LvcwAzmwRMAsjKyjqIQzQvf3j5U1ISQlw/rl/QUUREahRJr6GcxjqYu2/Y97OZ3Q88X8e2k4HJALm5ufU2RTVHc1Zu4dVPNnDj+MPp3DYp6DgiIjWKpNfQuWbWLvzzL83sGTMb0ZCDmVmPag/PApbWtm1L5+78/sV8uqclc8WxB3PRJCLStCK5R/D/3H2nmR1L1ZxDU4D76tvJzJ4APgCOMLNCM7sC+KOZLTGzj4ATgB8dQvZm7YUl61hUsI0fjz+clMRQ0HFERGoVyT2CffMKfQu4191nmNmt9e3k7ufX8PSDB5GtxSopr+CPL39G/+7tOHtkRtBxRETqFMkVwRoz+1/gPOBFM0uKcL+Y9eiHq1m9ZQ83n9pfg8em9m6nAAANAUlEQVREpNmL5IrgPKoGhv3J3beF2/l/Gt1YLVB40ZldPUZzzxu7OPawznzj8C5BpxIRqVekK5Q9U+3xOmBdNEO1ONUWnUmyePoU38Itp12Oma4GRKT506oojaHaojNW6VyRsYZBPbX8pIi0DGrrbwzhRWcqiKOMePKOnxB0IhGRiKkQNIbMPFac+jh3lZ3L9MH/oMvAsUEnEhGJmJqGGsmti1JZmnQOb3/7hKCjiIgcFF0RNIKFq7cya/kmrvpGX9KSE4KOIyJyUFQIGsHf3/yc9DYJXHR076CjiIgcNBWCQ/Tx2u28ll/E5WNyaJukljYRaXlUCA7R39/8nHZJ8Vx6THbQUUREGkSF4BAs37CTl5au55JjetM+RfcGRKRlUiE4BP94awXJ8SEuH9NoSzaIiDQ5FYIG+nLzbmYsWsOFR2XRSYvOiEgLpkLQQPe+tYL4UByTjtOiMyLSsqkQNMCabXt5ekEhE0dl0jUtOeg4IiKHJGqFwMweMrMiM1ta7bmOZvaqmS0Pf+8QreNH0+S3V+AOV32jb9BRREQOWTSvCB6mah2D6m4GXnf3fsDr4cctStHOYp6YW8DZIzPolZ4SdBwRkUMWtULg7u8AW/Z7egJVax4T/n5mtI4fLQ/MWkl5RSU/OF5XAyLSOjT1PYJu4YVt9i1w07WJj39Ituwu5dEPv+SMYT3J7pwadBwRkUbRbOdEMLNJwCSArKysYMOEl6F8uSiLPaVtuPaEw4LNIyLSiJq6EGwwsx7uvi689nFRbRu6+2RgMkBubq43VcADhJeh9IpSzqoMUdj3Lvp1axdYHBGRxtbUTUPPAZeGf74UmNHExz944WUozStIoJxLehYEnUhEpFFFs/voE8AHwBFmVmhmVwB3ACeb2XLg5PDj5i17LB5KpJw4KuIS6D705KATiYg0qqg1Dbn7+bW8dFK0jhkVmXk8N+xePvvwJSacdR5HZOYFnUhEpFFpZHE9yisq+cPSNBb2vpwjcscFHUdEpNGpENTj1U82sHZ7Md8bkx10FBGRqFAhqMfD768io0MKJw3oFnQUEZGoUCGoQ/66HcxeuYVLRvcmFGdBxxERiQoVgjpMeX8VyQlxnJebGXQUEZGoUSGoxdbdpfx70RrOGtGL9DaJQccREYkaFYJaPDmvgOKySi1KLyKtngpBDSoqnUc++JKj+3Skf/e0oOOIiESVCkENXsvfwJpte7lMVwMiEgNUCGow5f1V9EpPYZy6jIpIDFAh2M9n63fy/orNXHR0b+JD+vOISOunM91+pnywiqT4OCaOUpdREYkNKgTVbN9TxrML1jBheE86pKrLqIjEBhWCaqbNK2BvWYW6jIpITFEhCKuodP714SrysjsyqGf7oOOIiDQZFYKwNz8tomDLXl0NiEjMCWTxejNbBewEKoByd88NIkd1Uz5YRY/2yYwfpC6jIhJbAikEYSe4+6YAj/+Vz4t2Mmv5Jn76zSNIUJdREYkxQRaC5qFgDitfnE5efC8mjtIKZCISe4IqBA7MNDMH/tfdJweSomAOPuXbnFBWyjcSEkjcOhraak1iEYktQbWDjHH3kcCpwLVmdtz+G5jZJDObZ2bzNm7cGJ0Uq2bh5aXEWyUJlMOqWdE5johIMxZIIXD3teHvRcCzwAEfw919srvnuntuly5dopOj97GUEk85cVgoEbLHRuU4IiLNWZMXAjNLNbN2+34GxgNLmzoHwJK4I7ig5OcsPfy/4NLnIFPNQiISe4K4R9ANeNbM9h3/cXd/OYAcTJ9fyMeh/vT5zjhITggigohI4Jq8ELj7F8Cwpj7u/orLKpixaC2nDO5OmoqAiMSwmO00/3p+Edv3lnHOkRlBRxERCVTMFoLp8wvo0T6ZY/p2DjqKiEigYrIQbNhRzNvLNnL2yAxCcRZ0HBGRQMVkIXh24RoqHc5Ws5CISOwVAndn+vxCcnt3IKdzatBxREQCF3OFYHHhdj4v2qWbxCIiYTFXCJ6aV0ByQhzfGtoj6CgiIs1CTBWC4rIKnlu8llMH96Cdxg6IiAAxVghe/WQDO4vL1SwkIlJNTBWCp+YX0is9hdF9OgUdRUSk2YiZQrB+ezHvLt/I2SN7EaexAyIiX4mZQvDMwkKNHRARqUFMFIJ9YwfysjvSu5PGDoiIVBcThWDB6m18sXE35+TqakBEZH8xUQimzy8kJSHEaUM0dkBEZH+tvhAUl1Xw/OK1nDqkO22TgliHR0SkeQukEJjZKWb2mZl9bmY3R/NYr3y8np0lGjsgIlKbINYsDgF/B04FBgLnm9nAaB1v+vxCMjqkcHSOxg6IiNQkiCuCPOBzd//C3UuBqcCEaBxo4yezGLLyQX542BaNHRARqUUQjea9gIJqjwuBoxr9KAVzSJ9+Nj8OlRGXPwMKsiAzr9EPIyLS0gVxRVDTR3M/YCOzSWY2z8zmbdy48eCPsmoWocoy4q2SuIoyWDWrAVFFRFq/IApBIZBZ7XEGsHb/jdx9srvnuntuly5dDv4o2WOJi08CC0EoEbLHNjiwiEhrFkTT0Fygn5nlAGuAicAFjX6UzDy49LmqK4HssWoWEhGpRZMXAncvN7MfAq8AIeAhd/84KgfLzFMBEBGpRyAjrNz9ReDFII4tIiJf1+pHFouISN1UCEREYpwKgYhIjFMhEBGJcSoEIiIxztwPGNTb7JjZRuDLBu7eGdjUiHGagjJHX0vLC8rcVFpa5rry9nb3ekfktohCcCjMbJ675wad42Aoc/S1tLygzE2lpWVujLxqGhIRiXEqBCIiMS4WCsHkoAM0gDJHX0vLC8rcVFpa5kPO2+rvEYiISN1i4YpARETq0KoLgZmdYmafmdnnZnZz0HnqY2YPmVmRmS0NOkskzCzTzN40s3wz+9jMrg86U33MLNnM5pjZ4nDm3wSdKVJmFjKzhWb2fNBZImFmq8xsiZktMrN5Qeepj5mlm9l0M/s0/P/06KAz1cXMjgj/bfd97TCzGxr0Xq21acjMQsAy4GSqFsOZC5zv7p8EGqwOZnYcsAv4l7sPDjpPfcysB9DD3ReYWTtgPnBmM/8bG5Dq7rvMLAF4F7je3T8MOFq9zOzHQC6Q5u6nB52nPma2Csh19xbRJ9/MpgCz3P0BM0sE2rj7tqBzRSJ8vlsDHOXuBz3mqjVfEeQBn7v7F+5eCkwFJgScqU7u/g6wJegckXL3de6+IPzzTiCfqjWpmy2vsiv8MCH81ew/DZlZBvAt4IGgs7RGZpYGHAc8CODupS2lCISdBKxoSBGA1l0IegEF1R4X0sxPUi2ZmWUDI4DZwSapX7iJZRFQBLzq7s0+M/BX4GdAZdBBDoIDM81svplNCjpMPfoAG4F/hpvfHjCz1KBDHYSJwBMN3bk1FwKr4blm/8mvJTKztsDTwA3uviPoPPVx9wp3H07Vetl5Ztasm+HM7HSgyN3nB53lII1x95HAqcC14abP5ioeGAnc6+4jgN1As7+vCBBuxjoDeKqh79GaC0EhkFntcQawNqAsrVa4nf1p4DF3fyboPAcjfOn/FnBKwFHqMwY4I9zmPhU40cweDTZS/dx9bfh7EfAsVc21zVUhUFjt6nA6VYWhJTgVWODuGxr6Bq25EMwF+plZTrhiTgSeCzhTqxK+8fogkO/ufw46TyTMrIuZpYd/TgHGAZ8Gm6pu7n6Lu2e4ezZV/x+/4e4XBRyrTmaWGu5AQLiJZTzQbHvDuft6oMDMjgg/dRLQbDs97Od8DqFZCAJas7gpuHu5mf0QeAUIAQ+5+8cBx6qTmT0BHA90NrNC4Nfu/mCwqeo0BrgYWBJucwf4eXhN6uaqBzAl3MsiDpjm7i2iO2YL0w14tuqzAvHA4+7+crCR6vVfwGPhD45fAN8LOE+9zKwNVT0jrzqk92mt3UdFRCQyrblpSEREIqBCICIS41QIRERinAqBiEiMUyEQEYlxKgQiIjFOhUCkBuEpia+p4/UUM3s7PB6htm1eM7MO0Uko0nhUCERqlg7UWgiAy4Fn3L2ijm0eqec9RJoFFQKRmt0B9A0v+HFnDa9fCMyAqnUZzOyd8LZLzWxseJvnqBr+L9KsaWSxSA3C02o/X9MCQeEpCFa7e/fw458Aye7+u3BTUZvw+gyY2XLgaHff3GThRQ5Sq51rSCSKOgPVFy2ZCzwUnon13+6+qNprRUBPQIVAmi01DYkcvL1A8r4H4ZXljqNqqcBHzOySatsmh7cXabZUCERqthNoV9ML7r4VCJlZMoCZ9aZq4Zj7qZqWe2T4eQO6A6uaIrBIQ6kQiNQg3Kb/Xvjmb003i2cCx4Z/Ph5YZGYLgbOBu8PPHwl86O7l0c4rcih0s1ikAcxsBPBjd7+4jm3uBp5z99ebLpnIwdMVgUgDuPtC4M26BpQBS1UEpCXQFYGISIzTFYGISIxTIRARiXEqBCIiMU6FQEQkxqkQiIjEuP8PZdb76ik9wp0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Instantiate the symbolic equations to the example.\n",
    "v0_num = 0.0\n",
    "v_t_example = v_t.subs(m,1576.0).subs(vd,40).subs(v0,v0_num).subs(k,1000).subs(d,0.5)\n",
    "dv_dt_example = dv_dt.subs(m,1576.0).subs(vd,40).subs(v0,v0_num).subs(k,1000).subs(d,0.5)\n",
    "\n",
    "# Define simulation step size, time range, and total number of steps to be taken.\n",
    "h=0.2\n",
    "t_range = np.arange(0,7,h)\n",
    "num_steps = len(t_range)\n",
    "\n",
    "# Compute analytical trajectory\n",
    "xs_c = [v_t_example.subs(t, x) for x in t_range]\n",
    "\n",
    "# Compute approximate trajectories\n",
    "xs_ee = Explicit_Euler().simulate(lambda v_val: dv_dt_example.subs(v,v_val), \n",
    "                                      v0_num, 0.0, h, num_steps)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t_range,xs_c, label='v')\n",
    "plt.plot(t_range,xs_ee, 'o', markersize=3, label='v_exEuler')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('speed (m/s)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Implement the Midpoint Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The midpoint method approximates the solution of a given IVP in the form of\n",
    "$$ \\dot{x} = f(x) $$\n",
    "by the following:\n",
    "$$ \\mathit{aux} = x(t) + f(x(t))\\frac{h}{2}$$\n",
    "$$ x(t+h) = x(t) + f(aux)h$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPXV+PHPyWSyEAghLGFJQgIiuwQMQWRxQ1xKwb2IioqCVlu1ra3a+qtYbR9ba1vrU7W4VOqGiFosKoIroAgECIJGQSCSsIU1rNkm5/dHLjwBskxCZu4kc96vV16Te+d77z0T5Z6531VUFWOMMeErwu0AjDHGuMsSgTHGhDlLBMYYE+YsERhjTJizRGCMMWHOEoExxoQ5SwTGGBPmLBEYY0yYs0RgjDFhLtLtAPzRrl07TUtLczsMY4xpUpYvX75TVdvXVa5JJIK0tDSys7PdDsMYY5oUEfnen3JWNWSMMWHOEoExxoQ5SwTGGBPmAt5GICIeIBvYrKpjRCQdmAEkAiuA61S1tL7nLSsro6CggOLi4sYNOEzFxMSQnJyM1+t1OxRjTJAFo7H4TiAXiHe2/wj8VVVniMjTwE3AU/U9aUFBAa1atSItLQ0Rabxow5CqsmvXLgoKCkhPT3c7HGNMkAW0akhEkoEfAM862wKcC8xyikwHLmnIuYuLi2nbtq0lgUYgIrRt29aerowJU4FuI/gb8CugwtluC+xV1XJnuwDoUt2BIjJFRLJFJHvHjh3VntySQOOxv6UxISh/KSx8rPI1gAKWCERkDFCoqsur7q6maLVrZarqNFXNVNXM9u3rHA9hjDHB4c/NuTHK5C+F6WPho99XvgYwGQSyjWAYMFZELgZiqGwj+BuQICKRzlNBMrAlgDEYY0yl/KWQtxDSRkBKVsPKHLk5+0rBEwXXv31COd20BP49DnylqMfL9z+Ywe7EARws8XGotJyDJT5aFK5g1LKb8VSU4Yvw8sIpf2ettzeHynwUl/o4XObjor2vMKG8BA8VldfLW1hz3CcpYIlAVe8D7gMQkbOBu1X1GhF5HbiCyp5D1wOzAxWDMaYZcOkGvvaCl9kWfxr7i8vZX1zOvuIy+qx/k2HOzdlXXsKs11/hlegyDhSXcbDEx8GScib63uBnnhIipQJfmfL6rFd50nfomGvd5pnD+ZFlREgFFb4yDq/7lEVR7YmN8hDrrfzZEDcQ36EZRFCOeKIqP1uAuDHFxD3ADBF5GFgJPOdCDCftnnvuoWvXrtx2220ATJ06lVatWvGLX/zC5ciMaUIa4QZ+fBmdOJv9HQax52Apew6VsedgKe1XvUOf8hIinBv4e2/PZFYLKDpcRtHhMvYdLmN8yevcFfF/N/C3/zOTJ30lx1zq9IgkBkdF4qWcciJZ6elHfEwkXRJiiIuKJC46ki5l56Nf/YeKinLweBl57iUM6TSYuCgPLaIiiYv20GpnPJ7X30Z9pURGRnHH9Tdyxwmf/0zI71t3EmwEQUkEqvoJ8Inz+wagUT/Rg//9iq+37GvMU9KnczwP/LBvje+PHz+eu+6662gimDlzJnPnzm3UGIxp0hrjJp+3EPWVIupDfaVsXjmPNUUp7DxQyq4Dpew8UEJm/muMcb6ll5eX8Ldpz/G/5TuPOc0gSeBl5wZeJpG8s687uypKaR3rpXPrWOJjvXQoPx/9pvIGLpFezht9OeclD6ZVjJf4GC+tYiJpEXUxUjAM8hYSmTaCR6q9OfeFrK6Qt5CItBGcUV2ZtiMqP29dN/mUrIAmgCOaxKRzoWjgwIEUFhayZcsWduzYQZs2bUhNTXU7LGOCozGqYjYuBOcmX+ErZc2iOXzaoQ2F+0so3F9M4f4S2u+N5fEKD16UMvVwx+I4Vny+4ug5Elp42RvTgwvxAmVUiJeOA87n/o69SWgRRWKcl4QWUbRpcTble7KI3voFsekjeKram2t/yO969HOdfjI358YqEyTNIhHU9s09kK644gpmzZrFtm3bGD9+vCsxGNPoGvsmX15KzoL/Mi8xnm1Fh9lSVMy2omI67otkuuf/bvJTv2zDCl1L61gvHVpF0yE+mrjuZ/I6/6B36SoOdR7KXSlDaNsyinYto0mMi8LrcTo+5g+CvIV40kZwbU0313bDoMew2j97CN2cg6lZJAK3jB8/nsmTJ7Nz504+/fRTt8Mxpm4BqK5Zt+RdPt3YjoI9h8jfc5iCPYdou8fL8+Lc5PHw8JpEVkdsICk+hs6tY8lISaBTwmg+quhKj8M5VHQdzuPdz6R9q2hivJ7jgsoArqr9c4XpDbyxWCI4CX379mX//v106dKFTp06uR2OMbXz8yaPrxScm3zx2k/I1R7k7TxY+bPrENHb4vmdevBq5Tf5e5e3ZoXm0io6kuTEFnRtG0fyKefxgaTSszgHSRvB0z2H0S4umoiI44cS9aaBkwuYRmSJ4CStXr3a7RCMqVTXt/0qN/mq/dLLfRVs2n2IdYUH2L8rnbFEEuFU11zzgZcV8z8HIEKgc0Is6e1O48UOT5DhW01F12H87pRhpLRpQXxs5HEj1PtiN/mmwRKBMU1BI1TplKcOI8LjhXLwSSSPr+3A/OwFbNx5kFLfkVlgonm/1VRGxa5lb9IQLk4dwm1t40hrF0dKYizRkUeqbYYA1i7WXFgiMCbU1bNKB18pxd99ypflp5C7dR+5W/fx9dZ9fLttP31993JGRC5faG927U7mlPaxnN2rPae0b0mPpFZ0bx9HqxibijzcWCIwxm0NrNI5Yu+hUjZG9KOfRBLh1NtPmBfJCl0MQJsWXvp0jue6M7rSu1N/enZsxR0dWlbTKGvClSUCY9zkz7f9tBGV7/lKUU8U38YM4LNFG/myYC+r8veSt6ty+oJBch8XtVpPUdIQzus2lJ92iqd3p3iS4qNtdllTK0sExripjm/724qKWbo7mW3pfyOq4DPe2dedZW+UAF/TMT6GASmtuWpwChnJCfRLHk28VeuYBrBEYEwg1VXtc9y3/S0Jp7No2SaWbNzNsrzd5O8+DECLqAQGpl5DVkYCk5MTGJCSQFJ8TJA/jGmuLBEYEyh1VPuoKnmxfck9/Z+UrlvA7KJ0Pn75ILCaxLgoBqe14fqhaWSlJ9KnUzyRnkCvI2XClSWCEJOXl0fv3r3p2bPn0X0///nPmThxYo3HTJ06lZYtW3L33XcHI0RzRAMaefckZvDZ+p0sWreThet2snnvYSCSLgk/YMipifxPeiKD0xLp3j7O6vVN0IRXIvBnXvMQ0L17d3JycgJ2/vLyciIjw+s/faPzs5FXPVGorxSfRPLr7HhmvTsfVWgVE8mZ3dty69ndGXFKO7q2bWE3fuOa8Lkb+PMPtx7qsx7Bo48+ysyZMykpKeHSSy/lwQcfZNmyZdx0000sXboUn89HVlYWr732Gi1btqzxmi1btuTAgQMAzJo1izlz5vDCCy8cU2b9+vXcfvvt7NixgxYtWvDMM8/Qq1cvbrjhBhITE1m5ciWDBg3isccea/BnN9TayFt0uIxPvi1k/teR7C37DaeVr2YZfaBDX+46rz3De7RjQHJrq+oxISNgiUBEYoAFQLRznVmq+oCIvACcBRQ5RW9Q1cB9/T2ijt4Z9eXvegTz5s1j3bp1LF26FFVl7NixLFiwgJEjRzJ27Fjuv/9+Dh8+zLXXXku/fv3Iy8tj/fr1ZGRkHD3HE088wYgR/q1ONGXKFJ5++ml69OjBkiVLuO222/joo48AWLt2LR988AEej/Ufr1M9GnnxRFHYdjDvfraR+bnbWbJhN+UVSruWUZzX/ywyel/Jj7u3tYFaJmQF8omgBDhXVQ+IiBdYJCLvOe/9UlVnBfDaJzruH+7JLvvm73oE8+bNY968eQwcOBCAAwcOsG7dOkaOHMlvf/tbBg8eTExMDH//+9+PHtPQqqEDBw7w+eefc+WVVx7dV1LyfyssXXnllZYE/OHP02NKFt+PmcGmle/z5u503vr3fuBrureP4+YR3Ti/TxIZKQl4TphkzZjQE8g1ixU44Gx6nR8N1PXqlJLl34pA9eDPegSqyn333cctt9xywnu7d+/mwIEDlJWVUVxcTFxcXK3Xq1qHXFxcfML7FRUVJCQk1JhE6jq/cdTy9Ji/+xBvr9rCf1dt4ZttB4mQ4WR2TeTXgzswqncS3drXXLVnTKgKaCWliHhEJAcoBOar6hLnrd+LyJci8lcRiQ5kDMdIyYIRv2i0huLx48czY8YMZs2axRVXXFFtmQsuuIDnn3/+aN3+5s2bKSwsBCqrcR566CGuueYa7rnnnjqvl5SURG5uLhUVFbz11lsnvB8fH096ejqvv/46UJmEVq1a1dCP13zlL4WFj1W+VufI06N4wBPFng5DeOGzjVz65GeM+NPHPPr+t8RFR/K7cX1Z+ptRzLx1KFNGdrckYJqsgDYWq6oPyBCRBOAtEekH3AdsA6KAaVQuZv+7448VkSnAFCBkl4D0Zz2C0aNHk5uby9ChQ4HKBt+XXnqJuXPnEhkZyYQJE/D5fJx55pl89NFHdOvW7YQ2gkmTJnHHHXfwyCOPMGbMGFJSUujXr9/R5FLVyy+/zI9//GMefvhhysrKGD9+PAMGDAjMH6Ap8rPa59DVb/Hdsvd4c1c6/36hiAotolfHVtxzYS9+OKATyW1auBO/MQEglTU4QbiQyAPAQVX9c5V9ZwN3q+qY2o7NzMzU7OzsY/bl5ubSu3fvQIQatsLib7rwMfjo95XVPuKBc39T+ZSI8wRVUMSMpZt4e9UWDpX6SE1swbiMzowd0JkeSa1cDt6Y+hGR5aqaWVe5QPYaag+UqepeEYkFRgF/FJFOqrpVKiu8LwHWBCoGY05QTaeBosNlzM7ZzCtLNvHNtv3Eej38cEAnfjQ4lUGpCda/3zR7gawa6gRMFxEPlW0RM1V1joh85CQJAXKAWwMYQ9CsXr2a66677ph90dHRLFmypIYjTEDU1e3T6TSgGxfyTcwAnlkcxburP6C4rIJ+XeJ5+JJ+jMvobF09TVgJZK+hL4GB1ew/N1DXdFP//v0DOhrY+MGP+v/DpT7e2JzE9OxM1hUeoGX0di4blMzVg1Ppn9zapcCNcVf4jCw2zV8t3T537C/hxcV5vPjF9+w5VEb/Lq350+Wn8YPTOhEXbf8MTHizfwGm+aim/n/d9v08t2gjb67cTGl5BaN6d2DyiG5kpSda3b8xDksEpmmprQ2gSv3/am9//voBfPztAqIjI7ji9GRuGp5Od+vrb8wJLBGYpsOP+f3fL0rliZwsvtqyj7ZxRfxs1Klce0YqbVsGb9yiMU2NJYIQ8fTTT9OiRYsT1h3Iy8tjzJgxrFljvWxragNQVT76ppC/zF/LV1v2kd4ujkcu688lA7vYAu3G+CGsEkFOYQ7Z27PJTMoko0NG3QcE0a23NotetIF1XBuAdh3OgrU7+Mv8tazK30tqYgv+fOUALsnobFM8G1MPYZMIcgpzmDxvMqW+UqI8UTwz+pmTSgb+rkfwySef8MADD5CUlEROTg6XXXYZ/fv35/HHH+fw4cP85z//oXv37sesMrZ8+XImTZpEixYtGD58+El97malysSBX0b253fv+Mj+fildEmJ55LL+XH56Ml5LAMbUW9j8q8nenk2pr5QKKiirKCN7e3bdB9Vi/PjxvPbaa0e3Z86cecz0z1WtWrWKxx9/nNWrV/Piiy+ydu1ali5dys0338wTTzxxQvkbb7yRv//97yxevPikYmxy6poMDljmO4XxuWcydnYZ+XsO8dC4vnx091mMz0q1JGBMA4XNE0FmUiZRnijKKsrwRnjJTKpz+o1a+bseAcDgwYOPTkrXvXt3Ro8eDVQOQvv444+PKVtUVMTevXs566yzALjuuut47733aPbqaAjO332Ih9/5mve/2k67ltH8dkwfJgxJtTYAYxpB2CSCjA4ZPDP6mUZtI/BnPQKonGriiIiIiKPbERERlJeXH1NWVcOzf3sNDcEHS8p58pPveGbhRjwi/OL8U7l5RDdioywBGNNYwiYRQGUyaMxG4vHjxzN58mR27tzJp59+2ijnTEhIoHXr1ixatIjhw4fz8ssvN8p5Q141DcH/WVnAI+99w/Z9JVyS0Zl7L+pNx9YxbkdqTLMTVomgsfmzHkFD/Otf/zraWHzBBRc02nlDWpWG4HWxGdwzp5wVm1ZxWnJrnrxmEKd3TXQ7QmOaraCtR3AybD2C4HD7b1q4v5g/zf2WWcsLaNcyml9d2JMrBiUTYev+GtMgrq9HYMwxapkaoqJCeWXpJh557xtKyn3cclY3fnLOKTYVtDFBYomgkdh6BLWopUfQpl2HuOeNL1m8YRfDTmnLw5f0J71dnMsBGxNemnQiCKUeNk19PYKAVhFW0yOoostgXvg8j0ff/5bICOGRy/rzo8EpIfPf05hwEsilKmOABUC0c51ZqvqAiKQDM4BEYAVwnaqW1vf8MTEx7Nq1i7Zt29rN4ySpKrt27SImJkA9co7rEVSQcDp3/XMx2d/v4eye7fnDpf3pnBAbmGsbY+oUyCeCEuBcVT0gIl5gkYi8B/wc+KuqzhCRp4GbgKfqe/Lk5GQKCgrYsWNH40YdpmJiYkhOTg7MyZ0eQb4NC/hvUTfuee0wMd5S/nLVAC4d2MUSuTEuC+RSlQoccDa9zo8C5wITnP3Tgak0IBF4vV7S09NPPlATFGujevPL1aWsKihidJ/2PHxJPzrE25gAY0JBQNsInIXrlwOnAP8A1gN7VfXIcNoCoEsgYzDuUlX+vfh7fv9OLi1jInni6oGMOa2TPQUYE0ICmghU1QdkiEgC8BZQXSf1alspRWQKMAWocQ4fEyJq6BpadLiMe9/4kvfWbOPcXh149IrTbIEYY0JQUHoNqepeEfkEOANIEJFI56kgGdhSwzHTgGlQOaAsGHGaBqiha+iq/L385NUVbN1bzK8v7sXNw7vZwDBjQlTA5u0VkfbOkwAiEguMAnKBj4ErnGLXA7MDFYMJguO6hurGhTy/aCNXPP05FRXw2i1DmTKyuyUBY0JYIJ8IOgHTnXaCCGCmqs4Rka+BGSLyMLASeC6AMZhAq9I1VD1e/vhNO57e8DWjeifx5ytPI6FFlNsRGmPqEMheQ18CA6vZvwHIOvEI0yQ5XUO35Mxj6upEPsprz/0/6MVNw9OtQdiYJqJJjyw27lNVnstrxyOLB5EUH8Prtw5kYGobt8MyxtSDJQLTYKXlFdz35mreWFHA6D5JPHrFAFq3sInijGlqLBGY2tXQNXTPwVJueWk5Szfu5q5RPbjzvB5WFWRME2WJwNSshq6h63cc4KYXlrGlqJjHx2cwLsPGBBrTlFkiMDWrZtbQz0u7ceuLy/F6Inh18hBbOcyYZsASganZcbOGfnC4B7c+t5T0dnE8f8NgUhJbuB2hMaYRWCIwNXO6hlZsXMhLW1P47UceRvRoyz+uGUS8rR5mTLNhicDU6lDSIO76OIJ5X2/n2jNSmfrDvkR6AjYg3RjjAksEpka7DpRww7+W8dWWIn47pg83DkuznkHGNEOWCEy1CvcVc82zS9i0+xDPTMzkvN5JbodkjAkQSwTmBFv2HuaaZ5ewfV8x/7pxMGd2b+d2SMaYALJEEM6qGSy2adchJjz7BUWHynjxpizrHmpMGLBEEK6qGSy2PqYP1zyzhOJyH69MPoP+ya3djtIYEwSWCMLVcYPFdqz+kB+t2Ieq8urkM+jdKd7tCI0xQWKJIFxVGSxWEeHl7mUt8Xjh5ZuHckqHlm5HZ4wJIksE4coZLLY5Zx73LG/Nxpi+zJw8hK5t49yOzBgTZIFcqjJFRD4WkVwR+UpE7nT2TxWRzSKS4/xcHKgYTO2WlHVn9LJM8uP6MfPWoZYEjAlTgXwiKAd+oaorRKQVsFxE5jvv/VVV/xzAa5s6rNi0hxv+tYwubWJ5+eYhJMXHuB2SMcYldSYCEekADAM6A4eBNUC2qlbUdpyqbgW2Or/vF5FcwOYrDgHrtu9n0gvL6BAfzSuTh9ChlSUBY8JZjVVDInKOiLwPvANcROVi9H2A+4HVIvKgiPjVtURE0qhcv3iJs+snIvKliDwvIrauYRBt3nuYic8vxeuJ4MVJlgSMMbU/EVwMTFbVTce/ISKRwBjgfOCN2i4gIi2dMnep6j4ReQp4CFDn9TFgUjXHTQGmAKSmpvr1YUwV1QwW232wlInPLeFASTkzbxlKalubRtoYA6KqgTu5iBeYA7yvqn+p5v00YI6q9qvtPJmZmZqdnR2QGJulagaLHewwiAnPLuGbrft48aYhZKXbiGFjmjsRWa6qmXWVq7PXkIjcKSLxUuk5EVkhIqP9OE6A54DcqklARDpVKXYplW0OpjEdN1isfMMCbn1pOWs2F/G/EwZZEjDGHMOf7qOTVHUfMBpoD9wIPOLHccOA64Bzj+sq+icRWS0iXwLnAD9rYOymJkcGi4kH9UTx+HdJLFy3k/+5rD/n97FZRI0xx/Kn++iRCegvBv6lqqvEj0npVXVRlWOrerce8ZmGcAaL6caFPFvQhSe+bMW9F/XiqswUtyMzxoQgfxLBchGZB6QD9zljAmrtOmpCQEoWT6xtw1++XMuUkd249azubkdkjAlRNSYCEYlU1XLgJiAD2KCqh0SkLZXVQyaEvfTF9/xl/louH5TMvRf2cjscY0wIq+2J4AsRKQDmAnNVdS+Aqu4CdgUjONMwi9bt5Lez13Burw48cnl/IiJseUljTM1qTASqmikiXakcTPY3EekCLALeAz5V1ZIgxWjqYdOuQ/zk1RX06NCKJ64eiNcWmjfG1KHWu4Sqfq+qT6vqJcCZwH+BUcBCEXknGAEa/x0sKWfKi9lUVCjTJp5OXLRNLmuMqZvfdwpVLQM+cn5wnhBMiFBVfjlrFWu37+eFG7NsJlFjjN/8GVA2RkRWisgeEdknIvtFZJ+qbg5GgMY/T36ynndXb+Pei3ox8tT2bodjjGlC/Hki+BtwGbBaAzkfhfHfcfMIffxNIX+e9y1jB3Rm8ohubkdnjGli/EkE+cAaSwIh4rh5hDaPm8Edb5TQu2M8f7z8NPwY62eMMcfwJxH8CnhXRD4FjvYUqm4SORMEVeYRUl8pc+fMwusZx7SJpxMb5XE7OmNME+RPIvg9cACIAaICG46pkzOPkPpKKSOS9w6cwj8mDSK5jU0pbYxpGH8SQaKq1jnbqAkSZx6hzz54i7+s7cDYMeMY2r2t21EZY5owf0YbfeDPtNMmeN7fl8q13w6n26Bzuf7MNLfDMcY0cf4kgtuBuSJyuGr30UAHZqqXv/sQd89cxYDk1jx8ST9rHDbGnLQ6q4ZUtVUwAjF1K/dV8LPXcgD43wmDiPFa47Ax5uTVtnh9Wm0HOiuWJTd2QKZmT32ynuzv9/DQJf1ISbTGYWNM46itauhREXlDRCaKSF8R6SAiqSJyrog8BHwG9K7pYBFJEZGPRSRXRL4SkTud/YkiMl9E1jmvbRr5MzVLOfl7+duH6xiX0ZlLBtrsHsaYxlNjIlDVK4H/B/QE/gEsBGYDNwPfAueq6vxazl0O/EJVewNnALeLSB/gXuBDVe0BfOhsm1ocLCnnzhkr6Rgfw+/G9XM7HGNMM1NrG4Gqfg38piEnVtWtwFbn9/0ikgt0AcYBZzvFpgOfAPc05Brh4sH/fkX+7kPMmDKU1rFet8MxxjQzQZms3mlvGAgsAZKcJHEkWXQIRgxN1XurtzIzu4Afn92drPREt8MxxjRDAZ+wXkRaAm8Ad6nqPn+7O4rIFGAKQGpqauACDDVVJpTbFn8a9765mtOSW3PXqFPdjswY00wFNBGIiJfKJPCyqr7p7N4uIp1UdauIdAIKqztWVacB0wAyMzPDY8K7KhPKqSeKp9r8gdLyzjw+3lYaM8YETm2L1w+q7UBVXVHb+1L51f85IPe4CereBq4HHnFeZ/sdbXNXdUK58lLitnzBA+PuJ72dLTJjjAmc2p4IHnNeY4BMYBUgwGlU1vUPr+Pcw4DrgNUikuPs+zWVCWCmiNwEbAKubFjozVCVCeVKKjxo2jB+NDjF7aiMMc1cbYvXnwMgIjOAKaq62tnuB9xd14lVdRGViaM659U/1DCQkkXJhLd46bWX+by8N3+ecLVNIWGMCTh/2gh6HUkCAKq6RkQyAhhTWPvD6lZML7qIF2/Kok2czfptjAk8fxJBrog8C7wEKHAtkBvQqMLU8u/3MH3x99xwZhojeti6w8aY4PAnEdwI/Bi409leADwVsIjCVJmvgl+/uZrOrWP45QU93Q7HGBNG/Jl9tFhEngbeVdVvgxBTWHp24Ua+3b6fZydmEhcd8OEdxhhzVJ2d00VkLJADzHW2M0Tk7UAHFk427TrE4x+u5cK+HRnVJ8ntcIwxYcafUUoPAFnAXgBVzQHSAhhTWFFV7p+9hsiICKaO7et2OMaYMORPIihX1aKARxKm5ny5lQVrd3D36FPp2DrG7XCMMWHIn8roNSIyAfCISA/gDuDzwIYVHooOl/Hgf7/mtOTWXDc0ze1wjDFhyp8ngp8CfYES4BWgCLgrkEGFiz/N/YbdB0v4w6X98UTYwDFjjDv86TV0CPiNiPxBVQ8GIabmq8rMossrevDykk3cPDydfl1aux2ZMSaM1ZkIRORM4FmgJZAqIgOAW1T1tkAH16wcN7Poi9G/o3Pr7vzsfJte2hjjLn+qhv4KXADsAlDVVcDIQAbVLB03s2jnvdn8blw/GzNgjHGdX5Pcq2r+cbt8AYileTsys6h4KFEPmjbcxgwYY0KCP19H853qIRWRKCp7DdlcQ/WVkoVOnM3rb8xg9t5uPHbVVW5HZIwxgH9PBLcCt1O58PxmIMPZNvU0Z08Kv9o+ivNHj7ExA8aYkOFPr6GdwDVBiKVZO1zq4+F3vqZ/FxszYIwJLf7MNdRNRP4rIjtEpFBEZotINz+Oe94pv6bKvqkisllEcpyfi0/2AzQVzy7cwPZ9JTzwwz42ZsAYE1L8qRp6BZgJdAI6A68Dr/px3AvAhdXs/6uqZjg/7/obaFO2Y38JT3+6ngv6JpGZluh2OMYYcwx/EoGo6ouqWu78HFmgplaqugDYfdIRNgOPf7iW4vIDX337AAAPWklEQVQKfnVhL7dDMcaYE/iTCD4WkXtFJE1EuorIr4B3RCRRRBry9fYnIvKlU3XUpgHHNynfFR7g1aX5TMhKpXv7lm6HY4wxJ/Cn++iPnNdbjts/icongzrbC6p4CnjIOe4h4DHnPCcQkSnAFIDU1NR6XCK0/HHuN8R6Pdw5qofboRhjTLX86TWU3lgXU9XtR34XkWeAObWUnQZMA8jMzKyzKioULd24m/lfb+fu0afSrmW02+EYY0y1/Ok1dKWItHJ+v19E3hSRgQ25mIh0qrJ5KbCmprJNnaryh3dz6Rgfw03D6/PQZIwxweVPG8H/U9X9IjKcyjmHpgNP13WQiLwKLAZ6ikiBiNwE/ElEVovIl8A5wM9OIvaQ9s7qreTk7+Xno08lNsrjdjjGGFMjf9oIjswr9APgKVWdLSJT6zpIVa+uZvdz9YitySop9/Gnud/Sq2MrLh+U7HY4xhhTK3+eCDaLyD+Bq4B3RSTaz+PC1ktfbGLT7kPce1EvGzxmjAl5/jwRXEXlwLA/q+pep57/l4ENqwlyFp050GkoT3x0gOGntOOsU9u7HZUxxtTJ3xXK3qyyvRXYGsigmpwqi85ESyTdiu/jvosnIWJPA8aY0GerojSGKovOSIVyU/Jm+na25SeNMU2D1fU3BmfRGR8RlBFJ1tnj3I7IGGP8ZomgMaRksf6iV3is7Epm9XuS9n1GuB2RMcb4zaqGGsnUnDjWRF/Bpz88x+1QjDGmXuyJoBGs3LSHhet2cstZ3YmP8bodjjHG1Islgkbwj4+/I6GFl2vP6Op2KMYYU2+WCE7SV1uK+CC3kEnD0mkZbTVtxpimxxLBSfrHx9/RKjqS689MczsUY4xpEEsEJ2Hd9v28t2YbE8/sSutYaxswxjRNlghOwpOfrCcm0sOkYY22ZIMxxgSdJYIG+n7XQWbnbOaaIam0tUVnjDFNmCWCBnrqk/VEeiKYMtIWnTHGNG2WCBpg897DvLGigPGDU+gQH+N2OMYYc1IClghE5HkRKRSRNVX2JYrIfBFZ57y2CdT1A2nap+tRhVvO6u52KMYYc9IC+UTwApXrGFR1L/ChqvYAPnS2m5TC/cW8uiyfywcl0yUh1u1wjDHmpAUsEajqAmD3cbvHUbnmMc7rJYG6fqA8u3Aj5b4Kfny2PQ0YY5qHYLcRJDkL2xxZ4KZDkK9/UnYfLOWlL75n7IDOpLWLczscY4xpFCE7J4KITAGmAKSmprobjLMM5dzCVA6VtuD2c05xNx5jjGlEwU4E20Wkk6puddY+LqypoKpOA6YBZGZmarACPIGzDKX6Srm0wkNB98fokdTKtXCMMaaxBbtq6G3geuf364HZQb5+/TnLUIr68FLOxM75bkdkjDGNKpDdR18FFgM9RaRARG4CHgHOF5F1wPnOdmhLG4F6oignAl+El46nne92RMYY06gCVjWkqlfX8NZ5gbpmQKRk8faAp/j2i/cYd+lV9EzJcjsiY4xpVDayuA7lvgr+uCaelV0n0TNzlNvhGGNMo7NEUIf5X29nS1ExNw5LczsUY4wJCEsEdXjh8zyS28RyXu8kt0MxxpiAsERQi9yt+1iycTcTh3bFEyFuh2OMMQFhiaAW0z/PI8YbwVWZKW6HYowxAWOJoAZ7Dpbyn5zNXDqwCwktotwOxxhjAsYSQQ1ey86nuKzCFqU3xjR7lgiq4atQXlz8PWd0S6RXx3i3wzHGmICyRFCND3K3s3nvYW6wpwFjTBiwRFCN6Z/n0SUhllHWZdQYEwYsERzn2237+Xz9Lq49oyuRHvvzGGOaP7vTHWf64jyiIyMYP9i6jBpjwoMlgiqKDpXx1orNjMvoTJs46zJqjAkPlgiqmJmdz+Eyn3UZNcaEFUsEDl+F8u8v8shKS6Rv59Zuh2OMMUFjicDx8TeF5O8+bE8Dxpiw48ri9SKSB+wHfEC5qma6EUdV0xfn0al1DKP7WpdRY0x4cSUROM5R1Z0uXv+o7wr3s3DdTn55QU+81mXUGBNm3EwEoSF/KRvfnUVWZBfGD7YVyIwx4cetRKDAPBFR4J+qOs2VKPKXotN/yDllpZzl9RK1Zyi0tDWJjTHhxa16kGGqOgi4CLhdREYeX0BEpohItohk79ixIzBR5C1Ey0uJlAq8lEPewsBcxxhjQpgriUBVtzivhcBbwAlfw1V1mqpmqmpm+/btAxNH1+GUEkk5EYgnCtJGBOQ6xhgTyoKeCEQkTkRaHfkdGA2sCXYcAKsjejKh5NesOfWncP3bkGLVQsaY8ONGG0ES8JaIHLn+K6o614U4mLW8gK88veh22SiI8boRgjHGuC7oiUBVNwADgn3d4xWX+Zids4UL+3Uk3pKAMSaMhW2n+Q9zCyk6XMYVpye7HYoxxrgqbBPBrOX5dGodw5nd27kdijHGuCosE8H2fcV8unYHlw9KxhMhbodjjDGuCstE8NbKzVQoXG7VQsYYE36JQFWZtbyAzK5tSG8X53Y4xhjjurBLBKsKiviu8IA1EhtjjCPsEsHr2fnEeCP4wWmd3A7FGGNCQlglguIyH2+v2sJF/TrRysYOGGMMEGaJYP7X29lfXG7VQsYYU0VYJYLXlxfQJSGWod3auh2KMcaEjLBJBNuKilm0bgeXD+pChI0dMMaYo8ImEby5ssDGDhhjTDXCIhEcGTuQlZZI17Y2dsAYY6oKi0SwYtNeNuw4yBWZ9jRgjDHHC4tEMGt5AbFeDxf3t7EDxhhzvGafCIrLfMxZtYWL+nekZbQb6/AYY0xocyURiMiFIvKtiHwnIvcG8lrvf7WN/SU2dsAYY2rixprFHuAfwEVAH+BqEekTqOvNWl5AcptYzki3sQPGGFMdz9SpU4N6wQcffPAM4DRVfWLq1Km+Bx98sA3Qa+rUqYtqOmbatGlTp0yZUu9rfbJ0OivXPsnZXeMY1m9otWVyCnOYs2EOERJBx7iOVqaJlwmlWKyMlQlGmdo8+OCDW6dOnTqtrnJuVJp3AfKrbBcAQxr7IjlrXuHurx6lrC3kHnyOUWuSyOg34dgyhTlMnjeZUl8pUZ4onhn9DBkdMqxMEy0TSrFYGSsTjDKNxY02guqG9eoJhUSmiEi2iGTv2LGj3hfJ3vA+ZQIVIpRL5fYJZbZnU+orpYIKyirKyN6ebWWacJlQisXKWJlglGksbiSCAiClynYysOX4Qqo6TVUzVTWzffv29b5IZrcLiFLwqOLVyu0TyiRlEuWJwiMevBFeMpMyrUwTLhNKsVgZKxOMMo1FVE/4Mh5QIhIJrAXOAzYDy4AJqvpVTcdkZmZqdnb9s2HOmlfI3vA+md0uOKFa6GiZwhyyt2eTmZRZ42OXlWk6ZUIpFitjZYJRpjYislxV68wgQU8EACJyMfA3wAM8r6q/r618QxOBMcaEM38TgSsjrFT1XeBdN65tjDHmWM1+ZLExxpjaWSIwxpgwZ4nAGGPCnCUCY4wJc5YIjDEmzLnSfbS+RGQH8H0DD28H7GzEcILBYg68phYvWMzB0tRiri3erqpa54jcJpEIToaIZPvTjzaUWMyB19TiBYs5WJpazI0Rr1UNGWNMmLNEYIwxYS4cEkGdc3GHIIs58JpavGAxB0tTi/mk4232bQTGGGNqFw5PBMYYY2rRrBOBiFwoIt+KyHcicq/b8dRFRJ4XkUIRWeN2LP4QkRQR+VhEckXkKxG50+2Y6iIiMSKyVERWOTE/6HZM/hIRj4isFJE5bsfiDxHJE5HVIpIjIiE/fbCIJIjILBH5xvl/uvr1bUOEiPR0/rZHfvaJyF0NOldzrRoSEQ+V6x6cT+ViOMuAq1X1a1cDq4WIjAQOAP9W1X5ux1MXEekEdFLVFSLSClgOXBLif2MB4lT1gIh4gUXAnar6hcuh1UlEfg5kAvGqOsbteOoiInlApqo2iT75IjIdWKiqz4pIFNBCVfe6HZc/nPvdZmCIqtZ7zFVzfiLIAr5T1Q2qWgrMAMa5HFOtVHUBsNvtOPylqltVdYXz+34gl8o1qUOWVjrgbHqdn5D/NiQiycAPgGfdjqU5EpF4YCTwHICqljaVJOA4D1jfkCQAzTsRdAHyq2wXEOI3qaZMRNKAgcASdyOpm1PFkgMUAvNVNeRjpnIhp18BFW4HUg8KzBOR5SIyxe1g6tAN2AH8y6l+e1ZE4twOqh7GA6829ODmnAikmn0h/82vKRKRlsAbwF2qus/teOqiqj5VzaByvewsEQnpajgRGQMUqupyt2Opp2GqOgi4CLjdqfoMVZHAIOApVR0IHARCvl0RwKnGGgu83tBzNOdEUACkVNlOBra4FEuz5dSzvwG8rKpvuh1PfTiP/p8AF7ocSl2GAWOdOvcZwLki8pK7IdVNVbc4r4XAW1RW14aqAqCgytPhLCoTQ1NwEbBCVbc39ATNOREsA3qISLqTMccDb7scU7PiNLw+B+Sq6l/cjscfItJeRBKc32OBUcA37kZVO1W9T1WTVTWNyv+PP1LVa10Oq1YiEud0IMCpYhkNhGxvOFXdBuSLSE9n13lAyHZ6OM7VnES1ELi0ZnEwqGq5iPwEeB/wAM+r6lcuh1UrEXkVOBtoJyIFwAOq+py7UdVqGHAdsNqpcwf4tbMmdajqBEx3ellEADNVtUl0x2xikoC3Kr8rEAm8oqpz3Q2pTj8FXna+OG4AbnQ5njqJSAsqe0beclLnaa7dR40xxvinOVcNGWOM8YMlAmOMCXOWCIwxJsxZIjDGmDBnicAYY8KcJQJjjAlzlgiMqYYzJfFttbwfKyKfOuMRairzgYi0CUyExjQeSwTGVC8BqDERAJOAN1XVV0uZF+s4hzEhwRKBMdV7BOjuLPjxaDXvXwPMhsp1GURkgVN2jYiMcMq8TeXwf2NCmo0sNqYazrTac6pbIMiZgmCTqnZ0tn8BxKjq752qohbO+gyIyDrgDFXdFbTgjamnZjvXkDEB1A6oumjJMuB5ZybW/6hqTpX3CoHOgCUCE7KsasiY+jsMxBzZcFaWG0nlUoEvisjEKmVjnPLGhCxLBMZUbz/Qqro3VHUP4BGRGAAR6UrlwjHPUDkt9yBnvwAdgbxgBGxMQ1kiMKYaTp3+Z07jb3WNxfOA4c7vZwM5IrISuBx43Nl/OvCFqpYHOl5jToY1FhvTACIyEPi5ql5XS5nHgbdV9cPgRWZM/dkTgTENoKorgY9rG1AGrLEkYJoCeyIwxpgwZ08ExhgT5iwRGGNMmLNEYIwxYc4SgTHGhDlLBMYYE+b+P9g1atOJPZvSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Midpoint:\n",
    "    def simulate(self, f, x0, t0, h, max_steps):\n",
    "        x = x0\n",
    "        xs = [x]\n",
    "        for i in range(1, max_steps):\n",
    "            # TODO: implement the following \n",
    "            x = x\n",
    "            xs.append(x)\n",
    "        return xs\n",
    "\n",
    "# Run solver and plot\n",
    "xs_mid = Midpoint().simulate(lambda v_val: dv_dt_example.subs(v,v_val), \n",
    "                                v0_num, 0.0, h, num_steps)\n",
    "plt.figure(1)\n",
    "plt.plot(t_range, xs_c, label='v')\n",
    "plt.plot(t_range, xs_ee, 'o', markersize=3, label='v_exEuler')\n",
    "plt.plot(t_range, xs_mid, 'o', markersize=3, label='v_mid')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('speed (m/s)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}