{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dCdt(C, t):\n",
    "    \n",
    "    CA = C[0]\n",
    "    CB = C[1]\n",
    "    CC = C[2]\n",
    "    CD = C[3]\n",
    "    k1 = 0.9\n",
    "    k2 = 0.6\n",
    "    \n",
    "    dCAdt = -k1*CA*CB - 2*k2*CA**2\n",
    "    dCBdt = -k1*CA*CB\n",
    "    dCCdt = k1*CA*CB\n",
    "    dCDdt = k2*CA**2\n",
    "    \n",
    "    return [dCAdt, dCBdt, dCCdt, dCDdt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "CA = np.ones(1001)\n",
    "CB = 2*np.ones(1001)\n",
    "C = [CA, CB]\n",
    "\n",
    "t = np.linspace(0, 10, 1001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-31-219d3be89acf>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdCdt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-29-02ccf0753ff3>\u001b[0m in \u001b[0;36mdCdt\u001b[1;34m(C, t)\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[0mCA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mCB\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m     \u001b[0mCC\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m     \u001b[0mCD\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m     \u001b[0mk1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.9\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "dCdt(C, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = odeint(dCdt, [2, 1, 0, 0], t)\n",
    "CA = C[:, 0]\n",
    "CB = C[:, 1]\n",
    "CC = C[:, 2]\n",
    "CD = C[:, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1001, 4)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x18da6f44b70>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXwdVf34/9f77tn3pEuSNt3oAt2o0AoIgkARERGRggIfQbEs4vbVz0c/n48oLj/UjygfwYVNcKGIKFgUqXxkUyjQQgvdS/ekSZs0afbc5C7n98dMkpv0JrltbnKTe9/Px2Mec+acMzPnsrzn5MzMGTHGoJRSKnk5Et0ApZRSI0sDvVJKJTkN9EopleQ00CulVJLTQK+UUknOlegGRFNYWGimTp2a6GYopdS48eabbx4xxhRFKxuTgX7q1KmsX78+0c1QSqlxQ0T2D1SmQzdKKZXkNNArpVSS00CvlFJJbkyO0Sul1GgLBAJUVVXh9/sT3ZRB+Xw+SktLcbvdMe+jgV4ppYCqqiqysrKYOnUqIpLo5kRljKG+vp6qqioqKipi3m/IoRsRKRORF0Rkm4hsEZHPR6kjIvK/IrJLRN4RkcURZdeJyLv2cl3MLVNKqVHk9/spKCgYs0EeQEQoKCg47r86YunRB4EvG2PeEpEs4E0Rec4YszWizkXATHs5Hfg5cLqI5AO3A0sAY++72hhz9LhaqZRSo2AsB/luJ9LGIXv0xpgaY8xbdroF2AZM7lftUuDXxvIakCsiE4ELgeeMMQ12cH8OWH7crYzR//7jXV7aWTdSh1dKqXHpuJ66EZGpwCLg9X5Fk4HKiO0qO2+g/GjHvlFE1ovI+rq6EwvWv3xpNy/t0ECvlBq/nnzySUSE7du3x+2YMQd6EckE/gh8wRjT3L84yi5mkPxjM425zxizxBizpKgo6lu8Q8r0uWjrDJ7QvkopNRasWrWKM888k8ceeyxux4wp0IuIGyvI/84Y86coVaqAsojtUqB6kPwRkeF10aqBXik1TrW2tvLKK6/w4IMPxjXQD3kzVqyR/weBbcaYuwaothq4VUQew7oZ22SMqRGRNcD3RCTPrncB8LU4tDuqTA30Sqk4+NbTW9ha3X/gYnjmTsrm9kvmDVrnqaeeYvny5cyaNYv8/HzeeustFi9ePOg+sYilR38GcA1wrohstJcPishKEVlp13kG2APsAu4HbgYwxjQA3wbW2csddt6IyPTq0I1SavxatWoVK1asAGDFihWsWrUqLscdskdvjPkX0cfaI+sY4JYByh4CHjqh1h2nDK+Lhrb20TiVUiqJDdXzHgn19fU8//zzbN68GREhFAohIvzgBz8Y9mOfSTXXTZYO3SilxqknnniCa6+9lv3797Nv3z4qKyupqKjgX//617CPnVSBXm/GKqXGq1WrVnHZZZf1ybv88st59NFHh33spJrrRh+vVEqNVy+++OIxebfddltcjp1UPfpMr4tAyNAZDCW6KUopNWYkXaAHaPVrr14ppbolVaDPsAN9W6f26JVSqltSBfpMrxOAls5AgluilFJjR5IFeuuLK9qjV0qpXkkV6DPsHr0+eaOUUr2SKtBn+awx+hYN9EqpccjpdLJw4UIWLFjA4sWLefXVV+Ny3KR6jr73ZqwGeqXU+JOWlsbGjRsBWLNmDV/72td46aWXhn3cpOrR6+OVSqlk0dzcTF5e3tAVY5BcPXqPHei1R6+UGo6//Qcc2hTfY044BS66c9AqHR0dLFy4EL/fT01NDc8//3xcTp1Ugd7hENI9Th26UUqNS5FDN2vXruXaa6/tmc1yOJIq0IN+fEQpFQdD9LxHw7Jlyzhy5Ah1dXUUFxcP61hJNUYPGuiVUslh+/bthEIhCgoKhn2sWD4l+BDwIaDWGHNylPKvAJ+ION4coMgY0yAi+4AWIAQEjTFLht3iIWT6NNArpcan7jF6AGMMjzzyCE6nc9jHjWXo5mHgHuDX0QqNMT8EfgggIpcAX+z3ucD3G2OODLOdMcv0uvSpG6XUuBQKjcxb/UMO3RhjXgZi/c7rVUB8PnJ4grJ9bpr9OteNUkp1i9sYvYikA8uBP0ZkG+DvIvKmiNw4xP43ish6EVlfV1d3wu3ISXPT1KGBXimlusXzZuwlwCv9hm3OMMYsBi4CbhGR9w20szHmPmPMEmPMkqKiohNuRHaaSwO9UkpFiGegX0G/YRtjTLW9rgWeBE6L4/miyklz4w+E9StTSilli0ugF5Ec4GzgzxF5GSKS1Z0GLgA2x+N8g8lJs6Yqbu7QG7JKKQWxPV65CjgHKBSRKuB2wA1gjPmFXe0y4O/GmLaIXUuAJ+03ulzAo8aYZ+PX9Oiy7UDf1BGgKMs70qdTSqkxb8hAb4y5KoY6D2M9hhmZtwdYcKINO1GRgV4ppcaTQ4cO8YUvfIF169bh9XqZOnUqP/nJT5g1a9awjpt0b8b2Dt1ooFdKjR/GGC677DLOOeccdu/ezdatW/ne977H4cOHh33spJvrpifQ67P0Sqlx5IUXXsDtdrNy5cqevO63ZIcraQO9Dt0opU7U99/4Ptsbtsf1mLPzZ/Pvp/37gOWbN2/m1FNPjes5uyXd0E22zw707RrolVIKkrBH73E5SHM7tUevlDphg/W8R8q8efN44oknRuTYSdejB50GQSk1/px77rl0dnZy//339+StW7dOvxk7kJw0ndhMKTW+iAhPPvkkzz33HNOnT2fevHl885vfZNKkScM+dtIN3YD26JVS49OkSZN4/PHH437cpOzRWxOb6RQISikFSRvo3frClFJK2ZIy0OvQjVJK9UrKQJ+b5qG1M0ggFE50U5RSKuGSMtDnZ3oAONreleCWKKVU4iVnoE+3An1DmwZ6pZRKzkCfoYFeKTX+OJ1OFi5cyLx581iwYAF33XUX4fDwh6CT8jl6DfRKqfEoLS2NjRs3AlBbW8vVV19NU1MT3/rWt4Z13CF79CLykIjUikjUzwCKyDki0iQiG+3lGxFly0Vkh4jsEpH/GFZLj0N3oD+qgV4pNU4VFxdz3333cc8992CMGdaxYunRPwzcA/x6kDr/NMZ8KDJDRJzAvcD5QBWwTkRWG2O2nmBbY5abbs1gWa+BXil1Ag5973t0bovvNMXeObOZ8PWvH9c+06ZNIxwOU1tbS0lJyQmfe8gevTHmZaDhBI59GrDLGLPHGNMFPAZcegLHOW5up4OcNLf26JVS495we/MQvzH6ZSLyNlAN/D9jzBZgMlAZUacKOH2gA4jIjcCNAOXl5cNuUH6GR3v0SqkTcrw975GyZ88enE4nxcXFwzpOPJ66eQuYYoxZAPwUeMrOlyh1B7w0GWPuM8YsMcYsKSoqGnaj8jM8+hy9UmrcqqurY+XKldx6662IRAunsRt2j94Y0xyRfkZEfiYihVg9+LKIqqVYPf5RkZfuoepo+2idTimlhq2jo4OFCxcSCARwuVxcc801fOlLXxr2cYcd6EVkAnDYGGNE5DSsvxLqgUZgpohUAAeBFcDVwz1frAoyPGw62Dhap1NKqWELhUIjctwhA72IrALOAQpFpAq4HXADGGN+AXwMuElEgkAHsMJYdw+CInIrsAZwAg/ZY/ejIi/DQ0NbF8aYYf/Zo5RS49mQgd4Yc9UQ5fdgPX4ZrewZ4JkTa9rwFGR4CIQMrZ1BsuwPhiulVCpKyikQwOrRAxxt0+mKlVKxicejjCPtRNqYtIG+wJ7Bsq61M8EtUUqNBz6fj/r6+jEd7I0x1NfX4/P5jmu/pJzrBqAo0wtAXYsGeqXU0EpLS6mqqqKuri7RTRmUz+ejtLT0uPZJ2kBfnN0d6P0JbolSajxwu91UVFQkuhkjInmHbjK8OER79EoplbSB3ukQCjK91GqgV0qluKQN9GCN02uPXimV6pI60Bdna49eKaWSOtAXZXqp1ZuxSqkUl9SBvjjby5HWLsLhsftcrFJKjbTkDvRZPkJhQ4NOV6yUSmFJHeiLsqxn6WubdZxeKZW6kjrQF9uBXqdBUEqlsqQO9L09er0hq5RKXUkd6EuyrYl/DmugV0qlsCEDvYg8JCK1IrJ5gPJPiMg79vKqiCyIKNsnIptEZKOIrI9nw2PhczvJz/BwsFEDvVIqdcXSo38YWD5I+V7gbGPMfODbwH39yt9vjFlojFlyYk0cnok5PmqaOhJxaqWUGhOGDPTGmJeBhkHKXzXGHLU3X8P6CPiYMTEnjRrt0SulUli8x+hvAP4WsW2Av4vImyJy42A7isiNIrJeRNbHcz7oSbk+qrVHr5RKYXGbj15E3o8V6M+MyD7DGFMtIsXAcyKy3f4L4RjGmPuwh32WLFkSt1dZJ+ak0eIP0toZJNObtNPvK6XUgOLSoxeR+cADwKXGmPrufGNMtb2uBZ4ETovH+Y7HpFzryZuaRu3VK6VS07ADvYiUA38CrjHG7IzIzxCRrO40cAEQ9cmdkTQpNw2A6iYdp1dKpaYhxzJEZBVwDlAoIlXA7YAbwBjzC+AbQAHwMxEBCNpP2JQAT9p5LuBRY8yzI/AbBjUxR3v0SqnUNmSgN8ZcNUT5p4FPR8nfAyw4do/RVZLtQ0R79Eqp1JXUb8YCuJ0OirO8VGuPXimVopI+0AOU5aVT2dCe6GYopVRCpESgLy9I54AGeqVUikqNQJ+fzqFmP/5AKNFNUUqpUZcSgX5KQTrGQNVRHadXSqWelAj05fkZABxoaEtwS5RSavSlSKBPB+BAvY7TK6VST0oE+sJMD+keJ/v1hqxSKgWlRKAXEcrz9RFLpVRqSolAD9bwzX4dulFKpaCUCvQHGtoJh+M2A7JSSo0LKRPopxSk0xkMU9vSmeimKKXUqEqZQF9RmAnAnrrWBLdEKaVGV8oE+hnFVqDfpYFeKZViUibQl2R7yfS62FWrgV4plVpSJtCLCDOKM3n3sAZ6pVRqiSnQi8hDIlIrIlE/BSiW/xWRXSLyjogsjii7TkTetZfr4tXwEzGjOFOHbpRSKSfWHv3DwPJByi8CZtrLjcDPAUQkH+vTg6djfRj8dhHJO9HGDteM4kzqWjppag8kqglKKTXqYgr0xpiXgYZBqlwK/NpYXgNyRWQicCHwnDGmwRhzFHiOwS8YJy4UgD99Ft7+/YBVZhR135BtGZEmKKXUWBSvMfrJQGXEdpWdN1D+MUTkRhFZLyLr6+rqjr8FTjfsfRl2Pz9glZkldqDXG7JKqRQSr0AvUfLMIPnHZhpznzFmiTFmSVFR0Ym1YsIpcOidAYtL89LxuBwa6JVSKSVegb4KKIvYLgWqB8kfGRNOgbodEPBHLXY6hOlFmezQJ2+UUikkXoF+NXCt/fTNUqDJGFMDrAEuEJE8+ybsBXbeyJhwCpgQ1G0bsMrcidlsrW4esSYopdRYE+vjlauAtcBJIlIlIjeIyEoRWWlXeQbYA+wC7gduBjDGNADfBtbZyx123siYcIq1PrRpwCrzJmVzpLWT2ubovX6llEo2rlgqGWOuGqLcALcMUPYQ8NDxN+0E5FWAJwtqBh6nnzcpG4At1c0UZ/tGpVlKKZVIyfVmrMMBE06GmrcHrDLHDvRba3T4RimVGpIr0AOULoGajRCMPh1xts9NeX46W6qbRrlhSimVGMkX6MtOh1DXoL36eZP0hqxSKnUkX6AvPc1aV74+YJW5E7PZV99Oi1+nQlBKJb/kC/RZJZA3ddBAf/LkHAA2H9RevVIq+SVfoAerV1/5Bpjo34ddWJYLwIbKo6PZKqWUSojkDPRTz4DWw3BkZ9TivAwP0wozeGt/4yg3TCmlRl9yBvpp77fWg0xwtrA8lw0HjmIG6PUrpVSySM5AnzcF8qfD7hcGrLK4PI/6ti4qGzpGsWFKKTX6kjPQA0x/P+z7FwS7ohYvLre+f/LWAR2nV0oltyQO9OdCoA2q3ohaPKskk3SPUwO9UirpJW+gn3oWOFzw7t+jFrucDhaV5/LG3pGbY00ppcaC5A30vmyoeB9se3rAxyyXVhSw/VALDW3Rh3eUUioZJG+gB5jzYWjYA4c3Ry1+74wCAF7fUz+arVJKqVGV3IF+9odAHLB1ddTi+aW5pHucrNVAr5RKYskd6DOLYMoZsPWpqMM3bqeDJVPzWbtbA71SKnnF9OEREVkO3A04gQeMMXf2K/8xYL+lRDpQbIzJtctCQPcnnw4YYz4cj4bHbN5l8NcvWbNZTlp4TPGyaQV8/9nt1LV0UpTlHdWmKaWOnzGGkAkRDAd71oFwgFA4RNAErfxwiEA4QNBY6e66IRMiHA73pEMmRCgcImzCg26HTESeXR40wWPzwsE++0aWRz12dx27PNuTzc8+8LO4/zMbMtCLiBO4Fzgf62Pf60RktTFma8Q/+C9G1P8csCjiEB3GmGMj7Gg5+XJY83XY8Nuogf7MGYV8H/jXrjouW1Q6+u1TKkG6A1NXqItAOEAgHOiTDoQCUfP7p7sDbU9+KEBX2MrvDrCBcKAnAEcLxj1l4WDPduT+Pdt22VggCE5x4nQ4cYgDl7hwOBxWnlh53eU923a6Tx2HE5fDhUc8ZLgzRqStsfToTwN2GWP2AIjIY8ClwNYB6l8F3B6f5sVBWi7MuQQ2PQ4XfBvcaX2K503KpjDTyz+21WqgV6MqbML4g346Q510hjrpCnXhD/npCnVZeUErvzPcm+5fZ8B9IpbuANw/iAdN/AOmS1y4nW5cDhduR+/aKVYwczlcOMXZU+Z0OPGKtyftdrhxiasn+LkcLlzi6k3b+0crczqcPdtuh7vPdne6f7CNGnz7BWuXw9UnaHenHTJ+Rr5jCfSTgcqI7Srg9GgVRWQKUAFETjLjE5H1QBC40xjz1AD73gjcCFBeXh5Ds47Domtg0x9g219g/hV9ihwO4dzZRTy7+RCBUBi3c/z8y1MjLxQO0R5spy3QRlugjY5gR9TFH/THVNYnHRreB+rdDjc+pw+P04PX6cXr8uJ1evE4PficPtLd6T3bHocHt9ON2+Hum3Z6eoJudzoyv3+6e79odboDohp7Ygn0EiVvoJnAVgBPGGNCEXnlxphqEZkGPC8im4wxu485oDH3AfcBLFmyJL4zjU09y/pw+Bv3HRPoAc6dXczj66t4c/9Rlk4riOupVWIEwgFaulpo7mymuauZ1kAr7YF2WgOtPUE7coksj6zXETy+uZB8Th9prjR8Lmvdnc7z5THJNaknL7KO1+k9dukXtHuCubM3X4OqilUsgb4KKIvYLgWqB6i7ArglMsMYU22v94jIi1jj98cE+hHlcMDSm+FvX7HmqS87rU/xmTOLcDuFF7bXaqAfQ8ImTEtXC0f9R2nsbKSps4nmrubepbM33dLV0icvlgDtdrjJdGeS7k4n051JhjuDAl8B5VnlZLgzyHBnHFOe7k7vCebRAroGXzUWxRLo1wEzRaQCOIgVzK/uX0lETgLygLUReXlAuzGmU0QKgTOAH8Sj4cdt4dXwwnfg1Z/Clb/pU5TpdXF6RQH/2F7L1z44JyHNSwVhE6bB30B9Rz0N/gYaOxtp8Df0BPLuvKP+ozT4G2jqbCLU54/DvtJd6WR7s8n2ZJPlyaI0s5Ts/Ow+edkeK53pyTwmaHucnlH89UolzpCB3hgTFJFbgTVYj1c+ZIzZIiJ3AOuNMd1vI10FPGb6TvA+B/iliISxntm/M/JpnVHlzYQl18Mrd0PdTiia1af4vDnFfOvpreypa2VaUWZCmjheBcIB6jvqqWuvo66jjiMdR6jrqKOuvTd9pP0I9f76qIFbEHK8OeT58sjz5jElewoLixeS580jz5dHrjeXPF8eOZ6cniCe6cnE7XAn4NcqNf7IWPzwxpIlS8z69evjf+C2I3D3Aph5PlzxcJ+iQ01+lt35D774gVncdt7M+J97nDLG0BJooaa1hkNth6huq6amrYZDrYeoaauhuq2auvY6TL/bNoKQ58ujKK2IwvRCitKKrHRaIYVpheT58sj35fcEcKfDmaBfqFRyEJE3jTFLopXF9MJU0sgohKU3wcs/hDO/CBMX9BRNyPHxnqn5rH67ms+dOwORaPegk1PYhDncdpgDLQc40HKAyubKnnR1azVtgbY+9d0ONxMyJjAxYyJLJy5lYsZEitOLKUoroji9mMK0QvLT8rXHrdQYkVqBHmDZrdbTN89/Bz7xhz5Fl8yfyH//eQs7Drcwe0J2gho4ctoD7exu3M2uxl3sbtzN/ub9HGg5QFVLFV3h3hk83Q43ZVlllGWVcdqE05iYMbEnsE/KnES+L19vOio1jqReoE/LhbO+DM99A3augVkX9hRddMpEbl+9hb+8XTOuA30gHGBP4x62N2zvCeq7G3dT3db7sJTH4aE8u5yKnArOLj2bsuwyyrPKKc8qpzi9WIdSlEoiqRfoAU6/CTb8Dp75ijVnvf22bGGml/dOL2T129V86fxZOBxjf/imK9TFrsZdbK3fyrb6bWyt38rOozt7euhuh5uKnAoWFC/g8tzLmZ4znem50ynNKsXlSM1//UqlmtT8P93lgYv/Bx65BP75Izj3v3qKPrp4Ml96/G1e21vPe6cXJrCR0dV31LOxdiMbajewoW4D2+q3EQgHAMhyZzGnYA5Xz7maOflzmF0wm/Kscg3oSqW41I0AFe+D+Svgn3fBSRfB5FMB+OApE/nm6i38fl3lmAj0te21vFbzGm/UvMHGuo3sb94PWEMvJxeezCfnfpJ5BfOYmz+X0qzSlLqJrJSKTeoGeoCL7oR9/4Q/fRY++zJ40vG5nVy2aDKr1lXyrfYuctNH96Wa9kA7bx5+k7U1a1lbvZZdjbsAyPXmsrB4IZfPvJxFxYuYWzBXX/hRSsUktQN9Wh585Gfw60vh2X+HD/8UgCvfU84ja/fz5IaDfOqMihFvRoO/gZcqX+L5A8/zavWrdIW78Dg8LC5ZzCXTL2HZxGWclH9S0j3pYoyBYBATDGICAUwoBKHQsetwGBMMQjiMCYYgHFknDKGgtbbz++wfuU8wBCaMCYetD9GETcQ2EA7b28Yq7y4bdNvY+xmMscvCYcBE1A1bv7Wnbr9tjLVtev7B9C5WBj3vuxj6ltnpnvcY4nYcjinr3jYMco7uY/Qv75M+9t2dY97n6V/lmPLj3O5udzzPESVvuOdw5ucz7ck/HXueYUrtQA8w7RzrKZx//ggmLYIl1zN3UjYLynL5zWv7uW7Z1BG5KVvfUc+z+57luf3PsaF2A2ETZmLGRK446QreN/l9LC5ZjM/li/t5I5lQiHBbW+/i78R0+gn7/ZjOTozfH5Fnrzv8Vn6nH+PvxHR1WkE6YAfryCV4bB79t8cTEWveJBFriMzhiLot0FvmEISIMpHe40SWdR+/e+hNsIfhJEqZne6pKtHLu4/DQGXSW25viziOKbeS0t2oqOdApM9xop0n6nH6//MddLt/cf9jDLH/SJwjpmMOfo7IcmfGyLyVr4Ee4P3/CTXvwDNfheK5UL6U68+Yyucf28iLO2s5d3ZJXE7jD/p5sfJFnt7zNK8cfIWQCTEjdwafOeUznFt+LnPy5xzXGLsJBgk1NxNqbCTU2GSvI5bmJsKtEYG8tZVwWxshe9t0HN/MjIAVrHw+HF6vtfZ4EI8Hcbt7F58XR1Ym4o7Id7n61nG7EXdvHi4X4nQhLic4nL1rpwOcTsTp7F07HFZdpwO615H7uLrr9NvP6UTswNwToCO2JSL/mO3uYK7UOJRaUyAMpuMo3H8udDTC9WsI5M/grO+/wLSiDB79zNJhHXpv015+v+P3rN61mpZAC8XpxVw87WIumXYJM/OOnW7BhEIEDx0iUF1NsK6OQG0twbo6a6mt60mHm5sHPqnTiTM7G0dWFo6MDBwZ6TgzMu10Bo7MyHQGjvR0HD6fFbx9PsTrw+Hz9g3qXi+43RrwlBqDdAqEWKTlwSf/CA9eCL+5DPcNa/i3M6Zy59+2s6W6iXmTco7rcGET5qXKl3h0+6O8VvMaLoeL86ecz0dnfpT3lLwHh4FAVRUtb/4fXfv20VVZRaCykq6qKgLV1RDs+/Uf8XhwFRXhKi7GO306GUuX4szLw5mbizMnx1rn5uLMs7YdmZkakJVSgPboj1XzDjx8MWQW0/zxP7Ls3u2cM7uYe69eHNPuwXCQNfvW8MCmB9jVuIuS9BKunHoZHwrOwbNtP507ttP57i46d+/GdHb27OfMzcVdVoanrBT35FLcZaW4J03GXVKMq6gIR06OBm6l1IAG69FroI/mwOvwu4+BL5f7p/2Y767t5NkvnDXotAihcIi/7PkLv3znlzTUHuD8uhIubppK8d5GOnfs7Omhu0pK8M6YgXfmTLwzZ+CdMQPPtGk4s7JG69cppZKQBvoTUb0BfvNRwg43H2/9MoUzlvCLa06NWnXtwVd59M/fpWT9Xk6v9DH5oB8JGxzp6fjmzydtwQLSFswnbf58XIWJfwlLKZV8hj1GLyLLgbuxPjzygDHmzn7l/wb8EOsLVAD3GGMesMuuA7rnGPiOMeaR4/4FiTBpEXzqGRy/vZxVzm9wy7ab2HxwBidP7h2r3//2K7z84B2UvnGAmxvBOBykzz+JjA+fQcYZ7yXtlFOsJ0qUUiqBhuzRi4gT2Amcj/X92HXAVZFfirID/RJjzK399s0H1gNLsF4deBM41RhzdLBzjokefbeWQwQfvRpXzZs8mXMtH7nphzT+33PsfOinZG+tJCTQMn8qMz96HbkXXIgrLy/RLVZKpaDh9uhPA3YZY/bYB3sMuBSI5ZOAFwLPGWMa7H2fA5YDq2Jp+JiQNQHX9c+w45ef4uy1T7DtzOeR1jDtubD5IzO46JYfcHKZfmdWKTV2xRLoJwOVEdtVwOlR6l0uIu/D6v1/0RhTOcC+k6OdRERuBG4EKC8vj6FZo8OEwzT95VnkDzUcPpzLu5PDPLvcwUfOv5Lrz/ov64UbpZQaw2IJ9NGe6es/3vM0sMoY0ykiK4FHgHNj3NfKNOY+4D6whm5iaNeI8+/YwaFv3UHHW29RV5HHz69ysKOwjN92HWXuCz+Cg1vg4h9BTmmim6qUUgOKpTtaBZRFbJcC1ZEVjDH1xpjuh8LvB06Ndd+xyITD1D/4IHs/dgX+PbtZ/fFybr2ymfd95FYWZRsgoxwAABhASURBVN7Bxxu+QcOZt8Pel+He0+GVuyHgT3SzlVIqqlgC/TpgpohUiIgHWAGsjqwgIhMjNj8MbLPTa4ALRCRPRPKAC+y8MSvU1ETlpz9D7Q//BzljCV+9KY0/nnSUH597NzctuInvfOQUxOHipt3LCK98FaaeaX2W8J73wKYn7NkIlVJq7Bgy0BtjgsCtWAF6G/C4MWaLiNwhIh+2q90mIltE5G3gNuDf7H0bgG9jXSzWAXd035gdi7oOHGDflStoX7eO8Fc+y6fP3kFbmoPfXPQbzis/D4CJOWn814fm8PreBn63U+Dq38O1f4a0HPjjDXDf2bD1zxrwlVJjhr4wZfNv386BT10PxtD+7dtYWXc3ud5c7r/gfsqyyvrUNcZw7UNvsH7fUZ7+3BnMKM6yAvs7v4eXfwANe6DwJGv645M/Ck59ll4pNbIGe7xSHxkB/Dt3cuBT1yNeL6FffIeVdXdTlFbEw8sfPibIgzVH9f9csYB0j5Obf/cWHV0hazrbhVfBrevh8gfB4YQnb4SfnAIvfh9aDifglymllAZ6AjU1HLj+BsTtxn3v91i54w6yPFncf8H9TMiYMOB+Jdk+fnzlQt6tbeX21Zt7CxxOOOVjsPIVuOoxKJ4DL34PfjwX/vAp2P0ChEOj8MuUUsqS0tMUh9vbqbzlFkxHB/m/eYDrtv8nYRPml+f/ctAg3+19s4q49f0z+Onzu1hUnsdVp0U8/+9wWB8dP+kiOLIL1j8IG34HW/4EWRPhlCtgwQoomTeCv1AppVJ4jN4YQ/VXvkrzX//KpJ/dw7+HHuf1mtd56MKHWFS8KObjBENhrn9kPa/uOsIj15/GGTMGmbQs4Iedf4O3fw+7noNwEIrmwOyLYc6HYOLC6J8/U0qpIejslVE0Pf0Xqr/yFQpv+xyPLg3wwKYH+O+l/83HT/r4cR+r2R/gYz9/lZomP0/e/F7r5uxQ2o7A5j/BttWw/xUwYcgutYL+zAtgynvBk34Cv0wplYo00PcTOHyYPR+6BO/06dT88BZufP4mLp95Od987zdP+JhVR9v5yL2v4HU5eXzlMibnpsW+c1s97HwWtv8Vdv8Dgn5weqB8KUw/F6a9HybM7/2ItFJK9aOBvp+DX/oyLf/4B0VP/JYrN36BNFcaj1/yOGmu4wjOUWw+2MRV979GQYaHxz+7jOJs3/EfpKsdDrxq3bTd/QLUbrHy0/Kh7HQoPx3KllrTKLtP4PhKqaSkgT5C+7p17L/mWgpvuYW7Fh7k2b3P8tsP/pZ5hfG5Kfrm/qNc8+DrTMpN49HPnE5x1jCDccth2POiNd1C5WtQv8vKd3qsYF92mrWeuBDyKrTXr1SK0kBvM8aw7/KPEWw8SsOD3+TTL9/CjfNv5HOLPhfX86zdXc/1D6+jONvLb284nbL8OI61tx2BytfhwFrrk4c1GyHUZZV5s60hnkkLYeICK10wXV/YUioFaKC3tTz/AlU330zxd77Np1y/pivUxVOXPoXPFf8hkLcOHOVTv1qH1+Xg1zecNuj3Zocl2AV126B6I9S8bS2HN1vj/AAONxTOhKLZ1jP9RSdZT/rkTwNnSj9dq1RSGfanBJOBMYYjP/857tJSnp7RyN6393LvefeOSJAHWFyexx9WLuOaB1/nil+s5Z6rF3P2rKL4n8jlsXrvExf05oUCcGQn1LxjXQTqdsDBN61n+Ls5PVawz58O+RV22l5ySq0Xv5RSSSFlevRtr77KgetvIPf2r3M5P2dh0UJ+9oGfxfUc0VQdbefTj6xn5+EWvnbRHD59VgWSqGflu9qsoF+3HWq3Qf1ua16eo3t7/wIA66+AvKlW0M+bYgX+nFLIKbOWzBK9F6DUGKM9eqDh0Udx5ufzxLQjtGxv4bbFt43KeUvz0vnTze/l//3hbb77zDY2HWziu5edTJYvAePmngyYvNhaIoXD0FJjBf3+y4HXoLOpb32HG7In2YG/FHLLIGuC9cZv5gTIKoGMYuuvDaVUwqVEoA8crqX1hRdJv2YFv9m1igunXsjs/Nmjdv50j4t7r17Mz17czV3P7WRjZSN3r1jIovIx8iFxhwNyJltLxVnHlvuboOkgNFVB0wF7bS/7/gUt1dYLX/2lF/QG/sh1ZjFkFEJ6obVOy9f7BUqNoJT4v6vxj09AKMSz80N0Hunk5oU3j3obRIRb3j+DpdPyuW3VRq74xVo+f95MVp4zHbdzjA+D+HKspWRu9PJQENrqrL8KWg9Dy6Fj13U7rHQ4OMA5cu3gX2BfAAp60+kFVpkvF9Jye9f6NJFSMUn6MXpjDLvP+wDO8jI+ceFOlk5cyl3n3BWXY5+opo4A//XUZp5+u5rZE7L4/z56ytjp3Y+kcBg6GqyA33YE2uutpSdtr9si0gNdGADc6ccGf1+udVHqn+fNspdM6zFUTya4vDq3kEoawx6jF5HlwN2AE3jAGHNnv/IvAZ8GgkAdcL0xZr9dFgI22VUPGGM+zCjyv/MOgepq3v3YElq63uTaudeO5umjyklz89OrFvHhBZP476c289Gfv8p1y6byxfNnkZOWxL1Uh8PqmWcMMvFbJGOsYaPuC0JHI/gbrbzudGReYyX4N1l5XS0xtMdtB/4s8EReCLKsC4E329r2ZPZeKDwZ4E4Dd4Y1F1GfdLo+raTGpCEDvYg4gXuB87E+9r1ORFYbY7ZGVNsALDHGtIvITcAPgCvtsg5jzMI4tztmzX97FnG7uT9nIwvyFrCwOGFNOcb5c0tYOi2fH67ZwSNr9/HnjQf5/Hkz+cTSKWN/OGc0iFg98rRc68Wv4xEKQmczdBy1LgSdLfbSaq27+m+3WvXbG+Dofnvbzj8eTm9v0HenD5LuvmDYFwuXF1zda581vYXLd2x+ZJleVFSMYunRnwbsMsbsARCRx4BLgZ5Ab4x5IaL+a8An49nIE2WMoXnNGrqWzOXd4Ba+P/vziW7SMbJ8bu649GSufE8Z33tmG998eiu/Xrufry4/iQvmTsDh0KGFE+J0QXq+tQxHOGwF+57A3waBdgh0RKTbrTmKBkv7G6G5GgJt9r52PsMYOnW4Yrw4eK3F6bben3BGpF0eO6/fMmS+2z6OnXbZab34jEmxBPrJQGXEdhVw+iD1bwD+FrHtE5H1WMM6dxpjnoq2k4jcCNwIUF5eHq3Kcevcto1gTQ2vX1hClieL86acF5fjjoR5k3L47Q2n88KOWr77122s/O1bnFSSxefOm8FFJ0/EqQE/MRwO8GVbS7wZYwX9oB+Cnfba37sd6IghvxOCEduBiHr+xt78UJf1FnWoy3qhLtTZO3VGPIkj4gLgsobHnG7rAhA17bKWqGl37zEGrNOdto/bnY56TheIs1/aEZF2WWXiiEh35zui7O8YN/d4Ygn00X5J1G6IiHwSWAKcHZFdboypFpFpwPMisskYs/uYAxpzH3AfWDdjY2jXkFr/9QoAq3J2cHHF5Xid3ngcdsSICOfOLuF9M4v4yzs1/PT5d7n10Q1ML9rJyrOnc8mCSfjc2mNKGiLWcE6ivjtgjHWzO9jZ7wIQsC8MEekB8yOWYL/tUADCAevTmT3poDWsFpkO+q2/lsJBOy8QJR2w97PT0R7nTQTpf6Fw9LtodF8YnH3z+5RFXHTSC+CKh+PezFgCfRUQ+YXsUqC6fyUR+QDwn8DZxpjO7nxjTLW93iMiLwKLgGMC/Uho++c/6aiYQF36ES6bedlonDIuXE4HH1k0mUsWTOJvm2u45/ldfOWJd7jzb9u56rRyPrl0ChNydIpiNUwi9hDMOHwAIByOuFj0vzAMcHExIbteyE53L8G+22aw/O79w/2OFbTbFOy3f7jfsex60doS7LKG9EZALIF+HTBTRCqAg8AK4OrICiKyCPglsNwYUxuRnwe0G2M6RaQQOAPrRu2IC7W20b5hA++cXcy0nGnMLRjgGfAxzOkQPjR/EhefMpFXd9fzq1f2ce+Lu/jFS7u5YF4JVywp46wZhbj0xq1KNQ4HOLzA2P4rfawYMtAbY4IiciuwBuvxyoeMMVtE5A5gvTFmNfBDIBP4gz2PS/djlHOAX4pIGHBgjdFvjXqiOGtf9wYEg6wpqeX8KeOnNx+NiHDGjELOmFHIgfp2fr12H0+8VcUzmw5RlOXlo4smc/mppcwqieEThkqplJO0L0zV/uhH1P3qV1zzBVj10T9yUv5JcWrd2NAZDPHC9lqeePMgL+6oJRg2zJ6QxQdPmchFJ09gpgZ9pVJKSk5q1r5hA4cnpzMxv4hZebMS3Zy487qcLD95IstPnsiR1k5Wb6zmmU01/Pj/dnLXczuZXpTBRSdP5Lw5xcwvzdWndpRKYUkZ6E1XF/5Nm3lrQZDzppyXuGmBR0lhppfrz6zg+jMrONzsZ82WQ/xt0yF+9uIu7nlhF3npbs6aWcQ5JxVx1swiirJ0XFOpVJKUgd6/dSums5Ptkx18dnKU2RiTWEm2j2uXTeXaZVM52tbFy+/W8dLOOl7eWcfqt62HpeZOzOb0afmcXpHPe6bmU5CpgV+pZJaUgb59w0YA9k9JY2HR2JnyYLTlZXi4dOFkLl04mXDYsLWmmRd31PLq7npWvXGAX72yD4AZxZk9QX9BWS5TC9KT/q8gpVJJUgZ6/5YtNOY4mTXzdNzj8RnhEeBwCCdPzuHkyTnceu5MuoJhNh1s5PW9Dbyxt4E/b6zmd68fAKxJ1+aX5rCgNJf5pTksLMulOFuf21dqvErKQN+6dTO7i8Ism7gs0U0ZszwuB6dOyefUKfncfA4EQ2F2Hm7lnapG3q5q5O3KJn7+0m5CYeuprKIsL7MnZNlLNidNyGJGcaa+qavUOJB0gT7c2Ulo3wH2LYWrJ7030c0ZN1xOB3MnZTN3UjYrTrPmGuroCrG1pomNlU1sq2lm+6FmHlm7n66g9fq50yFUFGZYQb8ok2lFGVQUWktCPpWolIoq6QJ9565dSDhMQ2kuFTkViW7OuJbmcfb0+rsFQ2H21bez/VAzOw61sK2mhXeqGnlmUw2Rr2QUZXmZVpjBtKIMphVmMqUgnbL8dErz0vQioNQoS75Av30HADnz5usNxRHgcjqYUZzJjOJMPjS/N98fCHGgoZ09da3sOdLGnro29h5p49nNhzjaHuhzjJw0N2X5aZTmplvrPOsCUJqXzoQcH9k+l/67UyqOki7QN27egN8NU+csTXRTUorP7WRWSVbUaRiOtnVxoKGdqqMdVB5tp+qolX63toUXdtTSGew7E2Ga20lJtpeSbB8TcnxMyPb1pLvXRZlePC6d40epWCRfoN/6DpWFsGDCokQ3RdnyMjzkZXhYUJZ7TJkxhiOtXT3B/3Czn0NNfg41+znc7OetA0c53NRJV+jYaWlz090UZHgoyPRSlOmlINNDQYa1Lsy08gsyPBRmecny6l8JKnUlXaAPH6iiptTJ8nE4W2UqEhGKsrwUZXkH/EC6MYaj7QEONVnBv/siUN/aRX1bJ0dauth2qJn61i6aOgJRj+FxOsjLcJOb5iEn3U1umpvcdDe56R5yutNpHnLT3b3b6R4yPE69QKhxL6kCfbi9nbSj7YSXTsLj9CS6OSpORIT8DA/5GR7mThr8S09dwTBH27s40trZcyGob+2irrWTxrYAjR1dNLYHONDQzjtV1rY/MPBHLFwOISfNTZbPRabPRZbXba19LrK8LrJ8vduZXhfZ/bazfG4yvS6da0glVFIFev/+fQCkVxznh6RV0vC4HJTYY/qx8gdCNHUEaGwP0NjeRWNHgKb23otCY0eAVn+Q1s4gLf4AlQ3ttERsh2OYANbndpDhcZHmcfauvU7S3C7SI9IZXmffOh6r3FqsvHSPE5/bic/twOdy6neF1ZCSKtAf2r4BgMJZ84eoqVQvK2g6j+vi0M0YQ0cgRIs/2Cf4t9rbLfZ2e1eI9q4g7Z0h2rtCtHUF6egK0dDWQUdXkLauEB12/vHOHO5xOvC6HfjcTrwuR5+LQHfa63ba24OU2/t7nA48LntxOnDb2147r3vbKhMd2hoHkirQH96xkVxg6lx94kaNDhEh3eMi3eOiJA7fDzfG0BkM09YZtC8O9gUiIt3RFcIfCOEPhq11wFp3Bq1099ofsOo3tIXxB0N0Brrrh+gMhgnG8qdIDI65MLjEznPaedJTFnmh8LocuBwOXE7B7XTgcggupwN399opvXlOOaau22ltuxx2eb98t10/Mt29r9ORWheomAK9iCwH7sb6wtQDxpg7+5V7gV8DpwL1wJXGmH122deAG4AQcJsxZk3cWt9P295dmEw4dfIpI3UKpUaUiPT8hVEwwucKhsIRF4vei0QgZOgKhq0lFLLXEXnddUJhOnvywgRC3ftY686IPH8gTHNHsKdep10vGAoTDBkCYWsdr4tPLHouHg7B4ZCeC0D3EpnvEOuC4XQ4cAq4HI4+daPt43QITns/h3TnOXA6wDnAeTO8Lq5ZOiXuv3XIQC8iTuBe4HysD4WvE5HV/T4JeANw1BgzQ0RWAN8HrhSRuVjfmJ0HTAL+T0RmGWNC8f4hAFJZQ2NROl6nTrur1FBcTgeZTgeZ3rHzh70xVrDvE/xDYQJhex0yBO38QMj6qyQQ6r5I2OV90r37DnTMUBhC4TAhYwjZ5+5Jhw3hfutQxNIZDBEy1v7BkCFsBq7bJ9/05kcqyvImJtADpwG7jDF7AETkMeBSIDLQXwp8004/Adwj1t9FlwKPGWM6gb0isss+3tr4NL+vzMMttC0qG4lDK6VGgYjgdgpuJ6SR/BPmGWMIGwiGw4TCJqYb+ycillcLJwOVEdtVdl7UOsaYINAEFMS4LwAicqOIrBeR9XV1dbG1PkJXVwdH5peSvUwnMlNKjQ8i1pCN12U9VTVSf13FctRodyz6X3cGqhPLvlamMfcB94H1cfAY2tWHx5PGxb8aseF/pZQat2Lp0VcBkeMhpUD1QHVExAXkAA0x7quUUmoExRLo1wEzRaRCRDxYN1dX96uzGrjOTn8MeN4YY+z8FSLiFZEKYCbwRnyarpRSKhZDDt0YY4IiciuwBuvxyoeMMVtE5A5gvTFmNfAg8Bv7ZmsD1sUAu97jWDdug8AtI/XEjVJKqejEHO9reKNgyZIlZv369YluhlJKjRsi8qYxZkm0Mp3QWymlkpwGeqWUSnIa6JVSKslpoFdKqSQ3Jm/GikgdsP8Edy8EjsSxOeOB/ubkl2q/F/Q3H68pxpiiaAVjMtAPh4isH+jOc7LS35z8Uu33gv7meNKhG6WUSnIa6JVSKsklY6C/L9ENSAD9zckv1X4v6G+Om6Qbo1dKKdVXMvbolVJKRdBAr5RSSS5pAr2ILBeRHSKyS0T+I9HtGWkiUiYiL4jINhHZIiKfT3SbRouIOEVkg4j8JdFtGQ0ikisiT4jIdvvf97JEt2mkicgX7f+uN4vIKhHxJbpN8SYiD4lIrYhsjsjLF5HnRORde50Xj3MlRaCP+ID5RcBc4Cr7w+TJLAh82RgzB1gK3JICv7nb54FtiW7EKLobeNYYMxtYQJL/dhGZDNwGLDHGnIw1PfqKxLZqRDwMLO+X9x/AP4wxM4F/2NvDlhSBnogPmBtjuoDuD5gnLWNMjTHmLTvdgvU/f9Tv8SYTESkFLgYeSHRbRoOIZAPvw/rmA8aYLmNMY2JbNSpcQJr9xbp0kvDLdMaYl7G+3xHpUuARO/0I8JF4nCtZAn3MHyFPRiIyFVgEvJ7YloyKnwBfBcKJbsgomQbUAb+yh6seEJGMRDdqJBljDgL/AxwAaoAmY8zfE9uqUVNijKkBqzMHFMfjoMkS6GP+CHmyEZFM4I/AF4wxzYluz0gSkQ8BtcaYNxPdllHkAhYDPzfGLALaiNOf82OVPS59KVABTAIyROSTiW3V+JYsgT4lP0IuIm6sIP87Y8yfEt2eUXAG8GER2Yc1PHeuiPw2sU0acVVAlTGm+6+1J7ACfzL7ALDXGFNnjAkAfwLem+A2jZbDIjIRwF7XxuOgyRLoY/mAeVIREcEat91mjLkr0e0ZDcaYrxljSo0xU7H+HT9vjEnqnp4x5hBQKSIn2VnnYX2DOZkdAJaKSLr93/l5JPkN6Airgevs9HXAn+Nx0CE/Dj4eDPQB8wQ3a6SdAVwDbBKRjXbe140xzySwTWpkfA74nd2J2QN8KsHtGVHGmNdF5AngLaynyzaQhNMhiMgq4BygUESqgNuBO4HHReQGrAveFXE5l06BoJRSyS1Zhm6UUkoNQAO9UkolOQ30SimV5DTQK6VUktNAr5RSSU4DvVJKJTkN9EopleT+f+L++lvmrnxzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(t, CA, label='A')\n",
    "ax.plot(t, CB, label='B')\n",
    "ax.plot(t, CC, label='C')\n",
    "ax.plot(t, CD, label='D')\n",
    "\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1001, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import expand, factor, Symbol, evalf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x^{4} + 4 x^{3} + 6 x^{2} + 4 x + 1$"
      ],
      "text/plain": [
       "x**4 + 4*x**3 + 6*x**2 + 4*x + 1"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Symbol('x')\n",
    "sp.expand((x+1)**4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = factor(x**4 - x**3 + x**2 - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.9375$"
      ],
      "text/plain": [
       "2.93750000000000"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f.evalf(subs={x:1.5})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(x - 1\\right) \\left(x^{3} + x + 1\\right)$"
      ],
      "text/plain": [
       "(x - 1)*(x**3 + x + 1)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}