{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import necessary libraries and configure them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "\n",
    "from sympy import *\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg\n",
    "import numpy as np\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from collections import deque\n",
    "\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the Mass-spring-damper symbolic equations and print $\\dot{v_1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m1, c1, x1, d1, v1, Fc = symbols('m1 c1 x1 d1 v1 Fc', real=True)\n",
    "dx_dt = v1\n",
    "dv_dt = (1/m1)*(-c1*x1 - d1*v1 + Fc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABEAAAAMCAYAAACEJVa/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA9ElEQVQoFZWS4Q2CQAyFOcMAJG6gG5iwgW6AK+gGOoMj4AYGN3AEEjeAEcQN8HuEngfiD15S2r62D3pH1LZtBBLshOXy4szIK+xg+ZTvmmm6qAi2WBM2kkukGHEbuIdxsXNuBVFiwh57ddH3cSZMldKr4aNioDmPxBRhGmy8jr4usx55IK4yLtxdBXV40X5gIDolsmDIsCOoaXob0fvlKP9JQxHt+Aw7OIOM/BZyU3EoUtMgoQ4I6NpTvmwg3JcHzvU72mBBIDFdq87rKj8GL9D55dTXqsVhA6SueDbCdeYMa1WPwTqe/ROwhs5MP5/WUXzHyg/ZX7BOp2+vswAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle v_{1}$"
      ],
      "text/plain": [
       "v₁"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJsAAAAuCAYAAAA/ZmtKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGNklEQVR4Ae2dj3UUNxDGfTwXQKADpwMCFcR0AEkFwR3YJfCgA5IKgukAUgHBHeAOsN2B+X6ytEi70qG9u9U70Mx7Qv9Go9HMt5J2fRKr29vbg5hWq9Wp8k8Vjn35B8WXPv1A8X2Fxz4+U/vXvs6iPbeAfPtIKr5SwH9v5buTliqvxmALnUsxUHih+t9CWRyr/r3y56r/Oy639P5bwPv2uXz3rqW293Kd+SeAKma1Ep2rIsx4JR4r3zMLyLfxitVUuyzYpEFQiNmrRFeqIBj9WBZgi3SpWe2mtdolsP3pFfk/VkhPxVGcV9pmtpFBfoAsE8m6FWuxIWT3bH5NB/2/hp490I5r9mjiZRP6JbRV/EHtLqJ8k+S2eqg9L0MvvLJPFJ8p8MCx0X6o8EXj2tsXJO8zdGYWwx+A7JPCrP1aZAfw8Dkes+o+q+xVDS4OxJQENQT5vBy8ieuUZ492FJeN09Qr0PmzUKc0S/GnkG8Rq7+t9ZAMgIYRnX2UBnSM7ZQyn246rqBLTSz98OO1wuAzpfEFyt+vkRF4xO/soNjJDOXEImzCi2KCo1x+wqCGzEpI4Ql4owDIUPo6JyAuE09A+SDXtx/AF/Mvla7VQ3zMUO9zeqgcOwxO8fnBUcpj+MGRQYbKijIDz9KxdOBBwWcv4r6Ux5/ZB0TlWb1VzoPr/KeY9sxssX+fqWx4KKkTZWUdUjMijAj9roZMvwd+Gv2HdInEgyIo9jLmkYzncT5Oe7n/xWUV6b8ks7gk1+ghHowRvjGhc45ehvH7StqwHXA2UZzseyplTvpZwgbqxD0oit+OOsS3yeeOCr2vNNbQ5g+1T/yrPPb4qABO1ttVggaUkhbxzxi9PClrZyfV8yRlZ4lxH0vm5+ghXoyfjLWkm/iwi1tCSzyUi6plrpOzTZ10mPhCZfgQBdl3Jz6v0Zt2vv0w2/t2E5t43oldD1UxkJCJQCgg2WUkFPQmZa7C/+OfTgZTnHFi/qXSS+nhn1jUTmazpcaxjdw1NnC+lS83HUPpkwkvSlV0b8SFQGjd97U7juhfD0ZK2LNNSAZgXV+cdqkHOnvHoXdw1PAwqe40ql98bBt04Ja2qB2fs5z+fmyl7UPUJEnCP4yfGuQo+pd0DY3B5kCxIfr5s1UAq+sbZyiwl2v5PW5rPaQz4OLFKMz043Fh+JsI3MruB3mdmL0GMGk8vEmznwp+eCK+kK5VHP5YJisZchIArhO2UiWN2FA+VkAhiCXzSoG3jGqlNCjkQGGGY31v/h2qVg8PKj7xDN8TUV7l2IQXojA7AGAcBh9v6ex5sn8TLsmkTSvK6A/4bhTwD2Pib97Jclqjt3h4AMGD8+9cG2Q/6kpYF1Rj4LmGWELmXB024d+l3iVZ42V0Ez2tjVmgygK9g43lcte0hMxd65iTt0u987K07k6+ufzsZbI0G133NVwx357Yi0y+F82xg9rvXOac/jfl3aXe35PV9Z5NxjFqaIHel9GGprauDGyGgWYWMLA1M7V1ZGAzDDSzAH9BcD9VaNajddStBexttFvXtx+4LaPtbd5tjwa2bl3ffuAGtvY277ZHA1u3rm8/cANbe5t326OBrVvXtx+4ga29zbvt0cDWrevbD9zA1t7m3fZ42O3Io4H738xzgupIgYtYoHD8kAtlPurHia/Fx4EgjsRxSQvlHJZJDo6ozKhgAQPbnWGeCjRnAhPnZfkFLyf73akwlfET52vFcHIqyYFRecDIL3x/oSIm1QFKTqYlRwBjnh7T3YPNA4PjbRAzG+dB4+OHD1yNrshSeTyLUZ781t7L+t4dIl5cf1H3YJPLuYfuQkABOIBtPBsxS0HjC1Xgi8HHWQ4O7J5I1rFiglFkge5fEASQG28PB47R7EUVoMpdCwr/rGsqENYzdQ+2yPmAKneVAKBKLtXxMxczoSv3s2IkypI5CxjYvlmFu8eSZTFaWsczGHfO8bJwKR6WXlsyv9mxmDKwyTQeMMxUY1CVllbuRQnAPBHokpmvaO3OKwxsdwBgdmJfFgAUYEF57gIZXha4oYnLZvhUYlRhAftZeIWR5rIIhMyIk9uR5sr52fhtZlvGoyzJRiMLGNhGBtkmy95PgWWVe9BInyvwH88ZyQJfAe0YoP5eOx3GAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{Fc - c_{1} x_{1} - d_{1} v_{1}}{m_{1}}$"
      ],
      "text/plain": [
       "Fc - c₁⋅x₁ - d₁⋅v₁\n",
       "──────────────────\n",
       "        m₁        "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dv_dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace symbols by their values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx_dt_example = dx_dt\n",
    "dv_dt_example = dv_dt.subs(m1,1.0).subs(c1,1.0).subs(d1,1.0).subs(Fc,0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABEAAAAMCAYAAACEJVa/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA9ElEQVQoFZWS4Q2CQAyFOcMAJG6gG5iwgW6AK+gGOoMj4AYGN3AEEjeAEcQN8HuEngfiD15S2r62D3pH1LZtBBLshOXy4szIK+xg+ZTvmmm6qAi2WBM2kkukGHEbuIdxsXNuBVFiwh57ddH3cSZMldKr4aNioDmPxBRhGmy8jr4usx55IK4yLtxdBXV40X5gIDolsmDIsCOoaXob0fvlKP9JQxHt+Aw7OIOM/BZyU3EoUtMgoQ4I6NpTvmwg3JcHzvU72mBBIDFdq87rKj8GL9D55dTXqsVhA6SueDbCdeYMa1WPwTqe/ROwhs5MP5/WUXzHyg/ZX7BOp2+vswAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle v_{1}$"
      ],
      "text/plain": [
       "v₁"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx_dt_example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAARCAYAAAD+Mx8uAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADNElEQVRoBe2Zi1EbMRCGuYwLyIQOnA6cUEGSDkJSAdABjCtgSAckFTBOB0AFBDqADnDcgfN/svZG6KzLXXzmHmRn1pJW0mn170O6c7ZcLne6QFmWnUmPR/FC/FZ8Id3uVL54ejZscIamWFb7LP4tfl3nmRp/y1ybw3zxvXhssr6X7E/caWxGUnAjktdiuB/iB/FYTLsyaf6hBmP0nzZJ9YXktGfidybvW9k7bJqMOBkLw3LuVM4MGnspnsV6SEYk1XpW/IwutbWXzmPzSkq2TR8TCnB3gFL9q95h/6b2vhVsNj4mNrGFT6M8Yl7yHI6eHT+W6OJyea+o/4YcUh/3izPJvjvBAH7qYMN2A3xo7olPxGA3Ee+KH0PM1C5Q25nhjdfIPD1UMHaQqd8M94hpONDXP4UygTMRcwT1lSpj4x3B4eMxYt/wxLc5cr/GQMQYtZoZAuXKLp27UhoPv/Hj91XGjkIUEA1ECJFwRF3kssqq2tvfUmz8rgiO02CHZE/2bpkSPLjgO0phNFIHi137cVWLA3lcE98AYqOG61tk8O1hrvXsbeOL2uHGmUNmcc7i9TrSvjhvU2cucwrUMhaxPlWxYd6p9g0GRgTElclUXlkHZQqjkZ/Qyusba8sA6GeGp25kEfHgdSTqMS5y83gba+nQ2v9UtolFrHBVbJhn+ATPACeyZS1q+86AsnitGT5U3hwk9GruBblzBIO5IA2R6mDj9q+AIStAIW4ryV9+u+AMXAjfr9GTc+8u8nrOwSfHkzbP5ehizfwhiCphAwZiCyh3NAq3HCf1HQf9SVyadgZTqLAgyoiXYj495ySlSflzyTGqI8aqwmvkwUqS/9pXTifw4/bCjecju1fZCjbCAOPjNM4JVMZvVQTQIgoqiYo0KorqS6TQuWaR1k2ha8l+qX3rje3ONckwZn6rVd2IO8tU/SgO8WbwITay2icaMxPzpxbfFnhu7bORec9F0nXb2IAzl+ux1jpWydvWoV/XBZ7ZQPJSyjSwdEBfOwUGjnmu/XHc/Kc1CMQYNX1MrFmyNVEyLbemUfcWfoLR4DKDvJ2jhqODzECdFHqjDJF/vlb7RVMKoz/OqujLla6SYAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle - 1.0 v_{1} - 1.0 x_{1}$"
      ],
      "text/plain": [
       "-1.0⋅v₁ - 1.0⋅x₁"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dv_dt_example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above equation is a linear ODE, can can be put in the form $\\dot{x} = Ax$, where $x$ is a vector of position and velocity, and $A$ is a 2x2 matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.,  1.],\n",
       "        [-1., -1.]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_00 = diff(dx_dt_example, x1)\n",
    "A_01 = diff(dx_dt_example, v1)\n",
    "A_10 = diff(dv_dt_example, x1)\n",
    "A_11 = diff(dv_dt_example, v1)\n",
    "A_num = np.matrix([[A_00, A_01], [A_10, A_11]]).astype(np.float64)\n",
    "A_num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This matrix can be used for simulation, as shown below, or for computing the analytical solution of the system, as shown next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHGFJREFUeJzt3X2UVPWd5/H3Rx7SUTEaaJXYIjhhSIJhozQ+xBklJii6Htw9zsTnjHkiZ6MmHrMTHw9mzJkdNybZMTuOGQfZ+DQ+BDMj4eDKJqJms2Kg1YwgMjAsxFKzEEQnGFEavvtH3SZFUd11u7uq7q2qz+scDn1v/aj60qervv37fn+/exURmJmZpbFf1gGYmVnzcNIwM7PUnDTMzCw1Jw0zM0vNScPMzFJz0jAzs9ScNMzMLDUnDTMzS81Jw8zMUhuZdQC1Nm7cuJg4cWLWYZiZNZWenp7fRERntXEtlzQmTpzIypUrsw7DzKypSNqUZpzLU2ZmlpqThpmZpeakYWZmqbVcT8PMbOfOnRQKBXbs2JF1KLnT0dFBV1cXo0aNGtK/d9Iws5ZTKBQYM2YMEydORFLW4eRGRLB161YKhQKTJk0a0nNkWp6StEDSZkmr+nlckr4nab2kf5Z0XKNjNLPms2PHDsaOHeuEUUYSY8eOHdYMLOuexg+A2QM8fiYwOfkzF7i9nsH0bNrGbcvW07NpW8VjM2seThiVDff7kml5KiKekjRxgCHnAHdH8Z60yyUdLGl8RLxW61h6Nm3jovnLebd3N6NH7se8s6dy0+LVe47v+8KJTD/qkFq/rJlZU8l7T+MI4OWS40Jybq+kIWkuxZkIEyZMGNILLd+wlXd7d7M7YGfvbh5d9dpex8s3bN0z7sSjxzqBmFlbyro8VU2leVTscyLijojojojuzs6qu+ArOvHosYweuR8jBKNG7seZx4zf6/iQ/Udz0fzlfGfpWi6av9wlKzPLxKWXXsrChQsze/28zzQKwJElx13Aq/V4oelHHcJ9Xzhxr5nElMPH7Dkun4ks37DVsw2zFtKzaVtLVhJ6e3sZObJ2H/V5n2ksAj6TrKI6EXizHv2MPtOPOoTLPvHBPT8wpcflM5ETjx4LuFlu1gr6epq1qiSsWLGCadOmsWPHDt566y2mTp3KqlUVF4lyyy23MGPGDKZNm8aNN94IwMaNGznmmGP2jPn2t7/NN77xjX3j7unh1FNPZfr06Zxxxhm89lrx43HmzJlcd911nHrqqdx6663D+r+Uy3SmIel+YCYwTlIBuBEYBRAR3weWAGcB64HfAZ/NJtLKM5Hy5rmb5WbNqdaVhBkzZjBnzhxuuOEG3n77bS6++OK9kkCfpUuXsm7dOn7xi18QEcyZM4ennnoqVW92586dXHHFFTzyyCN0dnby4IMPcv3117NgwQIA3njjDZ588skh/x/6k/XqqQuqPB7AZQ0Kp6rpRx2y1w+SS1ZmraGvkrCzd/delYThmDdvHjNmzKCjo4Pvfe97FccsXbqUpUuXcuyxxwKwfft21q1blypprF27llWrVjFr1iwAdu3axfjx4/c8ft555w37/1BJ3nsauVaPHzQza7xKlYThev3119m+fTs7d+5kx44dHHDAAfuMiQiuvfZavvSlL+11vlAosHv37j3HlTbjRQRTp07l6aefrvj6lV6vFvLe08i1vh+0q06fsqc05R6HWXMq72kO19y5c/nmN7/JRRddxNVXX11xzBlnnMGCBQvYvn07AK+88gqbN2/msMMOY/PmzWzdupV33nmHxYsX7/Nvp0yZwpYtW/YkjZ07d7J69eqaxD4QzzSGqbRk5R6HmQHcfffdjBw5kgsvvJBdu3bx8Y9/nMcff5zTTjttr3Gnn346a9as4aSTTgLgwAMP5N577+XQQw9l3rx5nHDCCUyaNIkPfehD+7zG6NGjWbhwIV/5yld488036e3t5corr2Tq1Kl1/b+p2DZoHd3d3ZHVnftuW7ae7yxdy+6AEYKrTp/CZZ/4YCaxmLWzNWvW8OEPfzjrMHKr0vdHUk9EdFf7ty5P1VB/y3LNzFqFy1M1VI9mmpk1vxdeeIFLLrlkr3Pvec97eOaZZzKKaOicNGqsfFmumWUjInJzpduPfvSjPP/881mHARS/L8Ph8lQDeEWVWWN1dHSwdevWYX9Atpq+mzB1dHQM+Tk806gzr6gya7yuri4KhQJbtmzJOpTc6bvd61A5adSZd42bNd6oUaOGfDtTG5jLU3XmFVVm1ko806gzr6gys1bipNEAXlFlZq3C5SkzM0vNSSMDXoJrZs3K5akG8xJcM2tmnmk0WKUluGZmzcJJo8G8BNfMmpnLUw3mJbhm1sycNDLgJbhm1qxcnjIzs9ScNHLCy3DNrBm4PJUDXoZrZs3CM40c8DJcM2sWmSYNSbMlrZW0XtI1FR6fIGmZpOck/bOks7KIs968DNfMmkVm5SlJI4DbgFlAAVghaVFEvFgy7AbgoYi4XdJHgCXAxIYHW2dehmtmzSLLnsbxwPqI2AAg6QHgHKA0aQRwUPL1+4BXGxphA3kZrpk1gyzLU0cAL5ccF5Jzpb4BXCypQHGWcUWlJ5I0V9JKSSt9e0czs/rJMmmowrnyu8BfAPwgIrqAs4B7JO0Tc0TcERHdEdHd2dlZh1DNzAyyTRoF4MiS4y72LT99HngIICKeBjqAcQ2JLmPet2FmeZRlT2MFMFnSJOAV4HzgwrIxvwI+CfxA0ocpJo2Wrz9534aZ5VVmM42I6AUuBx4D1lBcJbVa0k2S5iTDvgZ8UdIvgfuBSyOivITVcrxvw8zyKtMd4RGxhGKDu/TcvJKvXwRObnRcWevbt7Gzd7f3bZhZrvgyIjnkfRtmlldOGjnlfRtmlke+9pSZmaXmpGFmZqk5aZiZWWpOGk3Cm/3MLA/cCG8C3uxnZnnhmUYT8GY/M8sLJ40m4Js0mVleuDzVBLzZz8zywkmjSXizn5nlgctTZmaWmpOGmZml5qRhZmapOWk0MW/4M7NGcyO8SXnDn5llwTONJuUNf2aWBSeNJuUNf2aWBZenmpQ3/JlZFpw0mpg3/JlZo7k8ZWZmqTlpmJlZak4aZmaWmpOGmZmllmnSkDRb0lpJ6yVd08+YT0t6UdJqSf/Q6BibiXeIm1m9ZbZ6StII4DZgFlAAVkhaFBEvloyZDFwLnBwR2yQdmk20+ecd4mbWCFnONI4H1kfEhoh4F3gAOKdszBeB2yJiG0BEbG5wjE3DO8TNrBGyTBpHAC+XHBeSc6X+EPhDST+XtFzS7EpPJGmupJWSVm7ZsqVO4eabd4ibWSNkublPFc5F2fFIYDIwE+gCfibpmIh4Y69/FHEHcAdAd3d3+XO0Be8QN7NGyDJpFIAjS467gFcrjFkeETuB/ytpLcUksqIxITYX7xA3s3rLsjy1ApgsaZKk0cD5wKKyMf8EfAJA0jiK5aoNDY3SzMz2SJU0JB0iaaqkoyXVJNFERC9wOfAYsAZ4KCJWS7pJ0pxk2GPAVkkvAsuAP48Id3jNzDKiiMotAEnvAy4DLgBGA1uADuAwYDnwtxGxrEFxptbd3R0rV67MOgwzs6YiqSciuquNG6insRC4G/jj8sazpOnAJZKOjog7hxeq1VPPpm1ujptZzfSbNCJi1gCP9QA9dYnIasYb/sys1lKtnpI0DZhYOj4iflSnmKxGKm34c9Iws+GomjQkLQCmAauB3cnpAJw0cq5vw9/O3t3e8GdmNZFmpnFiRHyk7pFYzXnDn5nVWpqk8bSkj5ReSNCahzf8mVktpUkad1FMHL8G3qF4+Y+IiGl1jczMzHInTdJYAFwCvMDvexpmZtaG0iSNX0VE+eU9zMysDaVJGi8ld8z7McXyFOAlt2Zm7ShN0ngvxWRxesk5L7ltUt4hbmbDUTVpRMRnGxGI1Z93iJvZcPV7xVpJN0h6/wCPnybp7PqEZfXgW8Ka2XANNNN4AfixpB3As/z+KreTgY8BPwH+S90jtJrxDnEzG65+L42+Z4A0GTgZGA+8TfHeF09FxNv1D2/wfGn0gbmnYWaV1OLS6ABExDpgXU2issx5h7iZDUeWt3s1M7Mm46RhZmapOWmYmVlqae6n0Ql8kX1vwvS5+oVljeLGuJkNRpod4Y8AP6O4xHZXfcOxRvJmPzMbrDRJY/+IuLrukVjD+XawZjZYaXoaiyWdVfdIrOH6NvuNEN7sZ2appNnc91vgAOBdYGdyOiLioDrHNiTe3Dc47mmYGdR2c9+Y2oS0L0mzgVuBEcD8iLi5n3F/AvwQmBERzgg15M1+ZjYYaXoaSJoDnJIcPhERi4f7wpJGALcBs4ACsELSovJ7kUsaA3wFeGa4r2lmZsNTtach6Wbgq8CLyZ+vJueG63hgfURsiIh3gQeAcyqM+ybwLWBHDV7TzMyGIU0j/CxgVkQsiIgFwOzk3HAdAbxcclxIzu0h6VjgyGozG0lzJa2UtHLLli01CM3MzCpJuyP84JKv31ej11aFc3u68pL2A/4b8LVqTxQRd0REd0R0d3Z21ig8MzMrl6an8VfAc5KWUfygPwW4tgavXQCOLDnuAl4tOR4DHAM8IQngcGCRpDluhteXV1SZWX/SrJ66X9ITwAyKSePqiPh1DV57BTBZ0iTgFeB84MKS130TGNd3nMTwn50w6su7xM1sIAPd7vVDyd/HUbwBU4FiD+IDyblhiYhe4HLgMYo3dnooIlZLuilZrWUZ8C1hzWwgA800rgLmAt+p8FgApw33xSNiCbCk7Ny8fsbOHO7rWXW+JayZDaTfpBERc5Mvz4yIvZa7Suqoa1SWmelHHcJ9XzjRPQ0zqyhNI/z/AOXlqErnrEV4l7iZ9affpCHpcIr7Jt6b7JfoWyJ7ELB/A2IzM7OcGWimcQZwKcWlsN8tOf9b4Lo6xmRmZjk1UE/jLuAuSedGxMMNjMnMzHJqoPLUxRFxLzBR0lXlj0fEdyv8M2tB3uxnZn0GKk8dkPx9YCMCsXzyZj8zKzVQeervkr//onHhWN74lrBmVirNpdG/JekgSaMk/VTSbyRd3IjgLHu+JayZlUqzT+P0iPi6pP9I8VIifwosA+6ta2SWC97sZ2al0iSNUcnfZwH3R8TryVVnrU14s5+Z9UmTNH4s6SXgbeDLkjrxXfTMzNpS1Z5GRFwDnAR0R8RO4C0q35bVzMxaXNWZhqRRwCXAKUlZ6kng+3WOy8zMcijN7V5vB6YDf5v8OS45Z22sZ9M2blu2np5N27IOxcwaKE1PY0ZE/LuS48cl/bJeAVn+ecOfWftKM9PYJekP+g4kHQ3sql9Ilne+u59Z+0oz0/hzYJmkDRQvj34U8Nm6RmW55rv7mbWvqkkjIn4qaTIwhWLSeCki3ql7ZJZb3vBn1r7SrJ7qAL4M/BHFe4P/TNL3y28Ba+3FG/7M2lOa8tTdFG+89N+T4wuAeyheTsTMzNpImqQxpWz11DKvnjIza09pVk89J+nEvgNJJwA/r19I1oy8b8OsPaSZaZwAfEbSr5LjCcAaSS8AERHT6hadNQXv2zBrH2mSxux6vbik2cCtwAhgfkTcXPb4VcAXgF5gC/C5iNhUr3hsaHyjJrP2kWbJbV0+pCWNAG4DZlG8T8cKSYsi4sWSYc9RvFDi7yT9J+BbwHn1iMeGzvs2zNpHmplGvRwPrI+IDQCSHqB49dw9SSMilpWMXw74joE55H0bZu0jy6RxBPByyXGBYv+kP58HHq30gKS5wFyACRMm1Co+GwTv2zBrD2lWT9VLpdv/RcWBxXuSdwO3VHo8Iu6IiO6I6O7s7KxhiGZmVirLmUYBOLLkuAt4tXyQpE8B1wOn+vIlZmbZynKmsQKYLGmSpNHA+cCi0gGSjgX+DpgTEZsziNHMzEpkljQiohe4HHgMWAM8FBGrJd0kaU4y7BbgQOCHkp6XtKifp7Oc8WY/s9aUZXmKiFgCLCk7N6/k6081PCgbNm/2M2tdWZanrEX5Jk1mrctJw2qub7PfCOHNfmYtJtPylLUmb/Yza11OGlYX3uxn1ppcnjIzs9ScNKxhvAzXrPm5PGUN4WW4Zq3BMw1rCC/DNWsNThrWEF6Ga9YaXJ6yhvAyXLPW4KRhDeNluGbNz+UpMzNLzUnDMuMluGbNx+Upy4SX4Jo1J880LBNegmvWnJw0LBNegmvWnFyeskx4Ca5Zc3LSsMx4Ca5Z83F5yszMUnPSsFzxMlyzfVV6X5Sfa9R7x+Upyw0vw7V20LNp2169vPLj8jHAPu+L8nPzzp7KTYtXN+S946RhuVFpGa6ThjWTwSaASh/25WPOPa6r4vL00nOPrnqtYe8dJw3Ljb5luDt7d3sZrjWFwc4IyhNApQ972DshBFR8X5SeO/OY8azY+HpD3jtOGpYbXoZreVaprDTYGUF5Aujvw750zLnHdXHucV37vC/K3ytTDh/TkPdOpklD0mzgVmAEMD8ibi57/D3A3cB0YCtwXkRsbHSc1jhehmt5UW0WUV5OTTMjqJQAKn3YV/rlqfx9Uf5eadR7J7OkIWkEcBswCygAKyQtiogXS4Z9HtgWER+UdD7wX4HzGh+tZaVSjdisHgZKEpVmEeXl1LQzAqDqh32ef3nKcqZxPLA+IjYASHoAOAcoTRrnAN9Ivl4I/I0kRUQ0MlDLhldTWb0MttRUaRbRXzm1mRLAUGSZNI4AXi45LgAn9DcmInolvQmMBX5TOkjSXGAuwIQJE+oVrzWYV1NZrQy31NTfLKLVEkIaWSYNVThXPoNIM4aIuAO4A6C7u9uzkBbh1VRWC2ka1mlLTe2WICrJMmkUgCNLjruAV/sZU5A0Engf8HpjwrOseTWVDUV56SlNwzptqcmyTRorgMmSJgGvAOcDF5aNWQT8GfA08CfA4+5ntJd2nP7b4FQrPQ1mFuGfteoySxpJj+Jy4DGKS24XRMRqSTcBKyNiEXAncI+k9RRnGOdnFa/lh1dUWZ80pafLPvFBzyJqKNN9GhGxBFhSdm5eydc7gD9tdFyWX15R1d6GUnoCzyJqyTvCral4RVV7qVXpyWrHScOaildUtQ+XnvLJScOaildUtbbSmYVLT/nkpGFNp/xDwo3x1lA+s5h39lSXnnLIScOamhvjzataU3vb79516SmHnDSsqbkx3pwqJftK/SqXnvLHScOamhvjzWOgfsVATW3LFycNa2qVGuPuceRPtX6Fm9rNw0nDml7pB417HPkw1H6F5Z+ThrUU9ziy535Fa3PSsJbiHkc23K9oH04a1lL62/znPkdtDXR5D/crWpuThrWcSpv/3OeonWqX93C/orU5aVjLc59jeIZ6UyN/j1uTk4a1PPc5hi5NU9uX92gvThrW8ryXY3CG2tT297E9OGlYW/BejnS8Cc+qcdKwtuMex+95E54NlpOGtZ3+ehztVrLyJjwbCicNazv99TjaoWTlTXg2XE4a1pbKf3uu9AHad76ZPzy9Cc9qzUnDjH1LVofsP7rpZh7l5TVvwrN6cNIwY9+SVbPNPCqV17wJz+rBScMsUf4BmmbmkWXzvFp/wpvwrB4ySRqS3g88CEwENgKfjohtZWM+BtwOHATsAv4yIh5sbKTWrtLOPMoTCdRmNlKp1DRQ6alSf6K/izc6WdhwZDXTuAb4aUTcLOma5PjqsjG/Az4TEeskfQDokfRYRLzR6GCtPQ0086iUSB5+tsCPni1UnY0MJSHctHj1gKWn/voTLj9ZrWWVNM4BZiZf3wU8QVnSiIh/Kfn6VUmbgU7AScMarr/f2ksTiaDqbKQ8AaRJCI+ueq1q6cn9CWuUrJLGYRHxGkBEvCbp0IEGSzoeGA38ayOCM6uk/EO5PJEAPPxsYcDZSHkCSJMQzjxmPCs2vp6q9GRWb3VLGpJ+Ahxe4aHrB/k844F7gD+LiN39jJkLzAWYMGHCICM1G7ryRFJtNlKeANImhCmHj3HpyXJBEdH4F5XWAjOTWcZ44ImImFJh3EEUS1d/FRE/TPPc3d3dsXLlyprGazYcg+1pmGVBUk9EdFcdl1HSuAXYWtIIf39EfL1szGjgUeDHEfHXaZ/bScPMbPDSJo39GhFMBTcDsyStA2Ylx0jqljQ/GfNp4BTgUknPJ38+lk24ZmYGGc006skzDTOzwcv7TMPMzJqQk4aZmaXmpGFmZqk5aZiZWWot1wiXtAXYNIynGAf8pkbhNILjrS/HW1+Ot74GE+9REdFZbVDLJY3hkrQyzQqCvHC89eV468vx1lc94nV5yszMUnPSMDOz1Jw09nVH1gEMkuOtL8dbX463vmoer3saZmaWmmcaZmaWmpNGQtJsSWslrU+uvJtrkhZI2ixpVdaxVCPpSEnLJK2RtFrSV7OOqRpJHZJ+IemXScx/kXVM1UgaIek5SYuzjiUNSRslvZBcjDT3F4yTdLCkhZJeSn6WT8o6pv5ImlJyodfnJf2bpCtr8twuTxXfbMC/ULzibgFYAVwQES9mGtgAJJ0CbAfujohjso5nIMk9U8ZHxLOSxgA9wH/I+fdXwAERsV3SKOB/A1+NiOUZh9YvSVcB3cBBEXF21vFUI2kj0B0RTbHvQdJdwM8iYn5y64b9IyL3t59OPt9eAU6IiOHsYQM80+hzPLA+IjZExLvAAxTvY55bEfEU8HrWcaQREa9FxLPJ178F1gBHZBvVwKJoe3I4KvmT29+wJHUB/x6YX22sDV5yQ7hTgDsBIuLdZkgYiU8C/1qLhAFOGn2OAF4uOS6Q8w+1ZiVpInAs8Ey2kVSXlHueBzYD/ysi8hzzXwNfByreEjmnAlgqqSe5ZXOeHQ1sAf5HUgKcL+mArINK6Xzg/lo9mZNGkSqcy+1vlc1K0oHAw8CVEfFvWcdTTUTsioiPAV3A8ZJyWQaUdDawOSJ6so5lkE6OiOOAM4HLkpJrXo0EjgNuj4hjgbeAZuh9jgbmAKlul52Gk0ZRATiy5LgLeDWjWFpS0hd4GLgvIn6UdTyDkZQhngBmZxxKf04G5iQ9ggeA0yTdm21I1UXEq8nfm4F/pFgmzqsCUCiZbS6kmETy7kzg2Yj4f7V6QieNohXAZEmTksx8PrAo45haRtJUvhNYExHfzTqeNCR1Sjo4+fq9wKeAl7KNqrKIuDYiuiJiIsWf3ccj4uKMwxqQpAOSRREkZZ7TgdyuBIyIXwMvS5qSnPokkNuFHCUuoIalKShOudpeRPRKuhx4DBgBLIiI1RmHNSBJ9wMzgXGSCsCNEXFntlH162TgEuCFpEcAcF1ELMkwpmrGA3clK0/2Ax6KiKZYytokDgP+sfj7BCOBf4iI/5ltSFVdAdyX/GK5AfhsxvEMSNL+FFeEfqmmz+slt2ZmlpbLU2ZmlpqThpmZpeakYWZmqTlpmJlZak4aZmaWmpOGmZml5qRhVgPJZbO/PMDj75X0ZLLvo78xP5F0SH0iNKsNJw2z2jgY6DdpAJ8DfhQRuwYYc0+V5zDLnJOGWW3cDPxBcsObWyo8fhHwCBTvLyLpqWTsKkl/nIxZRPGyD2a55R3hZjWQXPJ9caUbYiWXnfhVRByeHH8N6IiIv0zKVfsn9xlB0jrgxIjY2rDgzQbB154yq79xQOkNe1YAC5Ir//5TRDxf8thm4AOAk4blkstTZvX3NtDRd5DcdfEUirfgvEfSZ0rGdiTjzXLJScOsNn4LjKn0QERsA0ZI6gCQdBTFmyb9PcVLxh+XnBdwOLCxEQGbDYWThlkNJD2InyeN7UqN8KXAHyVfzwSel/QccC5wa3J+OrA8InrrHa/ZULkRbtYAko4FroqISwYYcyuwKCJ+2rjIzAbHMw2zBoiI54BlA23uA1Y5YVjeeaZhZmapeaZhZmapOWmYmVlqThpmZpaak4aZmaXmpGFmZqn9f+5XVKLJcDYPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Explicit Euler solver\n",
    "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, t0 + h*(i-1))*h\n",
    "            xs.append(x)\n",
    "        return xs\n",
    "\n",
    "# Simulation step size, time range, and total number of steps to be taken\n",
    "h = 0.1\n",
    "t_range = np.arange(0.0,7.0,h)\n",
    "num_steps = np.size(t_range)\n",
    "\n",
    "# Initial condition\n",
    "x0 = np.matrix([[1.0],[0.0]]) \n",
    "\n",
    "# Simulation\n",
    "xs_e = Explicit_Euler().simulate(lambda x,t: A_num*x, x0, 0.0, h, num_steps)\n",
    "\n",
    "# Plotting\n",
    "plt.figure(1)\n",
    "plt.plot(t_range,[x[0,0] for x in xs_e], 'o', markersize=3, label='x_euler')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('position (m)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytical solution of a linear IVP $\\dot{x} = Ax$, where $x_0$ is known, is given by \n",
    "$$x(t) = e^{tA}x_0$$\n",
    "where $e^{At}$ is the matrix exponencial.\n",
    "Given a simulation step size $h$, the above formula can be put in the iterative form\n",
    "$$x(t+h) = e^{hA}x(t)$$\n",
    "where $x(t)$ is given.\n",
    "The iterative form is more efficient for computation, as we only compute $e^{hA}$ once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_solution_matrix(A, x0, h, max_steps):\n",
    "    x = x0\n",
    "    exphA = linalg.expm(h*A) # Computes matrix exponencial.\n",
    "    xs = [x]\n",
    "    for i in range(1, max_steps):\n",
    "        x = exphA * x\n",
    "        xs.append(x)\n",
    "    return xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the analytical solution against the simulation computed above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VOX1wPHvyUxCwh5CWCSQgCCyCiRAEGRRQUAEW62oqFVBarUuPzdsa926aF1rrdoqoqKIWkVBigqKgAtREgTZBVkkgAIh7FuW8/vjDhjCJBnIzNyZ5Hye5z6Ze+fNnZNA5sy7i6pijDHGBCLG7QCMMcZED0saxhhjAmZJwxhjTMAsaRhjjAmYJQ1jjDEBs6RhjDEmYJY0jDHGBMyShjHGmIBZ0jDGGBMwr9sBBFvDhg01LS3N7TCMMSaq5OTkbFfV5IrKVbmkkZaWRnZ2ttthGGNMVBGRDYGUs+YpY4wxAbOkYYwxJmCWNIwxxgSsyvVpGGOqjoKCAnJzczl48KDboVQZ8fHxpKSkEBsbe1Lfb0nDGBOxcnNzqVOnDmlpaYiI2+FEPVUlLy+P3NxcWrZseVL3cLV5SkQmiMhWEVlaxvMiIv8UkTUi8q2IdAt3jMYY9xw8eJCkpCRLGEEiIiQlJVWq5uZ2n8bLwOBynh8CtPEdY4HnQhlMzoZ8nvl0DTkb8v2eG2PCzxJGcFX29+lq85SqzhORtHKKjAAmqrMnbZaI1BeRpqq6Jdix5KzPY9GEm1lXlMKsmDQuHTKQ+z9YzeHCYuK8MUwak0l6amKwX9YYY6JKpPdpNAM2ljjP9V07JmmIyFicmggtWrQ4qRdaumIFo+Qj4mMLACia+Qe6yil86enAC4UXkLU2D4CstXlktkqyBGKMqZYiPWn4q0fpcRdUnweeB8jIyDju+UB0bN+Bbl++QrOiTXTw5jKmzX62r8nmCs/HXOaZzboNl3Ht7D78VFjbah7GmEq7+uqrGTZsGBdffLHboZyQSE8auUDzEucpwOZQvFB6aiKvjjnzaE2iY2oiORvymbzsW4bueIW2a15jVsxbTPAM5rnCC8lam2dJw5gIlLMhv0q2CBQWFuL1uv+W7XZHeEWmAVf5RlFlArtC0Z9xRHpqIjcOaH30P1p6aiJXDe1HwysmsPwXM5mnXbjZ+x5vxP2Zs5oUAtZZbkwkydmQz6jxWTw+cxWjxmdV+u9ywYIFdO7cmYMHD7Jv3z46dOjA0qV+B3vy6KOP0r17dzp37sx9990HwPr16+nYsePRMo899hj333//8XHn5NCvXz/S09M577zz2LLFeZvr378/f/jDH+jXrx9PPfVUpX6WYHE1bYnIZKA/0FBEcoH7gFgAVf03MAMYCqwB9gPXuBMpdDyjO4fqv8GMr//LeavuwzPjQlbs+zejph6wznJjIkTW2jwOFxZTrFBQWFzpFoHu3bszfPhw7rnnHg4cOMAVV1xxTBI4YubMmaxevZqvv/4aVWX48OHMmzcvoD7WgoICbrrpJqZOnUpycjJvvvkmf/zjH5kwYQIAO3fuZO7cuSf9MwSb26OnLqvgeQVuDFM4FUpPTYTUsfBTb5h8KW3+9yuGFI/hXe0TlP+gxpjKyWyVRJw3hoLCYmK9MWS2Sqr0Pe+99166d+9OfHw8//znP/2WmTlzJjNnzqRr164A7N27l9WrVweUNFatWsXSpUsZOHAgAEVFRTRt2vTo8yNHjqz0zxBM7jeQRaPGHeC6ORx49XKe/PFZUuUn/i2/Csp/UGPMyUtPTWTSmMyg9mns2LGDvXv3UlBQwMGDB6lVq9ZxZVSV3//+9/zmN7855npubi7FxcVHz/1NqlNVOnTowPz58/2+vr/Xc1Ok92lErlpJ1LluOttbX8yt3neY1Wsp6b7Oc+vjMMY9pfsmK2vs2LH8+c9/ZtSoUYwbN85vmfPOO48JEyawd+9eADZt2sTWrVtp3LgxW7duJS8vj0OHDjF9+vTjvrdt27Zs27btaNIoKChg2bJlQYk9FKymURmeWBpe/jz89zDNv/4L6+LqM2puivVxGFNFTJw4Ea/Xy+WXX05RURFnnnkms2fP5uyzzz6m3KBBg1ixYgW9evUCoHbt2rz22ms0atSIe++9l549e9KyZUtOP/30414jLi6Ot99+m5tvvpldu3ZRWFjIrbfeSocOHcLyM54ocboNqo6MjAwN+859hYfg9UsoXvsZvym4lVlF6XgEbhvUlhsHtA5vLMZUIStWrKBdu3Zuh1Hl+Pu9ikiOqmZU9L3WPBUM3howchIHGnbkX95/0itmRdA64YwxJpJY81Sw1KhNrWvf48B/BvLy3idYPeJ9OlrTlDFVypIlS7jyyiuPuVajRg2++uorlyIKP0sawVSzAQnXvAf/6UvHz2+C0z+GuJpuR2WMCZJOnTqxaNEit8NwlTVPBVv95nDRC7B1Ocy4A7BZ48aYqsNqGqHQ+lzoeyfMe4T1tc5g1LwWNqLKGFMlWE0jVPrfDS37kTL/T7QqWnfMsgbGGBOtLGmESowHLnqR4vj6PBv7FPVkv42oMsZEPUsaoVQ7mbiRr5Aas403W7xrTVPGGADS0tLYvn17uWVefvllNm/+eSeIMWPGsHz58lCHViFLGqGW2gs56zZO/+l/pB+sPsPyjDGVUzppjB8/nvbt27sYkcM6wsOh752w8n8w/VZoMR8SrLZhzAn74G74cUlw79mkEwx5uMynFyxYwOjRo/n6668pKiqiR48evPnmm8ctj75lyxZGjhzJ7t27KSws5LnnnuOss85i8uTJ/O1vf0NVOf/88/n73/9+zPetX7+eYcOGHd2j47HHHmPv3r107NiR7OxsRo0aRUJCAvPnz2fIkCE89thjZGRklHnf2rVrc8sttzB9+nQSEhKYOnUqjRs3DuqvzGoa4eCtASOegb1b4aM/2hBcY6JEyf007rrrrjL303j99dc577zzWLRoEYsXL6ZLly5s3ryZcePGMXv2bBYtWsSCBQt47733Anrdiy++mIyMDCZNmsSiRYtISEg4+lx59923bx+ZmZksXryYvn378sILLwTnF1GC1TTCpVk36HMrfPY4/8lJ4ePCM2wIrjEnopwaQSgFsp9G9+7dufbaaykoKODCCy+kS5cuzJ49m/79+5OcnAzAqFGjmDdvHhdeeGGl4lmwYEGZ942Li2PYsGEApKenM2vWrEq9lj9W0winfuPIq3kqD8S8QG3dZ0NwjYkCR/bT2LNnj9/9MAD69u3LvHnzaNasGVdeeSUTJ04kkMVgvV5vhfttlFbefWNjYxERADweD4WFhRXe70RZ0ggnbw22nvMkyezknthJNgTXmCgQyH4aGzZsoFGjRlx33XWMHj2ahQsX0rNnT+bOncv27dspKipi8uTJ9OvX75jvK2+/jTp16rBnz57jXiuQ+4aSNU+FWbv0fvy4fiyXLPk3nYfeyOnWNGVMxAp0P405c+bw6KOPEhsbS+3atZk4cSJNmzbloYceYsCAAagqQ4cOZcSIEcd8X2xsbJn7bVx99dVcf/31RzvCjwjkvqFk+2m44fA++FcPSKgPY+eCx3K3Mf7YfhqhYftpRJu4WjD4b/DTUlgw3u1ojDEmYPYR1y3thkOrAfDpX6HjL8nJiyVrbR6ZrZJsNJUxEcr207Ck4R4RGPooPNuLvHfHMeq7S2wlXGP8UNWjI4LcVhX206hsl4Q1T7mpYRs483ckfT+FzkXLbSVcY0qJj48nLy+v0m90xqGq5OXlER8ff9L3cLWmISKDgacADzBeVR8u9XwL4BWgvq/M3ao6I+yBhlLfOzn8zZs8qC8z/PBfifHG2jBcY3xSUlLIzc1l27ZtbodSZcTHx5OSknLS3+9a0hARD/AMMBDIBRaIyDRVLbmM4z3AW6r6nIi0B2YAaWEPNpTiahE39CFO/++veaHDUmqf9VtrmjLGJzY2lpYtW7odhinBzeapHsAaVV2rqoeBN4DSg40VqOt7XA/YTFXUfgSk9qbf5hdJb+xxOxpjjCmTm0mjGbCxxHmu71pJ9wNXiEguTi3jJn83EpGxIpItItlRWY0VgUF/gf3b4fMn3Y7GGGPK5GbS8DcconRv12XAy6qaAgwFXhWR42JW1edVNUNVM44s4hV1mnWDTpdA1rOwc2PF5Y0xxgVuJo1coHmJ8xSOb34aDbwFoKrzgXigYViic8M594IqzP6zLZ9ujIlIbo6eWgC0EZGWwCbgUuDyUmV+AM4BXhaRdjhJIwrbnwJUvzn0ugE+f5KHvunMwsKWNm/DGBNRXKtpqGoh8DvgI2AFziipZSLyoIgM9xW7HbhORBYDk4GrtaoP2O5zG/tjE7lDXqNY1eZtGGMiiqvzNHxzLmaUunZvicfLgd7hjstV8XXZnn4bmVl/4jxPDnNjeti8DWNMxLAZ4RGoxcAbOFDvVB6qN4VJ12ZY05QxJmJY0ohEHi8J591Pg/3rSd/5kdvRGGPMUZY0IlW7C+CUbjDnYSg85HY0xhgDWNKIXCLOENxdGyF7gtvRGGMMYEkjsp06AFr2hXmPwaHj9wo2xphws6QR6c65D/ZvZ9OHT9hkP2OM6yxpRLqUDPJbDKLuwud4cWY2o8ZnWeIwxrjGkkYU+LDRGGpxkLGe922ynzHGVZY0osBpnXowTc/ias9HNPPutMl+xhjXWNKIAumpiZx6yV+JjSnm7Y5f2WQ/Y4xrLGlEiU4dO+PpegWNvpsMuza5HY4xppqypBFN+t7hLJ3+2eNuR2KMqaYsaUST+i2g21WwcCLs/MHtaIwx1ZAljWhz1u3ObPF5j7kdiTGmGrKkEW3qNYP0q2HRJJYsXWwT/owxYWVJIxr1uY1iieG7t+7l8ZmrbMKfMSZsLGlEo7pN+bbJxYyQeTTnR5vwZ4wJG0saUUr63EoBXm7xvkusN8Ym/BljwsKSRpQ6o11bdne8kgs9X/L2JU1swp8xJiwsaUSxxoPHEeONpeP3490OxRhTTVjSiGZ1GkP6NbB4MuxY53Y0xphqwJJGtOt9C8R4bZa4MSYsLGlEu7pNnXkbiydD/ga3ozHGVHGWNKqC3reAxFhtwxgTcq4mDREZLCKrRGSNiNxdRplLRGS5iCwTkdfDHWNUqNcMul2FfjOJiR98ZhP9jDEh41rSEBEP8AwwBGgPXCYi7UuVaQP8Huitqh2AW8MeaJT4Nu1aCooV75dP2gxxY0zIuFnT6AGsUdW1qnoYeAMYUarMdcAzqpoPoKpbwxxj1Phsaw3+W9yfi2Pm0LBwm80QN8aEhJtJoxmwscR5ru9aSacBp4nIFyKSJSKD/d1IRMaKSLaIZG/bti1E4Ua2zFZJvMiFCHB97HSbIW6MCQk3k4b4uaalzr1AG6A/cBkwXkTqH/dNqs+raoaqZiQnJwc90GiQnprIo2MuYHXTYVwe+ynpDQ65HZIxpgpyM2nkAs1LnKcAm/2UmaqqBaq6DliFk0SMH+mpibS/5H5iiovgi3+6HY4xpgpyM2ksANqISEsRiQMuBaaVKvMeMABARBriNFetDWuU0aZBK+h8CWRPgL3Vs6nOGBM6ASUNEUkUkQ4i0kpEgpJoVLUQ+B3wEbACeEtVl4nIgyIy3FfsIyBPRJYDnwJ3qqr18FbkrNuh8CDM/5fbkRhjqhhRLd2N4HtCpB5wI05fQhywDYgHGgNZwLOq+mmY4gxYRkaGZmdnux2G+96+FlZ9CP+3FGo2cDsaY0yEE5EcVc2oqFx5tYa3cUY3naWqbVW1j6+zuTnwMDBCREYHKV4TbH3vhIJ9bPnoCdsS1hgTNN6ynlDVgeU8lwPkhCQiExyN2pGfOoTai8bz/OHOPO2tw6QxmbbvhjGmUgLt0+gsIsNF5JdHjlAHZipvZsMrqSMH+HXMR7YlrDEmKMqsaRwhIhOAzsAyoNh3WYEpIYzLBEHrzmfySXY613o/YFLxUJvwZ4yptAqTBpCpqu0rLmYiTXpqIiuG30f994cztcdKUlIvcjskY0yUC6R5an7phQRN9GiX3g9an0vKihfh8D63wzHGRLlAksYrOIljlYh8KyJLROTbUAdmgqjfONif50z4M8aYSgikeWoCcCWwhJ/7NEw0ad4DWvZzlhbpPgZiE9yOyBgTpQKpafygqtNUdZ2qbjhyhDwyE1z97oJ9W2HhRLcjMcZEsUBqGit9O+a9DxxdOlVVbfRUNEnrA6m94fN/OHuKe2u4HZExJgoFUtNIwEkWg4ALfMewUAZlQqTvnbBnM3Pe/IfNEDfGnJQKaxqqek04AjGhl+M5gxhtQ5vvnmfwyva8PKaPzRA3xpyQMmsaInKPiJS50p2InC0iVuOIIlnrdvBUwS9pJts5v3iuzRA3xpyw8moaS4D3ReQgsJCfV7ltA3QBPgb+FvIITdBktkriaU8XFhe34kbvVLam3eV2SMaYKFNmTUNVp6pqb+B6nCVEPMBu4DWgh6r+n6raLj9RJD01kUljerGx8000l62k75zldkjGmCgTSJ/GamB1GGIxYZCemggtroG8l+Gzx6DzSPAEMojOGGPc3e7VuEXEmSW+Yy0sfdvtaIwxUcSSRnXVdig07gTzHoXiIrejMcZECUsa1ZUI9LsT8tbAUpunaYwJTCD7aSQD1wFpJcur6rWhC8uExekXcKB+Ww58+FfW1RlAesuGbkdkjIlwgdQ0pgL1cIbY/q/EYaJczsZd/H7HYBrsX8ekl56yWeLGmAoFMmympqqOC3kkJuyy1ubxfkF3bohtxm95h1nfX2EzxI0x5QqkpjFdRIaGPBITdpmtkoj1enm66CLaxGxisMx3OyRjTIQLJGncgpM4DorIHt+xO9SBmdBzJvtlcvo5V3Kg/mm0WvaMjaQyxpSrwqShqnVUNUZV432P66hq3WC8uIgM9u0IuEZE7i6n3MUioiKSEYzXNT9LT03kxrNPI2HgH2D7dzaSyhhTroCG3IrIcBF5zHcEZZFCEfEAzwBDgPbAZf72IheROsDNwFfBeF1ThnYjoFEHmPt3q20YY8pUYdIQkYdxmqiW+45bfNcqqwewRlXXquph4A1ghJ9yfwYeAQ4G4TVNWWJioP84yFsNS99xOxpjTIQKpKYxFBioqhNUdQIw2HetspoBG0uc5/quHSUiXYHmqjq9vBuJyFgRyRaR7G3bbA3Fk3b6BdC4o1PbKCp0OxpjTAQKdEZ4/RKP6wXptcXPNT36pEgM8CRwe0U3UtXnVTVDVTOSk5ODFF41FBPjrEmVt8bWpDLG+BVI0ngI+EZEXhaRV4AcgrOPRi7QvMR5CrC5xHkdoCMwR0TWA5nANOsMD7HTh7G/QTt2fvgXctZtdTsaY0yECWT01GScN+wpvqOXqr4RhNdeALQRkZYiEgdcCkwr8bq7VLWhqqapahqQBQxX1ewgvLYpQ87GXdy+bRj1D2xkykuP2ixxY8wxytvu9XTf125AU5yawUbgFN+1SlHVQuB3wEfACuAtVV0mIg+KyPDK3t+cnKy1eXxU2IVFxadyg0zh6zVb3A7JGBNByltG5DZgLPC4n+cUOLuyL66qM4AZpa7dW0bZ/pV9PVOxzFZJxHk9PFF4CRPjHuL8glk4I6KNMQZEVcsvIBKvqgcruhYpMjIyNDvbWrAqI2dDPlnfb+eqVTdQZ98GuHkRxNV0OyxjTAiJSI6qVthnHEhH+JcBXjNVhDNLvA11htwPe3+C7BfdDskYEyHKbJ4SkSY48yYSfPMljgyRrQvYx87qIK03tBoAnz8J6VdDjTpuR2SMcVl5fRrnAVfjDIV9osT1PcAfQhiTiSRn3wPjz4Gv/g1973Q7GmOMy8pMGqr6CvCKiFykqrauRHWVkgGnDYEvnobuYyDB9tswpjorb8jtFb6HaSJyW+kjTPGZCLC83c3ood38OCMYS44ZY6JZeR3htXxfa+PMzi59mGogZ0M+v5yyi/eKelP/2xf5dvkKt0MyxriovOap//i+PhC+cEykyVqbx+HCYp7gIs6Pm4/ns0eg/Utuh2WMcUkgS6M/IiJ1RSRWRD4Rke0lmq5MFedM9othM415S8+l/Y9TYftqt8MyxrgkkHkag1R1NzAMZymR0wAbRlNNHNkS9rZBbel0+V8QbzzM/ovbYRljXFLekNsjYn1fhwKTVXWHiL9VzU1VlZ6aSHqqb9RUrxth3iOwaSE0q/QSZMaYKBNITeN9EVkJZACfiEgytote9XXmTZDQAD6xri5jqqNAlka/G+gFZKhqAbAP/9uymuogvi70vQPWzoHvZ7sdjTEmzALpCI8FrgTeFJG3gdFAXqgDMxEsYzTUawEz74XiIrejMcaEUSDNU88B6cCzvqOb75qprmLjWdvldvhpCetmT3A7GmNMGAWSNLqr6q9VdbbvuAboHurATOTK2ZDP+bOTWVzcioTP/sbC7zdX/E3GmCohkKRRJCKnHjkRkVaAtUlUY1lr8zhUCH8tGEUT2cHhz552OyRjTJgEMuT2TuBTEVmLszx6KnBNSKMyEe3IhL+cwnbM0u6cnfsK7P0/qN3I7dCMMSFWYdJQ1U9EpA3QFidprFTVQyGPzESsIxP+stbm0TTp73jeGwRzHoJhT7odmjEmxCpMGiISD9wA9MHZG/wzEfl3pG73asLjmAl/uaNhwQvQ4zfQ6HR3AzPGhFQgfRoTgQ7A08C/gPbAq6EMykSZfuMgrjbM+pPbkRhjQiyQPo22qnpGifNPRWRxqAIyUahWkrOr36w/wXcz4bRBbkdkjAmRQGoa34hI5pETEekJfBG6kEw0Wth0JPkJqRycPg4KD7sdjjEmRAJJGj2BL0VkvYisB+YD/URkiYh8G9LoTFTI2ZDP5S8t5LbdI4nfvZbcD61D3JiqKpDmqcGhenERGQw8BXiA8ar6cKnnbwPGAIXANuBaVd0QqnjMyTmyUdOn2oXZMV3p881T0O9qqNPY7dCMMUEWyIKFG8o7TvaFRcQDPAMMwelcv0xE2pcq9g3OQomdgbeBR0729UzoHJm34RH4O1fh1cO2Cq4xVVQgNY1Q6QGsUdW1ACLyBs7qucuPFFDVT0uUzwJsx8AIVHLeRmarM4n5bj188ZSzsGFKutvhGWOCKJA+jVBpBmwscZ7ru1aW0cAH/p4QkbEiki0i2du2bQtiiCZQ6amJ3DigtTN3o++dULsxfHAnFBe7HZoxJojcTBr+tv9TvwWdPckzgEf9Pa+qz6tqhqpmJCcnBzFEc1Jq1IFzH4BNObDoNbejMcYEkZtJIxdoXuI8BThuuVQRORf4IzDcli+JIp1HQoszYda9sG+729EYY4LEzaSxAGgjIi1FJA64FJhWsoCIdAX+g5MwtroQozlZMTEw7Ak4tMdJHMaYKsG1pKGqhcDvgI+AFcBbqrpMRB4UkeG+Yo8CtYH/isgiEZlWxu1MhMnZkM8zy2LZ0uE6WDQJ1n/udkjGmCAQVb/dCFErIyNDs7Oz3Q6jWsvZkM+o8VkcLiymrreAr+r9kRrxteD6z8Eb53Z4xhg/RCRHVTMqKudm85Spoo5M9itW2FMYy6y0O2D7KphvmzUZE+0saZigKznZL9YbQ9PuF0K7C2DuI7BjndvhGWMqwc3JfaaKOnayX5Izd6P+3+H7T+F/t8EVU0D8jbg2xkQ6q2mYkDhmsh9AvWZwzn3w/WynY9wYE5UsaZjw6T7Gmbvx4R9g9xa3ozHGnARLGiZscjbu4rXGd1JceBCm/x9UsZF7xlQH1qdhwqLkMNxNsb9i3HevwZK3ofOv3A7NGHMCrKZhwqLkMNzxBYP5sU5HZ0HDvTbR35hoYknDhEXJYbger5e8c5+Ew/vgf7dbM5UxUcSap0xYlB6G2yE1EXbfDZ88CN++BWeMdDtEY0wALGmYsElPTfx5CC7AmbfAdzNhxh3QIhMSU90LzhgTEGueMu7xeOGX/3Gap979DRQXuR2RMaYCljSMa3I25PPMokLW9XwAfpgPnz/hdkjGmApY85RxRckhuE97m/Jl62E0mPMwnHo2NLN9xY2JVFbTMK4oOQS3oFCZ0vR2qN0E3rkODu11OzxjTBksaRhXlF4Jt2vblvCLf8OOtTDjThuGa0yEsuYp4wq/K+FyFvQbB3MfhtRe0O0qt8M0xpRiScO45rghuAD97oKNX8H/7oCmXaBpZ3eCM8b4Zc1TJrLEeOCi8VAzCd66Cg7sdDsiY0wJljRMRMnZkM8zX+9k5VlPw66NMPVG698w1V7Ohnye+XQNORvyy7zmr0woWPOUiRjHDsONYVav39P86z/Dl09D75vdDs+YoMjZkH9MX17p89JlgKN/F3HeGCaNyTzu2r3DOvDg9GXHlDmu6TdILGmYiHHsMNxipsWP4Mb2S+Dj+6BRe2hzrtshGlOuE00A/t7sS5e5qFvKMX8XWWvzAI659sHSLceVsaRhqrwjw3ALCouJ9caQeWpD6P0M7Pge3r4GRs+CRqe7HaYxR51ojaB0AvD3Zg/HJgSFY/8u0urjKdjLFO92ahTtp57nENc23k3i+lXsLoojy9P1aDyhYEnDRAz/w3CBy96A5wfA5JEwZjbUCt0fhDFl8desdKI1gtIJYEjHpnyzfiv1C3fQ1Lub87z7iD2wjaLYb6ivu0iSvZy1LYZ7kvMp2pdHzeJ9eF7ZA8AnHsDjCy4bBnjgp/odyL34tyGrZYDLSUNEBgNP4fzo41X14VLP1wAmAulAHjBSVdeHO04TPn6H4dZLgcsmw0tD4c0r4Kqp4I1zJ0BTbVRUiyjdnHpcjcD3fYneQzQq+ok0z3Z+W2Mlt3TKZe/W9TQhj1qf/cRlnp8Qj2+wxyfOl5tj4KC3DiQkER+TDA2bQ/POEF8fEupDfD2oURdq1IG4Wke/No6vT+N6oUsY4GLSEBEP8AwwEMgFFojINFVdXqLYaCBfVVuLyKXA3wHbeKEa+fkP91TSL3wW3hnt7C8+4l8g4nZ4pgopL0n4q0WUbE6N9yqXn3qYa5N2sW3dEtp4t9Lwk6cg73tyPFt/rhF8DcTWpHG95lCvGTTvhNRNgbpNnWV0ajeC2o2hVjLxEfrByM2aRg9gjaquBRCRN4ARQMmkMQK43/f4beBfIiKqNgazOihd/Z805hzS+94F8x6B+i2g/zi3QzRR6kRZyeryAAAQsklEQVSbmn6uRRTR3JvPebGLaL1hPV+2/gbduoLEAxuJmXIYgNbgvPE3OBVOGwQNWkFiGtRPc/aMqZkU1R943EwazYCNJc5zgZ5llVHVQhHZBSQB20sWEpGxwFiAFi1ahCpeE2alq/9Za/NI7/97Z/7GnL9BzQbQ4zq3wzRR4OSamoRTCrdwRux6bi6ay90pS4jbtoT4wj3wsXPfBvVbQLMOkHwBJLd1jqQ2EF/XvR82xNxMGv5SbekaRCBlUNXngecBMjIyrBZSRRw3mqpVEsTEwPCnnZniM+6EhETodLHboZoIFkiHdWarJJK9+2lftJp07xqu3L6dBxMW4T28y7nJijho3AE6/xKadHKORu2rdHIoi5tJIxdoXuI8BdhcRplcEfEC9YAd4QnPuK3M0VSeWPjVS/DaRc6Of/H1oM1Ad4M1EaN001NZtYhGhVvpGfsdl/w4neTlC/nKsxI8oBKDFLSDjhc6e7uc0gWS29ngCx83k8YCoI2ItAQ2AZcCl5cqMw34NTAfuBiYbf0Z1Yvf0VQAsQnOiKqXh8GbV8JV7zn7jJtqp6KmpyM11kaFP9EndgW3793KfXW/osa+Tc4N1teF5j2dGmvznsgpXZ3RSMYv15KGr4/id8BHOGMLJqjqMhF5EMhW1WnAi8CrIrIGp4ZxqVvxmshxzCfJK6bAS4OdWseo/0LqmW6HZ8KovKan2oU72Z41mfMSVrC43mxq7M11vmlLEqT1gdRbnSX4G7V3Fso0AZGq9sE9IyNDs7Oz3Q7DhMjxI6oySU88AK8Mh92bnImArfq5HaYJkdJNT898uobHZ66iWKGGFHJ7+z3o6pn0ZjHtZQMxolCjHrQ8C1r2c5JF8ulO35g5hojkqGpGReVsRriJKn5HVA1oDdfMcBLH65fApZOgta1TVRVU1PTUp0kRm2Pn0EcX0SdmCXW+P4B6PGyp25ktrX5Js25D4ZSu4LG3umCx36SJKn5HVIEzKerq/8GrI2DyZXDJRGg7xN1gTaX4b3oqoh3rGagLSXn7QRrvWc4ZMbA3rhEH0kZQp8v5SKt+nBJfz+3wqyxrnjJRx99Kokft3wGv/RJ+XALD/gHdrnQnSHNSSv7bZq3N4/GZq/BoIb1jlnF901WkbZ9HE9lBsQr7G3Wldqfz4bTBznDYKJ4wFwmsecpUWaVHVB2bRBo4a1O9dRVM+x3kr4MB91gbdhQoWbNI9B7iqfTtPB07hb6yiDpygKJdCexO7cvHCb1I7nYBZ7Rt7XbI1ZIlDRPV/HaMpybCqLedNao+exzy18OIZyE23u1wTQmla4zfrFrL8OLZDPRm0zdmCTUWFVBQM4nV9QcR13E4rXueT2JsAtZb5S5LGiaq+e0YT010JgAOf9pZ9+eTB2DXJrj0dVtWPUIcSfZ1CvP5KTaHNimrGL1lPhJbyCZtyGQ9l15Dr6Zt94G0t+GwEcWSholqZXaMg9PGfdZtzmJx714P/zkLLn4JWpRe4syEw5GaRZ8mRezNeYcJ8j4941bgEWVnfnOk902sSBzA7F2nkHlqQ9qGcE8Ic/KsI9xEvUD2XGbzIvjvr2FXLpxzH5x5k3WchtHiFauYOvk5BpFFD1lJjCjf6yl8UNSDmdKL+0ZfQnpaA7fDrNYC7Qi3pGGqlDL7OAAO7oKpN8KK9+G0IXDhs85KuSbocjbk8+3KVZyjX9Fiy0x0wxcIyuriZnxQ3JPkzJGc1rEHWet2+B8FZ8LORk+ZaqnMPg5wFja85FX46j8w8x74dx9nWO5pg9wNuirZ8yM/fP4GxVlv8GucGsWB+m3I73IzY3NSWFHYjFhvDJM69XRGwVntIupY0jBVSrl9HOA0SWVeD827w3s3wOu/gs4jYfDDVus4Wbs388MXb+FZOZVTdn1DC5SD2ox/Fv+CD4szueCMc7hxQGse6FrO/BoTNax5ylQ5/vo0/PZzFB5yhuR+9riz9/LQR6HDL6yvIwDfLlvK7oXv0GXPXGpvzQFgVXEKM8mkZd8ruGPuwaOJ+5gmQhOxrE/DGJ9y+zkAflzqTATc/A207AsDH3TWKzLH2rYKVrzPvsXvUStvCQArNJXcpgN55Ie2rC5uhkfgtkFtj87otlpF9LA+DWN8yu3nAGjSEUZ/DNkTYO7D8Hx/6PQrOPtPzp7O1VVxEWzKYcvXU6jx/Qc02L8egD11OvKvwkuZUdSDXJowsnELNm7KxaM/NwmWuQ+KiXqWNEyVV2E/BziroPYcC2dcCl88BfOfgeVTIf0a6HWDM9ejOji0F9bNhVUz4LuPYN82kjWGr7UdH3MNF146loJaTXlpfBYFOL/Pi7qlcFG3FKtZVBPWPGWqhYDmcpS0ezN8+jdY/AZoEbQbDmfeDCnp4Q8+lFRh+2pYPRPWzIINX0LRYQ55arGvxdlkx2dy16Jkdmrto01PNw5oXfHvz0Qd69MwpgwV9nGUtHuzM0Q3+yU4tAuaZ0LXK6D9cGcIbzTa8xOsmwdr5zjHbt+Odsmn82Pjvoxb3IiswtMQbxz3DuvAg9OXWad2NWB9GsaUocI+jpLqngIDH4C+d8DCV2HBC06n+Yw7nP06Oo+EU88Gb43w/hAnYudG+GG+U4v4YT5sW+lcT0gkv3EvFje+kqQu59OpQyfe+XQNnxU6O+F5CovJ33+YSWMyrVZhjrKkYaqdsvo4ym1yqVHH6dvI/C1syoFv34Sl78CydyG2prON6KlnO0fD09wbtntgJ2xZBJsWOqPBNi38uSZRoy407+n027TqT86h5oyasMCpca3MZdKYFL+/G+vUNiVZ85Splvz1cQTcZHVEUQF8/6nTF/D9bMhb41yv3QRO6QJNOkPTM6BpZ6ibErw9PYqLYd9W2PkDbP/OqTlsWwVbV8KuH34ul9jSGTrcIhNa9ILGHcjZuPu4TY6KFeuvMNY8ZUx5Sn969tdkdeR6mW+enlhnCZIjy5Dkr3eSyIYv4cdvnc5lLXaei4mFes2gXnPnqJ0McXUgrhbUqO3UVlShuNDpeC8uhMP74ED+z8e+bU5T0+5NUHS4RBw1nNpNi57Q6Go4pZuTrGo2+DkBHE6CjbuPSYz3Duvgt8ZlNQtTHksaxnB8k1VizbgTr3kkpkHGNc4BcHg/bF3uJJCdPzhv+LtynSGt+7ZD0aGKA5MYZ7Z6QiLUTHJqDu2HO4mnfgtIau28bozn5wThTSK95vG1J2eP7Z8To/VXmJNhScMYnE/XJd9AT6rmUVpcTUjJcA5/igrg8F5nbkTBfhCP04QlHojxQFxtpx8igGYtf81rpX8GBeuvMJVmScMYn9JvoIHUPCrV/u+JdWoQCSf3pl3ytf0ludK1J5uEZ4LBlaQhIg2AN4E0YD1wiarmlyrTBXgOqAsUAX9V1TfDG6mprgKteZROJHCCtZEyVDQZsXTNwl//ROmf4Ug8lixMZbhV07gb+ERVHxaRu33n40qV2Q9cpaqrReQUIEdEPlLVneEO1lRP5dU8/CWSdxbmMmVhboW1kZNJCA9OX1Zu01NZ/RPW/GSCza2kMQLo73v8CjCHUklDVb8r8XiziGwFkgFLGibsyvrUXjKRCFRYGymdAAJJCB8s3VJh05P1T5hwcStpNFbVLQCqukVEGpVXWER6AHHA9+EIzhh/Sr8pl04kAO8szC23NlI6AQSSEIZ0bMqC9TsCanoyJtRCljRE5GOgiZ+n/niC92kKvAr8WvXIoPfjyowFxgK0aNHiBCM15uSVTiQV1UZKJ4BAE0LbJnWs6clEBFdmhIvIKqC/r5bRFJijqm39lKuL03T1kKr+N5B724xwE2lOtE/DGDdE9Cq3IvIokFeiI7yBqt5Vqkwc8AHwvqr+I9B7W9IwxpgTF2jSCNJiOCfsYWCgiKwGBvrOEZEMERnvK3MJ0Be4WkQW+Y4u7oRrjDEGbMFCY4wxRH5NwxhjTBSypGGMMSZgljSMMcYEzJKGMcaYgFW5jnAR2QZsqMQtGgLbgxROOFi8oWXxhpbFG1onEm+qqiZXVKjKJY3KEpHsQEYQRAqLN7Qs3tCyeEMrFPFa85QxxpiAWdIwxhgTMEsax3ve7QBOkMUbWhZvaFm8oRX0eK1PwxhjTMCspmGMMSZgljR8RGSwiKwSkTW+lXcjmohMEJGtIrLU7VgqIiLNReRTEVkhIstE5Ba3Y6qIiMSLyNcistgX8wNux1QREfGIyDciMt3tWAIhIutFZIlvMdKIXzBOROqLyNsistL3f7mX2zGVRUTalljodZGI7BaRW4Nyb2uecv7YgO9wVtzNBRYAl6nqclcDK4eI9AX2AhNVtaPb8ZTHt2dKU1VdKCJ1gBzgwgj//QpQS1X3ikgs8Dlwi6pmuRxamUTkNiADqKuqw9yOpyIish7IUNWomPcgIq8An6nqeN/WDTVVNeK3n/a9v20CeqpqZeawAVbTOKIHsEZV16rqYeANnH3MI5aqzgN2uB1HIFR1i6ou9D3eA6wAmrkbVfnUsdd3Gus7IvYTloikAOcD4ysqa06cb0O4vsCLAKp6OBoShs85wPfBSBhgSeOIZsDGEue5RPibWrQSkTSgK/CVu5FUzNfcswjYCsxS1UiO+R/AXYDfLZEjlAIzRSTHt2VzJGsFbANe8jUBjheRWm4HFaBLgcnBupklDYf4uRaxnyqjlYjUBt4BblXV3W7HUxFVLVLVLkAK0ENEIrIZUESGAVtVNcftWE5Qb1XtBgwBbvQ1uUYqL9ANeE5VuwL7gGjo+4wDhgMBbZcdCEsajlygeYnzFGCzS7FUSb5+gXeASao6xe14ToSvGWIOMNjlUMrSGxju6yN4AzhbRF5zN6SKqepm39etwLs4zcSRKhfILVHbfBsniUS6IcBCVf0pWDe0pOFYALQRkZa+zHwpMM3lmKoMX6fyi8AKVX3C7XgCISLJIlLf9zgBOBdY6W5U/qnq71U1RVXTcP7vzlbVK1wOq1wiUss3KAJfM88gIGJHAqrqj8BGEWnru3QOELEDOUq4jCA2TYFT5ar2VLVQRH4HfAR4gAmquszlsMolIpOB/kBDEckF7lPVF92Nqky9gSuBJb4+AoA/qOoMF2OqSFPgFd/IkxjgLVWNiqGsUaIx8K7zeQIv8LqqfuhuSBW6CZjk+2C5FrjG5XjKJSI1cUaE/iao97Uht8YYYwJlzVPGGGMCZknDGGNMwCxpGGOMCZglDWOMMQGzpGGMMSZgljSMMcYEzJKGMUHgWzb7hnKeTxCRub55H2WV+VhEEkMToTHBYUnDmOCoD5SZNIBrgSmqWlROmVcruIcxrrOkYUxwPAyc6tvw5lE/z48CpoKzv4iIzPOVXSoiZ/nKTMNZ9sGYiGUzwo0JAt+S79P9bYjlW3biB1Vt4ju/HYhX1b/6mqtq+vYZQURWA5mqmhe24I05Abb2lDGh1xAouWHPAmCCb+Xf91R1UYnntgKnAJY0TESy5iljQu8AEH/kxLfrYl+cLThfFZGrSpSN95U3JiJZ0jAmOPYAdfw9oar5gEdE4gFEJBVn06QXcJaM7+a7LkATYH04AjbmZFjSMCYIfH0QX/g6tv11hM8E+vge9wcWicg3wEXAU77r6UCWqhaGOl5jTpZ1hBsTBiLSFbhNVa8sp8xTwDRV/SR8kRlzYqymYUwYqOo3wKflTe4DllrCMJHOahrGGGMCZjUNY4wxAbOkYYwxJmCWNIwxxgTMkoYxxpiAWdIwxhgTsP8H2syjTiub7bQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs_sol = analytical_solution_matrix(A_num, x0, h, num_steps)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t_range,[x[0,0] for x in xs_e], 'o', markersize=3, label='x_euler')\n",
    "plt.plot(t_range,[x[0,0] for x in xs_sol], '-', markersize=3, label='x_solution')\n",
    "plt.xlabel('t (s)')\n",
    "plt.ylabel('position (m)')\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
}