diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4b5c4fa2..0005d211 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ CellListMap.jl Changelog
 
 Version 0.10.0-DEV
 --------------
+- ![FEATURE][badge-feature] Automatic tracking of position mutations and updating of cell lists. 
 - ![FEATURE][badge-feature] `pairwise!(; reset=[true(default)/false]`: the optional `reset` keyword of `pairwise!`, when set to `false`, avoids resetting the initial value of `output` to `zero(typeof(output))`.
 - ![FEATURE][badge-feature] `update!(sys::AbstractParticleSystem;)` as cleaner and more convenient way to update system properties. 
 - ![FEATURE][badge-feature] `reduce_output!` is exposed API for advanced custom reductions. 
@@ -34,6 +35,7 @@ Version 0.10.0-DEV
 - ![INFO][badge-info] Internal implementation of neighborlist moved to ParticleSystem interface.
 - ![INFO][badge-info] The internal representation of coordinantes is done with `ParticleSystemPositions` array type - which is not a subtype of abstract array. 
 - ![INFO][badge-info] Add performance test.
+- ![INFO][badge-info] Update examples.
 - ![INFO][badge-info] Default name for output in ParticleSystem is `default_output_name`, but it can be accessed though `sys.output`, as before.
 - ![BUGFIX][badge-bugfix] Fix cross-computation (`pairwise!(f, x, sys)`) with `NonPeriodicCell` failing to find pairs when particle coordinates span negative values or large coordinate ranges / do not modify input coordinates in `NonPeriodicCell`. This bug was never released, it was created and fixed in the development version.
 - ![BUGFIX][badge-bugfix] If the number of particles is reduced by updating and becomes smaller than nbatches, there was a crash. Fixed. Also fixed automatic updating of nbatches when first set of positions change (bug never released).
diff --git a/docs/Project.toml b/docs/Project.toml
index 97a63f13..26881625 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,5 +7,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/docs/make.jl b/docs/make.jl
index 884b398f..873a019a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,6 +3,7 @@ using Documenter
 using CellListMap
 ENV["LINES"] = 10
 ENV["COLUMNS"] = 120
