diff --git a/Report.qmd b/Report.qmd
index 5b22064..204aaf8 100644
--- a/Report.qmd
+++ b/Report.qmd
@@ -37,3 +37,7 @@ simulate_condor_population(0.659, 0.099, 0.099)
 ```
 
 ## Part F.  Under the assumptions of Part d, suppose that a program of double-clutching  persists for five years.  Afterwards, all captive birds are returned to the wild.   Compare the number of wild condors for this scenario to that in Part d.  
+
+```{r}
+simulate_part_f(0.659, 0.138, 0.069)
+```
diff --git a/project.R b/project.R
index 571ac54..fce92aa 100644
--- a/project.R
+++ b/project.R
@@ -38,26 +38,11 @@ captive_immature_to_adult_rate <- 1 / (captive_begin_breeding - 2)
 # an approximation which assumes an equal distribution of immature bird ages
 wild_immature_to_adult_rate <- 1 / (wild_begin_breeding - 2)
 
-# Simulates the Condor population over time using the given parameters.
-# * `rd1` - death rate of year one birds
-# * `rdI` - death rate of immature birds
-# * `rdA` - death rate of adult birds
-simulate_condor_population <- function(rd1, rdI, rdA) {
-  # captive population groups
-  c_Y0 <- c(30)
-  c_Y1 <- c(30)
-  c_I <- c(30)
-  c_A <- c(30)
-
-  # wild population groups
-  w_Y0 <- c(100)
-  w_Y1 <- c(100)
-  w_I <- c(100)
-  w_A <- c(100)
-
-  years <- seq(1, simulation_years)
-
-  for (t in 2:simulation_years) {
+simulate_condor_population_inner <- function(rd1, rdI, rdA,
+                                             c_Y0, c_Y1, c_I, c_A,
+                                             w_Y0, w_Y1, w_I, w_A,
+                                             nyears, captive_dc) {
+  for (t in (length(c_Y0) + 1):nyears) {
     # model captive population changes (year 1, immature, and adult)
     c_Y1[t] <- c_Y0[t - 1] - # last year's population
       rd1 * c_Y0[t - 1] # outflow to deaths
@@ -109,6 +94,16 @@ simulate_condor_population <- function(rd1, rdI, rdA) {
       (1 - rd_carrying) # deaths due to carrying capacity
   }
 
+  list(
+    c_Y0 = c_Y0, c_Y1 = c_Y1, c_I = c_I, c_A = c_A,
+    w_Y0 = w_Y0, w_Y1 = w_Y1, w_I = w_I, w_A = w_A
+  )
+}
+
+plot_populations <- function(c_Y0, c_Y1, c_I, c_A,
+                             w_Y0, w_Y1, w_I, w_A) {
+  years <- seq(1, simulation_years)
+
   captive_population <- c_Y0 + c_Y1 + c_I + c_A
   plot(
     years,
@@ -139,3 +134,79 @@ simulate_condor_population <- function(rd1, rdI, rdA) {
   points(years, w_A, col = "blue")
   points(years, wild_population)
 }
+
+# Simulates the Condor population over time using the given parameters.
+# * `rd1` - death rate of year one birds
+# * `rdI` - death rate of immature birds
+# * `rdA` - death rate of adult birds
+simulate_condor_population <- function(rd1, rdI, rdA) {
+  # captive population groups
+  c_Y0 <- c(30)
+  c_Y1 <- c(30)
+  c_I <- c(30)
+  c_A <- c(30)
+
+  # wild population groups
+  w_Y0 <- c(100)
+  w_Y1 <- c(100)
+  w_I <- c(100)
+  w_A <- c(100)
+
+  list2env(simulate_condor_population_inner(
+    rd1, rdI, rdA,
+    c_Y0, c_Y1, c_I, c_A,
+    w_Y0, w_Y1, w_I, w_A,
+    simulation_years, FALSE
+  ), envir = environment())
+
+  plot_populations(
+    c_Y0, c_Y1, c_I, c_A,
+    w_Y0, w_Y1, w_I, w_A
+  )
+}
+
+simulate_part_f <- function(rd1, rdI, rdA) {
+  # captive population groups
+  c_Y0 <- c(30)
+  c_Y1 <- c(30)
+  c_I <- c(30)
+  c_A <- c(30)
+
+  # wild population groups
+  w_Y0 <- c(100)
+  w_Y1 <- c(100)
+  w_I <- c(100)
+  w_A <- c(100)
+
+  # simulate first 5 years
+  list2env(simulate_condor_population_inner(
+    rd1, rdI, rdA,
+    c_Y0, c_Y1, c_I, c_A,
+    w_Y0, w_Y1, w_I, w_A,
+    5, TRUE
+  ), envir = environment())
+
+  # release all in captivity
+  w_Y0[5] <- w_Y0[5] + c_Y0[5]
+  w_Y1[5] <- w_Y1[5] + c_Y1[5]
+  w_I[5] <- w_I[5] + c_I[5]
+  w_A[5] <- w_A[5] + c_A[5]
+
+  c_Y0[5] <- 0
+  c_Y1[5] <- 0
+  c_I[5] <- 0
+  c_A[5] <- 0
+
+  # continue simulation
+  list2env(simulate_condor_population_inner(
+    rd1, rdI, rdA,
+    c_Y0, c_Y1, c_I, c_A,
+    w_Y0, w_Y1, w_I, w_A,
+    simulation_years, FALSE
+  ), envir = environment())
+
+  plot_populations(
+    c_Y0, c_Y1, c_I, c_A,
+    w_Y0, w_Y1, w_I, w_A
+  )
+}