diff --git a/README.md b/README.md
index ffd63b9..5e8ea13 100644
--- a/README.md
+++ b/README.md
@@ -21,39 +21,50 @@ how this growth and mortality are implemented can be found below.
 
 ## Quick Start
 
+This plot was generated using the following script as an example of how CoralBlox can be used:
+
+![coral blox plot](docs/assets/plain_cb.jpg)
+
+It is worth noting that in the example script below there are no environmental disturbances
+and all parameters (initial coral cover, size classes bounds, linear extensions and
+survival rates), together with the number of recruits per year, were mocked and adjusted to
+generate the above plot.
+
 Before running the model, one needs to instantiate an Array to hold the coral cover for each
 timestep, `FunctionalGroup` and `SizeClass`, and a vector of `FunctionalGroup`s:
 
 ```julia
 using CoralBlox: FunctionalGroup
 
-n_timesteps::Int64 = 75
-n_functional_groups::Int64 = 2
+n_timesteps::Int64 = 100
+n_functional_groups::Int64 = 3
 n_size_classes::Int64 = 4
 
 # Coral cover cache
 C_cover::Array{Float64, 3} = zeros(n_timesteps, n_functional_groups, n_size_classes)
 
-# Habitable area is the maximum possible cover that a location can hold
-habitable_area::Float64 = 1e7
+# Habitable area is the maximum possible cover
+habitable_area::Float64 = 1e6
 
 # Each row contains SizeClass' bounds
 size_class_bounds::Matrix{Float64} = [
-    0 0.1 0.3 0.5 0.7;
-    0 0.1 0.2 0.4 0.5
+    0 0.05 0.8 1.4 1.5;
+    0 0.05 0.5 0.9 1.0;
+    0 0.05 0.5 0.9 1.0
 ]
 
 # Mock initial coral cover
-C_cover[1,:,:] = [
-    0.07 0.08 0.06 0.05;
-    0.06 0.07 0.04 0.02;
+C_cover[1, :, :] = [
+    0.08 0.05 0.02 0.005;
+    0.07 0.06 0.03 0.007;
+    0.05 0.05 0.02 0.003
 ] .* habitable_area
 
 # Create functional groups
 functional_groups::Vector{FunctionalGroup} = FunctionalGroup.(
     eachrow(size_class_bounds[:, 1:end-1]),      # lower bounds
     eachrow(size_class_bounds[:, 2:end]),        # upper bounds
-    eachrow(initial_coral_covers)                # initial coral covers
+    eachrow(C_cover[1, :, :])                # initial coral covers
 )
 ```
 
@@ -62,46 +73,72 @@ entering the system) and matrices containing growth and survival rates for each
 `FunctionalGroup` and `SizeClass` is needed:
 
 ```julia
-using CoralBlox: timestep
-
-# Mock recruits
-recruits::Vector{Float64} = [1e5, 2e5]
+using CoralBlox: linear_extension_scale_factors, max_projected_cover,
+using CoralBlox: timestep!, coral_cover
 
 # Mock linear extensions
 linear_extensions::Matrix{Float64} = [
-    15 25 35 25;
-    10 20 25 15;
+    0.004 0.025 0.1 0.0;
+    0.003 0.007 0.005 0.0;
+    0.002 0.004 0.01 0.0
 ]
 
 # Mock survival rate
-survival_rate::Matrix{Float64} = 1 .- [
-    0.05 0.05 0.03 0.01;
-    0.1 0.08 0.04 0.02
+survival_rate::Matrix{Float64} = [
+    0.5 0.6 0.65 0.8;
+    0.6 0.7 0.8 0.8;
+    0.5 0.8 0.8 0.8
 ]
 
-for tstep::Int64 in 2:n_timesteps
-    # Linear extension scale factor
-    lin_ext_scale_factors::Vector{Float64} = linear_extension_scale_factors(
-        C_cover[tstep, :, :],
-        available_area,
-        linear_extensions,
-        size_class_bounds,
-        habitable_max_projected_cover,
-    )
+# Caclulate maximum projected cover
+habitable_max_projected_cover = max_projected_cover(
+    linear_extensions,
+    size_class_bounds,
+    habitable_area
+)
 
-    # Only re-scale if total cover is above a given threshold
-    scale_threshold = 0.7 .* available_area
-    lin_ext_scale_factors[
-        dropdims(sum(C_cover[tstep, :, :], dims=(1,2)), dims=(1,2)).< scale_threshold
-    ] .= 1
+# Only apply linear extension scale factor when cover is above the scale_threshold
+# Assuming the effects of competition for space are only relevant when population density
+# is high enough. Note that this is not a requirement of CoralBlox.
+scale_threshold = 0.9 * habitable_area
 
-    # Use scale factor to calculate growth rate
-    growth_rate::Matrix{Float64} = linear_extensions .* lin_ext_scale_factors
+# Linear extension scale factor
+local linear_extension_scale_factors::Float64
+
+for tstep::Int64 in 2:n_timesteps
+    # Apply scale factor to linear_extension when cover is above scale_threshold to account
+    # for spatial competition when population density is high
+    linear_extension_scale_factors = if sum(C_cover[tstep-1, :, :]) < scale_threshold
+        1
+    else
+        linear_extension_scale_factors(
+            C_cover[tstep-1, :, :],
+            habitable_area,
+            linear_extensions,
+            size_class_bounds,
+            habitable_max_projected_cover,
+        )
+    end
+
+    # Use scale factor to calculate growth to account for spatial competition, as explained
+    growth_rate::Matrix{Float64} = linear_extensions .* linear_extension_scale_factors
+
+    # Mock recruits proportional to each functional group's cover and available space
+    # This mock is for example purposes only
+    available_space::Float64 = habitable_area - sum(C_cover[tstep-1, :, :])
+    available_proportion::Float64 = available_space / habitable_area
+    adults_cover::Vector{Float64} = dropdims(sum(C_cover[tstep-1, :, 2:end], dims=2), dims=2)
+
+    recruits_weights::Vector{Float64} = [0.6, 0.9, 1.5]
+    availability_weight::Float64 = log(2.5, 1.5 + available_proportion)
+
+    # In reality, recruits cover for each functional group would come from a fecundity model
+    recruits::Vector{Float64} = adults_cover .* recruits_weights .* availability_weight
 
     # Perform timestep
     timestep!(
         functional_groups,
-        recruitment,
+        recruits,
         growth_rate,
         survival_rate
     )
@@ -111,7 +148,7 @@ for tstep::Int64 in 2:n_timesteps
 end
 ```
 
