zom

n <- 200

U <- numeric(n)
I <- numeric(n)
E <- numeric(n)
IK <-numeric(n)

U[1] <- 300   # survivors
I[1] <- 2     # infected
E[1] <- 0     # exposed 
IK[1] <- 0    # number of infected killed
P <- 3000


for (x in 1:(n)) {
  # The beginning of our for loop for N=#days
 
  E[x+1] <- E[x] + ((I[x] / P) * U[x]) - (E[x] / 2)
  I[x+1] <- I[x] + (E[x] / 2)
  U[x+1] <- U[x] - ((I[x] / P) * U[x])
  # Angry survivors mechanic
  if (U[x] < 150) {
   
    IK[x+1]=(U[x] * (1 / (( U[x]/2))))
    if (I[x] < IK[x]){
      I[x] <- 0
      break
     } else{
        I[x+1] <- I[x] + (E[x] / 2)-IK[x]
      }
  
    }
  }


# Plots

par(mfrow = c(2,2))
plot(U, type="l", col="green", main="Survivors")
plot(IK, type="l", col="yellow", main="Infectd Killed")
plot(I, type="l", col="red", main="Infected")
print(I)

##   [1]   2.000000   2.000000   2.100000   2.249933   2.429760   2.631941
##   [7]   2.854182   3.096428   3.359623   3.645204   3.954895   4.290625
##  [13]   4.654503   5.048815   5.476025   5.938787   6.439951   6.982582
##  [19]   7.569965   8.205620   8.893314   9.637075  10.441198  11.310264
##  [25]  12.249143  13.263006  14.357331  15.537908  16.810844  18.182556
##  [31]  19.659771  21.249513  22.959089  24.796071  26.768257  28.883644
##  [37]  31.150375  33.576682  36.170820  38.940983  41.895207  45.041269
##  [43]  48.386556  51.937931  55.701587  59.682873  63.886130  68.314502
##  [49]  72.969746  77.852046  82.959823  88.289557  93.835628  99.590177
##  [55] 105.542997 111.681464 117.990518 124.452685 131.048173 137.755014
##  [61] 144.549283 151.405364 156.296289 161.194116 166.027054 170.751764
##  [67] 175.338462 179.764390 184.010992 188.062727 191.906586 195.531904
##  [73] 198.930273 202.095485 205.023455 207.712127 210.161343 212.372698
##  [79] 214.349370 216.095939 217.618194 218.922951 220.017859 220.911223
##  [85] 221.611838 222.128835 222.471539 222.649348 222.671625 222.547612
##  [91] 222.286348 221.896610 221.386870 220.765251 220.039507 219.217003
##  [97] 218.304709 217.309195 216.236634 215.092812 213.883133 212.612640
## [103] 211.286021 209.907634 208.481518 207.011414 205.500781 203.952816
## [109] 202.370469 200.756457 199.113290 197.443273 195.748531 194.031019
## [115] 192.292534 190.534728 188.759121 186.967108 185.159973 183.338894
## [121] 181.504955 179.659151 177.802397 175.935534 174.059335 172.174511
## [127] 170.281714 168.381548 166.474563 164.561270 162.642136 160.717593
## [133] 158.788038 156.853838 154.915329 152.972823 151.026607 149.076944
## [139] 147.124080 145.168240 143.209633 141.248450 139.284871 137.319059
## [145] 135.351167 133.381335 131.409695 129.436367 127.461464 125.485088
## [151] 123.507337 121.528300 119.548059 117.566693 115.584272 113.600863
## [157] 111.616528 109.631324 107.645305 105.658522 103.671019 101.682841
## [163]  99.694028  97.704617  95.714643  93.724140  91.733138  89.741664
## [169]  87.749747  85.757411  83.764678  81.771571  79.778111  77.784316
## [175]  75.790203  73.795791  71.801093  69.806126  67.810902  65.815434
## [181]  63.819735  61.823816  59.827687  57.831358  55.834838  53.838137
## [187]  51.841262  49.844221  47.847020  45.849668  43.852169  41.854531
## [193]  39.856757  37.858854  35.860826  33.862678  31.864414  29.866037
## [199]  27.867551  25.868960  23.870267

plot(E, type="l", col="blue", main="Exposed")

```