{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SEIR model in R using deSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Author*: Lloyd Chapman @LloydChapman\n",
    "\n",
    "*Date*: 2018-10-01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load deSolve library\n",
    "library(deSolve)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to return derivatives of SEIR model\n",
    "seir_ode<-function(t,Y,par){\n",
    "  S<-Y[1]\n",
    "  E<-Y[2]\n",
    "  I<-Y[3]\n",
    "  R<-Y[4]\n",
    "  \n",
    "  beta<-par[1]\n",
    "  sigma<-par[2]\n",
    "  gamma<-par[3]\n",
    "  mu<-par[4]\n",
    "  \n",
    "  dYdt<-vector(length=3)\n",
    "  dYdt[1]=mu-beta*I*S-mu*S\n",
    "  dYdt[2]=beta*I*S-(sigma+mu)*E\n",
    "  dYdt[3]=sigma*E-(gamma+mu)*I\n",
    "  \n",
    "  return(list(dYdt))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameter values\n",
    "beta<-520/365;\n",
    "sigma<-1/60;\n",
    "gamma<-1/30;\n",
    "mu<-774835/(65640000*365) # UK birth and population figures 2016\n",
    "init<-c(0.8,0.1,0.1)\n",
    "t<-seq(0,365)\n",
    "par<-c(beta,sigma,gamma,mu)\n",
    "# Solve system using lsoda\n",
    "sol<-lsoda(init,t,seir_ode,par)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAElBMVEUAAAAAAP8A/wD/AAD/\npQD////SDJfdAAAACXBIWXMAABJ0AAASdAHeZh94AAAbVUlEQVR4nO3d7XrithaAUWYy3P8t\ntyEkwWAbf2zL2tJaP6Y5fYql0/BWRjZwuQK7Xc6eALRASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBCgQEgX\nSGbDszw+nBOGgEhCggBCggBCggBCggBCggBCggBCggBCggBCggBlQ1p6N4WQSKZkSMtvTBIS\nyRQM6beftyUJiWSKhjT2Y+gQcBIhQQCndhDAZgMEsP0NAVyQhQBCggBCggDnhGTXjsbUE9LO\nzzaCMzm1gwBCggBCggBCggBl77VbupsgJM7zd87Ug4REh2ZbmYxlTtmbVg8fAh6EpjJPSDSj\nRDBThERa52Xzyq4diVRUzhMhUbtK0xkSEjWqdeGZJCQqkiueR0LifHn7+SEkzpO/nx9C4gTN\n9PNDSJTUXkF3QqKIZgu6ExIHa7ygOyFxlHZ2EhYQEgfoqKA7IRGrv4ZuhESYThu6ERIh+k3o\ni5DYq+eF6IeQ2ENDd0JiIwvRIyGxhYaeCIm1RDRCSKzgdG6KkFhIQnOExBIqekNIvCOiBYTE\nPBUtIiRmWIyWEhITbNCtISTGiGglIfFKRasJiScWoy2ExAMRbSUkfqhoOyHxxWK0i5D4pKKd\nhITFKICQuieiCELqnMUohpB65pwujJD6JaJAQuqVxSiUkLrknC6akDokonhC6o7F6AhC6oyM\njiGknnhpdBgh9UNFBxJSL2R0KCF1QkbHElIXZHQ0IXXAWd3xhNQ8GZUgpNbJqAghtU1GhQip\nZc7qihFSu2RUkJCaJaOShNQoGZUlpCY5qytNSA2SUXlCao+MTiCk1sjoFEJqi7O6kwipJTI6\njZAaIqPzCKkdOjqRkJqhozMJqRFeHp1LSE2Q0dmE1AIZnU5IDdDR+YSUn44qIKTsvDyqgpBy\nk1ElhJSajGohpMx0VA0hJaajeggpLS+PaiKkrGRUFSElJaO6CCknHVVGSCnpqDZCykhH1RFS\nPrYZKiSkdHRUIyFlI6MqCSkZHdVJSLnoqFJCSkVHtRJSIrYZ6iWkPHRUMSGlIaOaCSkJy1Hd\nhJSDjipXNKTL5XL74+0RhPRER7UrGdLl11FDNEpH1SsY0vd6dP/jiCEapaP6FQ3p59FCWkNG\nCQipejrKwKld7XSUgs2GyukoB9vfddNREi7IVk1HWQipZjpKQ0gV01Ee54Rk124JHSVST0iX\ny4JdvT/7x05DR5kkO7XrKCQdpSKkSrm/LpeyIS27HCskHaWT7M6GXkLSUTal77V7+mntEJ2E\nJKN0it/9/fzjuiH6CElH+QipPjpKyKlddXSUkc2G2thnSMn2d210lJILspXRUU5CqouOkhJS\nVXSUlZBqoqO0soXUdEk6yitZSE0vSTpKTEjVcAEpMyFVQ0eZCakWOkpNSJXQUW5CqoOOkhNS\nFXSUnZBqoKP0hFQBHeUnpPO5gNQAIZ1PRw0Q0ul01AIhnU5ILRDS2XTUBCGdTEdtENK5dNQI\nIZ1KR60Q0pl01AwhnciV2HYI6UQ6aoeQzqOjhgjpPEJqiJBOo6OWCOksOmqKkE6io7YI6Rw6\naoyQTqGj1gjpFEJqjZDOoKPmpAupiZKE1JxsITWxJOmoPUIqT0cNElJxOmqRkIoTUouaCenf\n/w6ZSzQdNamVkG4VZShJR21qJKR/g79UTEeNaiuk+ksSUqPaCOmnn9pfJ+moVY2FVPuSJKRW\nNRdS1SXpqFlNhPRYT80h6SjYJdz2qRR5SOAQoyE9/lxvSTqKFv5UFNLE/6iKkKIJafsQ70Oq\ntSQdhRPS9iHehlTrkqSjeELaPsRISP9m/2cthBRPSNuHWBBSlSXp6ABC2j7E+5CqXJJ88cQR\nhLR9iKwhnT2BJglp+xCvIb10U2FIOjqEkLYPsSCkCksS0iGEtH2IlCHp6BhC2j7EkpBqK0lH\nBxHS9iFeQhqLRkh9ENL2IRKGpKOj7Hwqvt7vLaQlf/MkOjrM3pBuLcUcUEiHE9Jh9j0VL6/H\n6Cikl5LGm6mnJB0dR0g7hkgWknuDDrQzpNeHdxzSRDH1hHT2BFoWsNkQdcBWQ6qlJB0d6el5\n8mfe2AH63bUTEj8Cnord7trlCklHh7LZsGOIhSFVUZKOjiWkHUMIiW+7Nxuc2n2bzqWCkHR0\nsL1PxZ5vEUoUko6O5qbVHUMsDen8koR0NCHtGCJNSDo6nJB2DCEkvglpxxBZQtLR8YS0Y4hB\nSLOxnFuSkI4npB1DJAlJRwUIaccQQuKbkHYMsTykM0vSUQlC2jFEipB0VMTeN/a9fPqJkEYJ\nqXERN60GHbDpkE4rSUdlCGnHEAlC0lEhQtoxhJD4JqQdQ6wJ6ZySdFSKkHYMISS+PT1P/s17\neXTXH1lcf0g6KsaKtGeIh5Led3JCSUIqRkh7hqg8JB+tWo6Q9gyxKqTyJemonE5D+nph9/o5\nsSuHqDskC1JBnd4idJvzyFbJyiEqD6nweF3r86bVe0c/P24douqQdFSSkAqGVLYkIZUkpHIh\nlV2SdFRUryFd7o9++yIpa0g6KqvPkL53ScpuNgipYZ2G9JvSriHqDUlHhXUbUsgQK0MqWJKQ\nChPSniGqDUlHpXV6Z0PQEL8hLUtESM3qPqSxV0mXy8LNiLUhlSrJzUHFCSlm+7uykMoMw6/u\nQ9o1RKUhWZDKE9KeIVaHVKYkHZXXbUgLLyNlDElHJ3h6nnzMe/Poib+1bSoHPeT+wOg7GxYH\nUqAkIZ2gzxXpt5+oe+0qCklHZ+g1pLEf1w9RY0g6OoWQCod0eElCOkWvIQWf2q2o4+CQdHSO\nPkMK32wQUu+6vWk1ZPv7p6Q1dRxakmuxJ+k2pJgh6gvpyIMzTUi7hhASX9oJaemLnh1DvKot\nJB2dpZmQFu8ebB9ixJaQjixJSGdpKKTNAy8dYkRlIenoNEJaPsSIukKyZXceIS0fYkRlIR11\nYN5qJqSjSloU0soPIz6oJAvSiXbe2VDPV1+eutmwsoyjQjrmsCwRcItQHV/rIiQdnSniXrsq\nQjrKESEdU5KQziSkXUNUFJKOTtVSSPHndS9DPBMSX5p5jbT4U/F3DPFqY0gHlGTL7lxPz5O/\n814eXdWu3cNf4qQJKfyIrNHMqd13QDlCii9JSOcS0vIhRlQTkjO7kwlp+RAj6gkp+His1FBI\nmV4jRZdkQTpbMyGduWu38la7r8dsmcwkHZ2tnZtWT7yOtCmKyJJ0dLqWQjqEkFhCSPuGuJV0\ndkg6Ol8TIX2e0Z1z93cVIdlpqICQ9s1qe0hxJemoAk2EdKT6Q7Ig1aCZkM65IFtFSEHHYQ8h\nLR9ijJC4aSOky6PNE5gdYsKOkIJKcmZXhfZC2jz+/BATKggp5Cjs1EZI1/hTupEhxuwJKaYk\nIVVBSMuHGHN6SDqqg5CWDzHmM6Qt96zeCKkdQlo+xJhbSJsPvr8kHVXiEm77VHY+5JiS6g7J\nlh0vdr+x76zt7x057C1JR7wQ0moWJF6lvdfuxJD2PZwm9RjSvpIsSIzI+Fbzk0Pa82BalfHD\nT3aHtOfBFiTGZPw4rnND2jMuzUr7NgohURMhrePMjlEpQ/osad8e9uZH64hRKV8jCYnapNy1\n+z+kzTd/f9n6aGd2jEt5HekzpH0DbA5p37A0K+WdDftD2vh4CxIThLSGjpjQ6andtpIsSExJ\nu9mwd4hNIe0dlGZl3f7eH9L6I1iQmJT1guz+D17YENLuMWlWxyGtPYQFiWn9hrT6EDpi2u7N\nht8/A9UYkgWJGVEffhK7dVckpJXH0BEzhLSQBYk5vd7ZcLPmIDpijpCWsSAxq9tbhD6tCSli\nPNrV7S1Cn1ZcShISs7q9Rehm8VGc2TEv5wXZ65/SIYUMR7uEtIiQmNd3SEtLcmbHGzlfIxUP\nKWY02pVz1y4spGUlWZB4J+d1pL2fxvVwoCX/kI54J+edDUG739dlRVqQeCtmsyFauZCWHElH\nvCWkt4eyIPHe3lO7Y0oqGNL7Q+mI96LejxQ3o6chxkWG9O5YQuI9Ib07lo5YoPddu+vbJUlI\nLCCkNwfTEUvsCemQa7HDIaYIibrsCOmQm4OGQ0yKDWnuaEJiie0h3SI6qKSyIc0dzUUkFtkT\n0u+f0SoKKXQgmiWk+aMJiUWENHs4Z3YsI6TZw+mIZYQ0dzwLEgvt2rU75HO/F82qVEjBw9As\nIX0db/SAFiSWcovQzAF1xFJCuh9w5IgWJBYT0vQRdcRiJ4S04AXVGSG9HlJILFYwpBWbEyeE\n9HpIHbFc1pD+bJzEzCGfSxISy5U8tbvnE3JqFx/S85KkI1Yo+hpp8c0Q7/6Jf1chUZWymw33\n9zDtHuKYkIYlCYkVCu/afZZUa0jDJclFJNYovv296H6iGkI6YADaVf46Ur0hDUoSEmvkvLPh\n+JCc2bGKkIbH/aYjVjknpL0XZA8L6bskCxLr1BPSmrc3HRXSz5KkI9Zxajc88L0kIbGOkJ6P\n/ElHrFT8zoaKryPdjnwrSUisVPam1aWvgc4L6bYk6Yi1yr6N4uWnjUMcGdI/IbFe0ZDGftw0\nxIEh/X9se9+sJqTXY+uI1ZzavRzbgsR6Nhte6Ij1Um5//7sKibqkvCB7cEgHfEQRrRPSs79H\nfNYXrRPSs7+HfGoejRPSs79vvuYcRgjpyW3zW0islDako0q67dlZklgpa0hHLUn3q7FCYh0h\nDX1fRFISqwhpSEhsIqSB37salMQaQhoQEtsI6dHDbXY27lhDSI8e71cVEisI6cFfIbGRkB4M\n30ChJJYT0q+/QmIrIf16fkefklhMSL9eQlISSwnp18t7zIXEUhlDuj2/40N6/awGSxJLCenH\nyIeeCImF8oYUXtLYpwcpiWXShhS+JI1+LqSQWEZI38Y/zk5JLCKku4kPKhYSiwjpburzVZXE\nEkK6ExJ7COnL9Ad+K4kFhPRFSOwipJu5b6BQEu8J6Wb2q1yUxFtCuhES+wjp05vvFlMS7wjp\n07sv6VMSbwjpuuDLLoXEG0K6LvnWWCUxT0hXIbFf3pDiSpq4X3VsTBgnpCULkpB4Q0jLQlIS\ns4QkJAIkDOn7KR0V0qKOlMQsIS0NSUnMENLSjoTEDCEtDklJTBPS8pCUxKTuQ1rRkZCYJKQ1\n/7CSmNB7SKs6UhJThLSKkBjXeUhL7lcdHRwGEocUUdLajoTEuL5DWr0gKYlxnYe0/iFCYoyQ\n1lISI4S0mpJ41XVImzpSEiOEtJ6QeCGkDZTEMyFtoSSe9BzS5o6ExDMhbaIkhjoOaUdH139K\nYkBI2wiJgXwh/T6F94W04T678WlA7pD2lbSvIyEx0G1IOxckJTHQb0g7HvsyEbonpO2UxI9e\nQwroSEn8EtIOQuJbpyGFdKQkfghpDzc4cCekXYTEFyHtoyRu+gwprCMh8UVIOymJT52G9HGz\n9eEDSuKaPKQtJX0GdL/P7uMjpCcl0VtI92z+jv3N7ZRERyH9Lj6vN37vXJiURDchPYQyvtWw\noyUh0UlIg0Ym9+w2p6Sk7nUR0seyjl7+ycWE1L0eQnqKY/4i0rYzPCX1roOQnrp4fzF2S0pK\n6lz7IT1HseSuhg0pKalvzYf0UsSy24OUxCqth7Sxo+uGlJTUs7ZDGolh+f2qq8/vlNSxdCEN\nnq1vQtrV0XV1St4v27GWQxqrYOUbKFaWtO7gNKRoSJfL5fuHzUMsD2m0gdXvRFqVkpK6VTKk\ny+W7pBIhjRaw4YOKV53fKalXBUP6rOdeUpGQxv7mprfGrklJSZ0qHNK9pLCQpksaf+5vfI+5\nknijdEhfJUWFNLMkhYa0JiUldal4SLeSjg8puKM153dK6lHZzYb7Xy6HhzTxpN/16UFLU3I5\nqUcnhHQ9PqQjOrouTklIHcp9QXYipKmn+/6Ps1MS43oKafe3Xc4ce0hJ3WkxpOMWpLmjDyip\nN+eEdOhrpKnTr5AFaeb4A0rqTD0hXR5NP3RBSBOPDPvA70V7DkrqS3undpNP8sBPzlcST/oJ\nKbCj65KUXE7qStmQ3p64vR/iXUglFqTbOG9TElJPir+NYlFK20Mq1dF1waKkpI4Uvdfu8vLT\n+iHqCen9TriS+lH+ptWnH9cOMR/S9FP7iJCUxLfsIQ1LmjnbOiSkt6d3dhx6kfzU7mlJKt3R\n+z0HJXUi+WbDMKTiC9LXoLMpKakPybe/n0KafthxIb07vxNSF5JfkF0aUtR9dhOU1L2WQjpp\nQXoztLO7LmQL6fk5uSykgxek65vTOyW1r6GQTtlqeBh9NqXjJ8CpWgpp+ngFQnq3KJWYAedp\nJ6RzF6R3U3B217j0If2UdPaCdJvDdEpKalv2kH6WpBo6us4tSkpqmpCizZzeFZwFhbUS0lxH\nZUOaK0lKzWokpLkds8IdOb3rUishTR+r9IJ0ndtzUFKr2gipqgXp02RKSmpUIyHNHOuUkKbb\nVlKbmg/phDO7u8mUyk6DIpoIqbozuy9T53cWpQalD+mzpArP7L5MpKSk9gjpWErqRAshVdyR\nknohpKMpqQsNhFR3R9MvlErPgyPlD2n2/XQ1hDSVkkWpJQ2ENHOc8y4iPVFS69KHlGBBuhlP\nqfg0OEj+kCa+2fxTNQvSp9HzOyW1IntIH7MhBU0nyFhKTu8akT6k8S9kvqlqQboZLUlKLUge\n0sdsSGHziWNRalT2kK7TIdW3IH0aPb87YR7EajmksOnEGklJSenlDun2jMwW0tj5ndO77NoN\nqc4zu7uXlOw5JJcspKeOvp6O4yXV3NHY+Z2SUssd0tdfMoYkpcZkDukjdUjO75qSOqT7X0dD\nqr+jq5Qakjiknydh3pCc3zWj1ZBydHQdSUlJKeUN6ff5lzqk1/M7JWWUNqTHZ99rSVVfRHrx\nMVyWvFBKKGtIg/+Kj4QUOp0CpJRcmyHlWpC+DM/wlJRM0pCGLyteQ4qdTiFSSqzJkHJ25Pwu\ns5whPe1ztRKSlPLKGNLLRcynkPJ2dHV+l1WykD4+PkbeF9dSSM+b4VLKIVlI44YhZdyye/KY\nkvO7FNoIaVBS/o6uwzO8f1qqXxMhXdsLyYulZNoLqYEzuy/Dl4NWpbo1GFLkRM72dIZ33kR4\no7mQmlmQ7qSUQ3shhU6kAoMzPCXVqrWQWluQvgxfK4mpQs2FFDuRajjDq1xjIbW5IH1yhle3\nNkL6KanZjm4GKYmpKm2F1HZH1+F9eFKqiZCycYZXpaZC6qGj62BZcoZXCyHl5AyvMo2EdCup\n3S27EQ+7eJalCjQVUvxE6jZo6dypdE9IuVmWKtFQSD12dH3aejh1Jl1rJqQ/vYb06Tumf9al\ns7QS0vVPxx19+lmYtHSKZkLqvKMbC9N5hNSWh82HcyfSm1ZC+vtn/Jtk++Mc7xSNhPR36iuZ\nu/TxeI539mQ6IaRG/XwmrZiKaCOkzxdIQnqlpWLaCUlJ48RURBMhfe3YCWmKlo4npE7cXzP9\nc43pGC2EdL+EJKS3LE2HaSgkJS1hZTpGAyH93NMgpMUsTdHyh/R7b5CQ1rA0hRJS19QUJX1I\njzerKmkbOe2XPaTBTd9C2kFNuyQPafjmCSHtpaathMSzjzs1rZA7pOd38ykpkJ7WEBLz9LRI\n6pBe316upKMMehLUi8whjXxMg5AONuxJUD8ShzT6cSdKKuNDUUN5Qxr/2CAhlfWhqC9lQ7p8\niRhi6uO3lHSOj86TKhnS5dfuIaY/xk5JZ+uyqYIh/fbztqS3Q8x8FZKQ6vHx0U1URUMa+3Hb\nEHOfq6qkKn18NJ1VzpBmP5/YZ64m8PEse1gpT+3efM63kvJ5Ceu7rCxtZdxsePt5+UpqwGta\nVfeVcPt7yfdOKKlVE30NIzshs3wXZJd9f4tFqUNzkY3VFhlcupCWfg/Snz9a4tWi2GYCnDpu\ntpDWfJ/Ynz9qSuFPAt89Tf1/OCeksZdJl8uCzYgN38t39m+At9b/UutTT0jBQ0BJ2U7toEpC\nggAJryNBfTLe2QDVSXmvHdQm593fUBkhQQCndhDAZgMEsP0NAVyQhQBCggBCggBCggBCggBC\nggBCggBCggBCggBCggBCggCVhgTJbHiWx4eTYuwtks032XTTzXdISMslm2+y6aab75CQlks2\n32TTTTffISEtl2y+yaabbr5DQlou2XyTTTfdfIeEtFyy+Sabbrr5DglpuWTzTTbddPMdEtJy\nyeabbLrp5jskpOWSzTfZdNPNd0hIyyWbb7LpppvvkJCWSzbfZNNNN98hIS2XbL7JpptuvkO5\nZw+VEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEEBIEODGk\nbZ/6X9z3JH+nW/HExyZZ8XSvv1/9kGO+086b8tbvzyhq+Iv++rHeiV/GJlnvdB8nnGO+M84M\n6fe/9tV6DOl7uvVO/Hd+Kaabbr5zTpvxw7/Euj1OdPjrrtTYc7LO6Y7Nsub5zhHSO9lCuiR7\nYv4sSEnmO0FI7yQL6ZLriXlJ9q93kpDeyfWbzvbEvCQLf5KQ3kn1zHzcGkkw3U/ZTkUnCOmd\nTCE9XIv5/rPm6X4R0r6Bs/wLSxTS7/WXFNO1axcz8iXHv6+Hk6XaL3Q8zCrDdL837L5fJFU/\n3zlnhpTjCvZPSPVfer9ku1Mg250YM06ccZJ/X78vO15/qsxDSBmmex2fZc3znZZwylAfIUEA\nIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEA\nIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIeWV8atWm+V3\nkZeQKuJ3kZeQKuJ3kdXl5uxZcOc3kZWQquI3kZeMKuJ3kZeQKuJ3kZeQKuJ3kZeQKuJ3kZeQ\nKuJ3kZeQKuJ3kZft74r4TSQmpHr4TUAAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEA\nIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEA\nIUEAIUEAIUGA/wAkwVE7e7A4ZAAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot solution\n",
    "plot(t,sol[,2],type=\"l\",col=\"blue\",ylim=c(0,1),ylab=\"Proportion\")\n",
    "lines(t,sol[,3],col=\"orange\")\n",
    "lines(t,sol[,4],col=\"red\")  \n",
    "lines(t,1-rowSums(sol[,2:4]),col=\"green\")\n",
    "legend(300,0.7,legend=c(\"S\",\"E\",\"I\",\"R\"),col=c(\"blue\",\"orange\",\"red\",\"green\"), lty=1, cex=0.8)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}