-The `lin_ext_scale_factors` calculated at each timestep prevents the corals from outgrowing
+The `linear_extension_scale_factors` calculated at each timestep prevents the corals from outgrowing
 the habitable area, taking into account the simultaneous growth across all functional groups
 and size classes. Details about how this scale factor is calculated and used can be
 found below.
diff --git a/docs/assets/plain_cb.jpg b/docs/assets/plain_cb.jpg
new file mode 100644
index 0000000..bba13e9
Binary files /dev/null and b/docs/assets/plain_cb.jpg differ
diff --git a/figures/cover_species_sizeclasses.png b/figures/cover_species_sizeclasses.png
deleted file mode 100644
index 731cb0c..0000000
Binary files a/figures/cover_species_sizeclasses.png and /dev/null differ
diff --git a/figures/cover_species_tot.png b/figures/cover_species_tot.png
deleted file mode 100644
index 4631e16..0000000
Binary files a/figures/cover_species_tot.png and /dev/null differ
diff --git a/src/linear_extension.jl b/src/linear_extension.jl
index cb65255..6669d4b 100644
--- a/src/linear_extension.jl
+++ b/src/linear_extension.jl
@@ -1,5 +1,6 @@
 """
-    max_projected_cover(linear_extensions::Matrix{Float64}, C_bins::Matrix{Float64}, loc_habitable_areas::Vector{Float64})
+    max_projected_cover(linear_extensions::Matrix{Float64}, bin_edges::Matrix{Float64}, habitable_area::Float64)::Float64
+    max_projected_cover(linear_extensions::Matrix{Float64}, bin_edges::Matrix{Float64}, habitable_areas::Vector{Float64})
 
 Finds the maximum projected cover for each location considering that all habitable area is
 concentrated in the in the last size class of the "worst case scenario" functional group,
@@ -7,22 +8,33 @@ and that it grows by that functional group's linear extension.
 
 # Arguments
 - `linear_extensions` : Functional Groups x Size Classes
-- `C_bins` : Functional Groups x Bin Edges
-- `loc_habitable_areas` : Locations
+- `bin_edges` : Functional Groups x Bin Edges
+- `habitable_areas` : Habitable area
 """
 function max_projected_cover(
     linear_extensions::Matrix{Float64},
     bin_edges::Matrix{Float64},
-    loc_habitable_areas::Vector{Float64}
+    habitable_area::Float64
+)::Float64
+    return _diameter_coef(linear_extensions, bin_edges) * habitable_area
+end
+function max_projected_cover(
+    linear_extensions::Matrix{Float64},
+    bin_edges::Matrix{Float64},
+    habitable_areas::Vector{Float64}
 )::Vector{Float64}
