diff --git a/CovidDat.csv b/CovidDat.csv
new file mode 100644
index 0000000..6c98840
--- /dev/null
+++ b/CovidDat.csv
@@ -0,0 +1,61 @@
+"Time","C","R","D","I","delC","delD"
+1,1,0,0,1,NA,NA
+2,1,0,0,1,0,0
+3,2,0,0,2,1,0
+4,2,0,0,2,0,0
+5,3,0,0,3,1,0
+6,4,0,0,4,1,0
+7,4,0,0,4,0,0
+8,4,0,0,4,0,0
+9,4,0,0,4,0,0
+10,11,0,0,11,7,0
+11,12,0,0,12,1,0
+12,15,0,0,15,3,0
+13,15,0,0,15,0,0
+14,16,0,0,16,1,0
+15,19,0,0,19,3,0
+16,23,0,0,23,4,0
+17,24,1,0,23,1,0
+18,24,1,0,23,0,0
+19,25,3,0,22,1,0
+20,27,3,0,24,2,0
+21,28,3,0,25,1,0
+22,28,7,0,21,0,0
+23,28,7,0,21,0,0
+24,28,7,0,21,0,0
+25,28,9,0,19,0,0
+26,29,9,0,20,1,0
+27,30,10,0,20,1,0
+28,31,12,0,19,1,0
+29,31,12,0,19,0,0
+30,104,16,1,87,73,1
+31,204,16,2,186,100,1
+32,433,16,2,415,229,0
+33,602,18,6,578,169,4
+34,833,18,8,807,231,2
+35,977,22,10,945,144,2
+36,1261,22,12,1227,284,2
+37,1766,22,13,1731,505,1
+38,2337,22,13,2302,571,0
+39,3150,27,16,3107,813,3
+40,3736,30,17,3689,586,1
+41,4335,30,28,4277,599,11
+42,5186,30,28,5128,851,0
+43,5621,41,35,5545,435,7
+44,6088,41,35,6012,467,0
+45,6593,135,42,6416,505,7
+46,7041,135,44,6862,448,2
+47,7314,118,50,7146,273,6
+48,7478,118,53,7307,164,3
+49,7513,247,54,7212,35,1
+50,7755,288,60,7407,242,6
+51,7869,333,66,7470,114,6
+52,7979,510,66,7403,110,0
+53,8086,510,72,7504,107,6
+54,8162,510,75,7577,76,3
+55,8236,1137,75,7024,74,0
+56,8320,1407,81,6832,84,6
+57,8413,1540,84,6789,93,3
+58,8565,1540,91,6934,152,7
+59,8652,1540,94,7018,87,3
+60,8799,1540,102,7157,147,8
diff --git a/SIRD.R b/SIRD.R
new file mode 100644
index 0000000..05adec7
--- /dev/null
+++ b/SIRD.R
@@ -0,0 +1,89 @@
+### Numerical solution of the SIRD model: Korea
+{#########  Solve the SIRD model for given parameter values
+  dat <- read.table('CovidDat.csv', sep=",", header=TRUE, stringsAsFactors = 	FALSE)
+  dur <- nrow(dat)
+  
+  ### Initial conditions and parameter valuse
+  N <- 50000000  # total population = 50 mln (S. Korea)
+  S0 <- N
+  I0 <- 0.033
+  D0 <- R0 <- 0
+  beta0 <- 0.432
+  theta <- 0.234
+  gamma <- 0.0078
+  delta <- 0.00044
+  q1 <- -0.36
+  q2 <- 0.68
+  b_date <- 36.1
+  q_date <- 33
+  theta_q <- 0.3
+  
+  ### Calculate trajectory
+  model <- data.frame(Time = 1:dur, S=S0, I=I0, R=R0, D=D0)
+  pred <- data.frame(Time = 1:dur, I=NA, R=NA, D=NA, C=NA, delC=NA,delD=NA)
+  
+  for(t in 1:(dur-1)){
+    beta <- beta0*(1-1/(1+exp(-theta*(t-b_date))))
+    model$S[t+1] <- model$S[t] - beta*model$S[t]*model$I[t]/N
+    model$I[t+1] <- model$I[t] + beta*model$S[t]*model$I[t]/N - (gamma+delta)*model$I[t]
+    model$R[t+1] <- model$R[t] + gamma*model$I[t]
+    model$D[t+1] <- model$D[t] + delta*model$I[t]
+  }
+  
+  ### Data prediction
+  pred <- data.frame(Time = 1:dur, I=NA, R=NA, D=NA, C=NA, delC=NA,delD=NA)
+  for(t in 1:(dur)){
+    q <- q1 + (q2-q1)/(1+exp(-theta_q*(t-q_date)))
+    pred$I[t] <- q*model$I[t]
+    pred$R[t] <- q*model$R[t]
+    pred$D[t] <- 1*model$D[t]
+  }
+  pred$C <- pred$I + pred$R + pred$D
+  pred$delC[2:dur] <- pred$C[2:dur] - pred$C[1:(dur-1)]
+  pred$delD[2:dur] <- pred$D[2:dur] - pred$D[1:(dur-1)]
+  
+  ### Plot results and calculate the fit
+  par(mfrow=c(3,2))
+  plot(dat$Time,dat$I, type="p", pch=16, col="dark green", main="Infected (Active Cases)", xlab="",ylab="")
+  lines(pred$Time, pred$I, lwd=3, col="dark green")
+  plot(dat$Time,dat$D, type="p", pch=16, col="dark red", main="Total Deaths", xlab="",ylab="")
+  lines(pred$Time, pred$D, lwd=3, col="dark red")
+  plot(dat$Time,dat$delC, type="p", pch=16, col="dark green", main="New Cases", xlab="",ylab="")
+  lines(pred$Time, pred$delC, lwd=3, col="dark green")
+  plot(dat$Time,dat$delD, type="p", pch=16, col="dark red", main="Daily Deaths", xlab="",ylab="")
+  lines(pred$Time, pred$delD, lwd=3, col="dark red")
+  plot(dat$Time,dat$C, type="p", pch=16, col="blue", main="Total Confirmed Cases", xlab="Days from Jan. 22",ylab="")
+  lines(pred$Time, pred$C, lwd=3, col="blue")
+  plot(dat$Time,dat$R, type="p", pch=16, col="blue", main="Recovered", xlab="Days from Jan. 22",ylab="")
+  lines(pred$Time, pred$R, lwd=3, col="blue")
+  par(mfrow=c(1,1))
+  
+  predR2 <- data.frame(Var=c("I","D","delC","delD","C","R","Mean"), R2=NA)
+  pred <- pred[2:dur,]
+  dat <- dat[2:dur,]
+  predR2$R2[1] <- 1 - sum((pred$I - dat$I)^2)/sum((dat$I - mean(dat$I))^2)
+  predR2$R2[2] <- 1 - sum((pred$D - dat$D)^2)/sum((dat$D - mean(dat$D))^2)
+  predR2$R2[3] <- 1 - sum((pred$delC - dat$delC)^2)/sum((dat$delC - mean(dat$delC))^2)
+  predR2$R2[4] <- 1 - sum((pred$delD - dat$delD)^2)/sum((dat$delD - mean(dat$delD))^2)
+  predR2$R2[5] <- 1 - sum((pred$C - dat$C)^2)/sum((dat$C - mean(dat$C))^2)
+  predR2$R2[6] <- 1 - sum((pred$R - dat$R)^2)/sum((dat$R - mean(dat$R))^2)
+  predR2$R2[7] <- mean(predR2$R2[1:6])
+  print(predR2, digits=3)
+}
+
+{ ### Plot b(t) and q(t)
+  par(mfrow=c(1,2))
+  t <- 1:dur
+  beta <- beta0*(1-1/(1+exp(-theta*(t-b_date))))
+  plot(t,beta, type="l", lwd=3, col="dark red", main="(a)")
+  abline(v=seq(0,60, by=10), h=seq(0,1,by=0.05), col="grey")
+  
+  q <- q1 + (q2-q1)/(1+exp(-theta_q*(t-q_date)))
+  q[q < 0] <- 0
+  plot(t,q, type="l", lwd=3, col="dark red", main="(b)")
+  abline(v=seq(0,60, by=10), h=seq(0,1,by=0.2), col="grey")
+  par(mfrow=c(1,1))
+}
+
+ 
+ 
\ No newline at end of file
diff --git a/Turchin_2020_Covid19.pdf b/Turchin_2020_Covid19.pdf
new file mode 100644
index 0000000..3a1a813
Binary files /dev/null and b/Turchin_2020_Covid19.pdf differ