+ENV["GKSwstype"] = "nul"
 makedocs(
     modules = [CellListMap],
     sitename = "CellListMap.jl",
diff --git a/docs/src/ParticleSystem/examples.md b/docs/src/ParticleSystem/examples.md
index 6923dcab..576c0a6f 100644
--- a/docs/src/ParticleSystem/examples.md
+++ b/docs/src/ParticleSystem/examples.md
@@ -8,7 +8,7 @@ and run directly.
 In this example, a simple potential energy defined as the sum of the
 inverse of the distance between the particles is computed.
 
-```julia
+```@example
 using CellListMap
 using StaticArrays
 system = ParticleSystem(
@@ -19,6 +19,7 @@ system = ParticleSystem(
     output_name = :energy
 )
 pairwise!((pair, energy) -> energy += 1 / pair.d, system)
+println(" Energy: $(system.energy)")
 ```
 
 ## Force computation
@@ -26,7 +27,7 @@ pairwise!((pair, energy) -> energy += 1 / pair.d, system)
 Here we compute the force vector associated to the potential energy
 function of the previous example.
 
-```julia
+```@example
 using CellListMap
 using StaticArrays
 positions = rand(SVector{3,Float64},1000)
@@ -45,6 +46,10 @@ function update_forces!(pair, forces)
     return forces
 end
 pairwise!(update_forces!, system)
+println(""" 
+    Force on particle 1: $(system.forces[1]))
+    Force on particle 2: $(system.forces[1]))
+""")
 ```
 
 ## Energy and forces
@@ -52,10 +57,10 @@ pairwise!(update_forces!, system)
 In this example, the potential energy and the forces are computed in a single
 run, and a custom data structure is defined to store both values.
 
-```julia
+```@example 
 using CellListMap
 using StaticArrays
-# Define custom type
+# Define custom type to store energy and force vector
 mutable struct EnergyAndForces
     energy::Float64
     forces::Vector{SVector{3,Float64}}
@@ -63,17 +68,15 @@ end
 # Custom copy, reset and reducer functions
 import CellListMap: copy_output, reset_output!, reducer
 copy_output(x::EnergyAndForces) = EnergyAndForces(copy(x.energy), copy(x.forces))
-function reset_output!(output::EnergyAndForces)
-    output.energy = 0.0
-    for i in eachindex(output.forces)
-        output.forces[i] = SVector(0.0, 0.0, 0.0)
-    end
-    return output
+function reset_output!(x::EnergyAndForces)
+    x.energy = 0.0
+    fill!(x.forces, SVector(0.0, 0.0, 0.0))
+    return x
 end
 function reducer(x::EnergyAndForces, y::EnergyAndForces)
-    e_tot = x.energy + y.energy
+    x.energy += y.energy
     x.forces .+= y.forces
-    return EnergyAndForces(e_tot, x.forces)
+    return x
 end
 # Function that updates energy and forces for each pair
 function energy_and_forces!(pair, output::EnergyAndForces)
@@ -95,6 +98,11 @@ system = ParticleSystem(
 )
 # Compute energy and forces
 pairwise!(energy_and_forces!, system)
+# Print some results
+println("""
+    Energy: $(system.energy_and_forces.energy)
+    Force on particle 1: $(system.energy_and_forces.forces[1])
+""")
 ```
 
 ## Two sets of particles
@@ -102,7 +110,7 @@ pairwise!(energy_and_forces!, system)
 In this example we illustrate the interface for the computation of properties
 of two sets of particles, by computing the minimum distance between the two sets.
 
-```julia
+```@example
 using CellListMap
 using StaticArrays
 # Custom structure to store the minimum distance pair
@@ -114,15 +122,13 @@ end
 # Function that updates the minimum distance found
 function minimum_distance(pair, md)
     (; i, j, d) = pair
-    if d < md.d
-        md = MinimumDistance(i, j, d)
-    end
+    md = d < md.d ? MinimumDistance(i, j, d) : md
     return md
 end
 # Define appropriate methods for copy, reset and reduce
 import CellListMap: copy_output, reset_output!, reducer!
 copy_output(md::MinimumDistance) = md
-reset_output!(md::MinimumDistance) = MinimumDistance(0, 0, +Inf)
+reset_output!(::MinimumDistance) = MinimumDistance(0, 0, +Inf)
 reducer!(md1::MinimumDistance, md2::MinimumDistance) = md1.d < md2.d ? md1 : md2
 # Build system
 xpositions = rand(SVector{3,Float64},100);
@@ -145,14 +151,16 @@ In this example, a complete particle simulation is illustrated, with a simple po
 simulation with:
 
 ```julia-repl
-julia> trajectory = simulate(200)
+julia> trajectory = simulate(200) # 200 particles
 
 julia> animate(trajectory)
 ```
 
-One important characteristic of this example is that the `system` is built outside the function that performs the simulation. This is done because the construction of the system is type-unstable (it is dimension, geometry and output-type dependent). Adding a function barrier avoids type-instabilities to propagate to the simulation causing possible performance problems.
+Forces are computed similarly to what is done in the [Force computation](@ref) example. Additionally, we use
+the public (but not exported) `CellListMap.wrap_relative_to` function, to keep the saved coordinates within
+the minimal periodic coordinates relative to the origin, to produce a nice animation at the end.
 
-```julia
+```@example ex_simulation
 using StaticArrays
 using CellListMap
 import CellListMap.wrap_relative_to
@@ -211,7 +219,11 @@ function simulate(N; nsteps::Int=100, isave=1)
     end
     return trajectory
 end
+```
 
+The following function will create an animation from the resulting trajectory of the simulation:
+
+```@example ex_simulation
 using Plots
 function animate(trajectory)
     anim = @animate for step in trajectory
@@ -221,8 +233,17 @@ function animate(trajectory)
             lims=(-0.5, 0.5),
             aspect_ratio=1,
             framestyle=:box,
+            size=(400,400),
+            ticks=nothing,
         )
     end
     gif(anim, "simulation.gif", fps=10)
 end
 ```
+
+Now, running the simulation and creating an animation:
+
+```@example ex_simulation
+trajectory = simulate(200) # 200 particles
+animate(trajectory)
+```
diff --git a/docs/src/migrating.md b/docs/src/migrating.md
index 6fd84ab2..7d1e6900 100644
--- a/docs/src/migrating.md
+++ b/docs/src/migrating.md
@@ -66,10 +66,10 @@ The non-mutating `pairwise` function has been removed. This change emphasizes th
 
 ```julia
 # Before (0.9.x)
-u = pairwise(f, system)
+u = map_pairwise(f, system)
 
 # After (0.10.0)
-u = pairwise!(f, system)  # always use the mutating version
+u = pairwise!(f, system) # always use the mutating syntax
 ```
 
 ### Simplified mapped function signature