+    return [_diameter_coef(linear_extensions, bin_edges)] .* habitable_areas
+end
+
+function _diameter_coef(
+    linear_extensions::Matrix{Float64}, bin_edges::Matrix{Float64}
+)::Float64
     last_linear_extensions::Vector{Float64} = linear_extensions[:, end-1]
     last_C_bins::Vector{Float64} = bin_edges[:, end-1]
     max_diameter_idx::Int64 = findmax(last_linear_extensions .+ last_C_bins)[2]
     last_upper_bound::Float64 = last_C_bins[max_diameter_idx]
     last_linear_extension::Float64 = last_linear_extensions[max_diameter_idx]
-    diameter_coef::Float64 = ((last_upper_bound + last_linear_extension)^2) / (last_upper_bound^2)
-    max_projected_cover = [diameter_coef] .* loc_habitable_areas
-    return max_projected_cover
+    return ((last_upper_bound + last_linear_extension)^2) / (last_upper_bound^2)
 end
 
 """
@@ -31,12 +43,42 @@ end
 Adjusted linear extension.
 
 # Arguments
-- `C_cover_t` : dimensions functional_groups x size_classes x habitable_locations
+- `C_cover_t` : Cover at previous timestep
 - `loc_habitable_areas` : dimensions 1 x 1 x habitable_locations
 - `linear_extensions` : dimensions functional_groups x size_classes
 - `bin_edges` : dimensions functional_groups x bin_edges
 - `max_projected_cover` : dimensions habitable_locations
 """
+function linear_extension_scale_factors(
+    C_cover_t::AbstractMatrix{Float64},
+    habitable_area::Float64,
+    linear_extensions::AbstractMatrix{Float64},
+    bin_edges::AbstractMatrix{Float64},
+    max_projected_cover::Float64
+)::Float64
+    total_cover::Float64 = sum(C_cover_t)
+
+    # Average density for each functional group and size class
+    size_class_densities::Matrix{Float64} = _size_class_densities(C_cover_t, bin_edges)
+
+    # Projected coral cover at t+1
+    projected_cover::Float64 = _projected_cover(
+        size_class_densities, linear_extensions, bin_edges
+    )
+
+    adjusted_projected_cover::Float64 = _adjusted_projected_cover(
+        total_cover, projected_cover, max_projected_cover, habitable_area
+    )
+
+    # Solve quadratic equation
+    a::Float64 = _quadratic_coeff(size_class_densities, linear_extensions, Δd(bin_edges, 1))
+    b::Float64 = _linear_coeff(size_class_densities, linear_extensions, Δd(bin_edges, 2))
+    c::Float64 = _constant_coeff(
+        size_class_densities, adjusted_projected_cover, Δd(bin_edges, 3)
+    )
+
+    return (sqrt((b^2) - (4 * a * c)) - (b)) / (2 * a)
+end
 function linear_extension_scale_factors(
     C_cover_t::AbstractArray{Float64,3},
     loc_habitable_areas::AbstractVector{Float64},
@@ -44,12 +86,13 @@ function linear_extension_scale_factors(
     bin_edges::AbstractMatrix{Float64},
     max_projected_cover::AbstractVector{Float64}
 )::AbstractVector{Float64}
+    # Total cover for each location
     loc_C_cover::Vector{Float64} = dropdims(sum(C_cover_t, dims=(1, 2)); dims=(1, 2))
 
-    # Calculate average density for FG i and SC j. This average is not weighted in any sense
+    # Average density for each functional group, size class and location
     size_class_densities::Array{Float64,3} = _size_class_densities(C_cover_t, bin_edges)
 
-    # Calculate projected_coral_cover at t+1
+    # Projected coral cover for each location at t+1
     projected_cover::Vector{Float64} = _projected_cover(
         size_class_densities, linear_extensions, bin_edges
     )
@@ -59,9 +102,15 @@ function linear_extension_scale_factors(
     )
 
     # Solve quadratic equation
-    a::Vector{Float64} = _quadratic_coeff(size_class_densities, linear_extensions, Δd(bin_edges, 1))
-    b::Vector{Float64} = _linear_coeff(size_class_densities, linear_extensions, Δd(bin_edges, 2))
-    c::Vector{Float64} = _constant_coeff(size_class_densities, adjusted_projected_cover, Δd(bin_edges, 3))
+    a::Vector{Float64} = _quadratic_coeff(
+        size_class_densities, linear_extensions, Δd(bin_edges, 1)
+    )
+    b::Vector{Float64} = _linear_coeff(
+        size_class_densities, linear_extensions, Δd(bin_edges, 2)
+    )
+    c::Vector{Float64} = _constant_coeff(
+        size_class_densities, adjusted_projected_cover, Δd(bin_edges, 3)
+    )
 
     return (sqrt.((b .^ 2) .- (4 .* a .* c)) .- (b)) ./ (2 .* a)
 end
@@ -71,11 +120,29 @@ end
     n::Int64
 )::Matrix{Float64} = ((bin_edges[:, 2:end] .^ n) .- (bin_edges[:, 1:end-1] .^ n))
 
