diff --git a/Dean Barcelona model.R b/Dean Barcelona model.R
new file mode 100644
index 0000000..fbd3191
--- /dev/null
+++ b/Dean Barcelona model.R	
@@ -0,0 +1,27 @@
+library(rstan)
+
+barcelona_1490 <- c(1,0,1,1,0,1,5,3,1,0,1,1,2,3,5,0,6,3,6,3,8,1,5,2,
+                    1,1,2,2,2,5,7,12,4,3,5,3,8,5,8,8,6,12,11,22,15,14,
+                    24,14,15,20,20,13,11,25,28,30,24,28,42,24,32,24,27,
+                    31,34,33,29,31,38,40,42,38,53,44,66,52,53,56,63,49,
+                    60,57,65,55,55,47,67,62,65,57,47,46,62,54,52,48,49,
+                    64,46,67,52,50,56,46,41,38,36,39,31,32,41,25,32,35,
+                    36,36,33,26,42,31,19,27,23,22,15,24,32,19,10,16,12,
+                    15,14,13,12,13,12,6,12,15,5,9,3,5,12,6,7,3,3,3,3,2,
+                    3,3,0,3,2,3,3,1,1,4,2,3,0,2,3,2,0,1,1,4,1,2,2,1,1,2,
+                    0,1,1,2)
+
+init <- 25000.0
+
+dat_list <- list(N = length(barcelona_1490),
+                 init = init,
+                 D = barcelona_1490)
+
+n = 1000
+
+model <- stan(file = "C:\\Users\\bios3\\OneDrive\\Documents\\humanplague.stan",
+               data = c(dat_list),iter = n)
+
+model
+post <- rstan::extract(model)
+hist(post$beta)
diff --git a/humanplague.stan b/humanplague.stan
new file mode 100644
index 0000000..2aa31fa
--- /dev/null
+++ b/humanplague.stan
@@ -0,0 +1,49 @@
+functions {
+  array[] real dpop_dt(real t,
+                 array[] real init,
+                 array[] real beta,
+                 array[] real x_r, array[] int x_i) {
+    real H = init[1];
+    real d = beta[1];
+    real dH_dt = -d * H;
+    return { dH_dt };
+  }
+}
+data {
+  int<lower=0> N;
+  array[N] int<lower=0> D;
+  real<lower=0> init;
+}
+transformed data {
+  array[N+1] real times_measured;
+  for (i in 1:(N+1)) {
+    times_measured[i] = i;
+  }
+}
+parameters {
+  array[1] real beta;
+}
+transformed parameters {
+  array[N+1] real pop;                       // trajectory of H
+  array[1] real y0;                          // initial state
+  array[N+1, 1] real ode_sol;          // solver output
+  y0[1] = init;
+  ode_sol = integrate_ode_rk45(
+    dpop_dt, y0, 0, times_measured, beta,
+    rep_array(0.0,0), rep_array(0,0),
+    1e-5, 1e-3, 5e2
+  );
+  for (t in 1:(N+1)) {
+    pop[t] = ode_sol[t,1];             // extract single state variable
+  }
+}
+model {
+  array[N] real Deaths;
+
+  beta ~ normal(0, 0.25);
+
+  for (t in 2:(N+1)) {
+    Deaths[t-1] = fmax(pop[t-1] - pop[t], 0);
+    D[t-1] ~ poisson(Deaths[t-1]);
+  }
+}