+_size_class_densities(
+    C_cover_t::AbstractMatrix{Float64},
+    bin_edges::AbstractMatrix{Float64}
+)::Matrix{Float64} = (12 / π) .* (C_cover_t ./ Δd(bin_edges, 3))
 _size_class_densities(
     C_cover_t::AbstractArray{Float64,3},
     bin_edges::AbstractMatrix{Float64}
 )::Array{Float64,3} = (12 / π) .* (C_cover_t ./ Δd(bin_edges, 3))
 
+function _projected_cover(
+    size_class_densities::AbstractMatrix{Float64},
+    linear_extensions::AbstractMatrix{Float64},
+    bin_edges::AbstractMatrix{Float64}
+)::Float64
+    return sum(
+        (π / 12) .*
+        size_class_densities .* (
+            (3 .* (linear_extensions .^ 2) .* Δd(bin_edges, 1)) .+
+            (3 .* linear_extensions .* Δd(bin_edges, 2)) .+
+            Δd(bin_edges, 3)
+        )
+    )
+end
 function _projected_cover(
     size_class_densities::AbstractArray{Float64,3},
     linear_extensions::AbstractMatrix{Float64},
@@ -83,15 +150,29 @@ function _projected_cover(
 )::Vector{Float64}
     return dropdims(sum(
             ((π / 12) .*
-             size_class_densities .*
-             ((3 .* (linear_extensions .^ 2) .* Δd(bin_edges, 1)) .+
-              (3 .* linear_extensions .* Δd(bin_edges, 2)) .+
-              Δd(bin_edges, 3))
+             size_class_densities .* (
+                (3 .* (linear_extensions .^ 2) .* Δd(bin_edges, 1)) .+
+                (3 .* linear_extensions .* Δd(bin_edges, 2)) .+
+                Δd(bin_edges, 3)
+            )
             ),
             dims=(2, 1)
         ); dims=(2, 1))
 end
 
+function _adjusted_projected_cover(
+    total_cover::Float64,
+    projected_cover::Float64,
+    max_projected_cover::Float64,
+    habitable_area::Float64
+)::Float64
+    return total_cover + (
+        (projected_cover - total_cover) *
+        (
+            (habitable_area - total_cover) / (max_projected_cover - total_cover)
+        )
+    )
+end
 function _adjusted_projected_cover(
     loc_cover::AbstractVector{Float64},
     projected_cover::AbstractVector{Float64},
@@ -100,10 +181,19 @@ function _adjusted_projected_cover(
 )::Vector{Float64}
     return loc_cover .+ (
         (projected_cover .- loc_cover) .*
-        ((loc_habitable_areas .- loc_cover) ./ (max_projected_cover .- loc_cover))
+        (
+            (loc_habitable_areas .- loc_cover) ./ (max_projected_cover .- loc_cover)
+        )
     )
 end
 
+function _quadratic_coeff(
+    size_class_densities::AbstractMatrix{Float64},
+    linear_extensions::AbstractMatrix{Float64},
+    Δ¹d::AbstractMatrix{Float64}
+)::Float64
+    return sum(3 .* (size_class_densities .* (linear_extensions .^ 2) .* Δ¹d))
+end
 function _quadratic_coeff(
     size_class_densities::Array{Float64,3},
     linear_extensions::AbstractMatrix{Float64},
@@ -115,6 +205,13 @@ function _quadratic_coeff(
     )
 end
 
+function _linear_coeff(
+    size_class_densities::AbstractMatrix{Float64},
+    linear_extensions::AbstractMatrix{Float64},
+    Δ²d::AbstractMatrix{Float64}
+)::Float64
+    return sum(3 .* (size_class_densities .* linear_extensions .* Δ²d))
+end
 function _linear_coeff(
     size_class_densities::Array{Float64,3},
     linear_extensions::AbstractMatrix{Float64},
@@ -126,6 +223,13 @@ function _linear_coeff(
     )
 end
 
+function _constant_coeff(
+    size_class_densities::AbstractMatrix{Float64},
+    adjusted_projected_cover::Float64,
+    Δ³d::AbstractMatrix{Float64},
+)::Float64
+    return sum((size_class_densities .* Δ³d)) - ((12 / π) * adjusted_projected_cover)
+end
 function _constant_coeff(
     size_class_densities::Array{Float64,3},
     adjusted_projected_cover::Vector{Float64},