diff --git a/docs/src/example/paradox_enrichment.md b/docs/src/example/paradox_enrichment.md
index b4ff64c11..414d3f969 100644
--- a/docs/src/example/paradox_enrichment.md
+++ b/docs/src/example/paradox_enrichment.md
@@ -19,7 +19,7 @@ and check that our simulations fit with analytical predictions.
 
 We consider a 2-species system with one resource (``R``) and one consumer (``C``).
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0; 1 0]); # 2 eats 1
 ```
 
@@ -29,7 +29,7 @@ Here we choose the [`ClassicResponse`](@ref)
 with a handling time (`hₜ`), an attack rate (`aᵣ`), and a hill exponent (`h`) equal to one
 to simplify analytical derivations.
 
-```@repl befwm2
+```@example befwm2
 response = ClassicResponse(foodweb, aᵣ=1, hₜ=1, h=1);
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 3fb5a1851..26129effa 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -25,7 +25,7 @@ but can also be easily customized if wanted.
     Simulations of biomass dynamics in community food webs
 - [Miele et al., 2019, PLOS](https://doi.org/10.1371/journal.pcbi.1007269) -
     Effect of non-trophic interactions on diversity-productivity relationship
-- [Kefi et al. 2018, Ecology Letters]( https://doi.org/10.1111/j.1461-0248.2011.01732.x) -
+- [Kefi et al. 2018, Ecology Letters](https://doi.org/10.1111/j.1461-0248.2011.01732.x) -
     More than a meal... integrating non-feeding interactions into food webs
 
 !!! note "Todo"
diff --git a/docs/src/man/foodwebs.md b/docs/src/man/foodwebs.md
index c3651a4a0..535a99e90 100644
--- a/docs/src/man/foodwebs.md
+++ b/docs/src/man/foodwebs.md
@@ -28,40 +28,40 @@ that will return you the corresponding `FoodWeb` object.
 using BEFWM2
 ```
 
-```@repl befwm2
+```@example befwm2
 A = [0 0 0; 1 0 0; 0 1 0] # 1 producer ⋅ 2 eats 1 ⋅ 3 eats 2
 foodweb = FoodWeb(A)
 ```
 
 We can check that adjacency matrix stored in `foodweb` corresponds to the one we provided.
 
-```@repl befwm2
+```@example befwm2
 foodweb.A
 ```
 
 As we did not use a method (e.g. the `nichemodel`) to create the foodweb,
 the method is said to be `unspecified`.
 
-```@repl befwm2
+```@example befwm2
 foodweb.method
 ```
 
 Moreover, by default:
 - the consumers are assumed to be inverterbrates
 
-```@repl befwm2
+```@example befwm2
 foodweb.metabolic_class
 ```
 
 - all body-mass are set to 1
 
-```@repl befwm2
+```@example befwm2
 foodweb.M
 ```
 
 - the `species` vector stores only species indexes as no identities were provided.
 
-```@repl befwm2
+```@example befwm2
 foodweb.species
 ```
 
@@ -76,7 +76,7 @@ it is more suited to create the foodweb using structural models.
 implements various structural models to build foodwebs.
 You can pass any of those models, with the adequate arguments, to generate foodwebs.
 
-```@repl befwm2
+```@example befwm2
 using EcologicalNetworks
 S = 20; # number of species
 C = 0.2; # connectance
@@ -100,7 +100,7 @@ so you can directly give a `UnipartiteNetwork` object to the [`FoodWeb`](@ref) m
     This function is not yet able to attribute metabolic classes or a mass to species,
     it just pass the adjacency matrix.
 
-```@repl befwm2
+```@example befwm2
 unipartite_network = EcologicalNetworks.nz_stream_foodweb()[1] # load network
 foodweb = FoodWeb(unipartite_network, method="NZ stream")
 ```
@@ -110,7 +110,7 @@ foodweb = FoodWeb(unipartite_network, method="NZ stream")
 By default all species mass are set to 1.
 However, you cange that either by giving your own body-mass vector (`M`).
 
-```@repl befwm2
+```@example befwm2
 A = [0 0 0; 1 0 0; 0 1 0]; # define adjacency matrix
 M = rand(3) # body-mass are drawn randomly in [0,1]
 foodweb = FoodWeb(A, M=M)
@@ -121,7 +121,7 @@ Or by using a consumer-resource mass ratio `Z`,
 then mass will be computed using species trophic levels (``t_l``)
 such that: ``M = Z^{t_l - 1}``.
 
-```@repl befwm2
+```@example befwm2
 A = [0 0 0; 1 0 0; 0 1 0]; # trophic levels are respectively 1, 2 and 3
 foodweb = FoodWeb(A, Z=10)
 foodweb.M
diff --git a/docs/src/man/functionalresponse.md b/docs/src/man/functionalresponse.md
index a79c7adce..5a3c2e852 100644
--- a/docs/src/man/functionalresponse.md
+++ b/docs/src/man/functionalresponse.md
@@ -41,7 +41,7 @@ The linear response and its parameters
 can be accessed by calling the [`LinearResponse`](@ref) method
 with the [`FoodWeb`](@ref) as a mandatory argument.
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]); # 1 producer ⋅ 2 eats 1 ⋅ 3 eats 2
 f = LinearResponse(foodweb);
 f.ω # preferency
@@ -51,7 +51,7 @@ f.α # consumption rate
 Above parameters take default values, but you can specify custom values.
 For instance if you want to double the attack rate of predator 3, you can do:
 
-```@repl befwm2
+```@example befwm2
 f = LinearResponse(foodweb, α=[0.0,1.0,2.0]);
 f.α # custom attack rates
 ```
@@ -62,7 +62,7 @@ corresponding to the linear functional response.
 To do so, you just need to provide the species biomass vector (`B`),
 where `B[i]` is the biomass of species ``i``.
 
-```@repl befwm2
+```@example befwm2
 f = LinearResponse(foodweb);
 B = [1, 1, 1]; # defining species biomass
 f(B) # F matrix, F[i,j] = Fᵢⱼ
@@ -118,7 +118,7 @@ The bioenergetic response and its parameters
 can be accessed by calling the [`BioenergeticResponse`](@ref) method
 with the [`FoodWeb`](@ref) as a mandatory argument.
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]); # 1 producer ⋅ 2 eats 1 ⋅ 3 eats 2
 f = BioenergeticResponse(foodweb);
 f.ω # preferency
@@ -130,7 +130,7 @@ f.h # hill exponent
 Above parameters take default values, but you can specify custom values.
 For instance if you want to set the hill exponent (`h`) to 1 instead of 2, you can do:
 
-```@repl befwm2
+```@example befwm2
 f = BioenergeticResponse(foodweb, h=1);
 f.h # custom hill exponent
 ```
@@ -141,7 +141,7 @@ corresponding to the bioenergetic functional response.
 To do so, you just need to provide the species biomass vector (`B`),
 where `B[i]` is the biomass of species ``i``.
 
-```@repl befwm2
+```@example befwm2
 f = BioenergeticResponse(foodweb);
 B = [1, 1, 1]; # defining species biomass
 f(B) # F matrix, F[i,j] = Fᵢⱼ
@@ -192,7 +192,7 @@ The classic response and its parameters
 can be accessed by calling the [`ClassicResponse`](@ref) method
 with the [`FoodWeb`](@ref) as a mandatory argument.
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]); # 1 producer ⋅ 2 eats 1 ⋅ 3 eats 2
 f = ClassicResponse(foodweb);
 f.ω # preferency
@@ -205,7 +205,7 @@ f.hₜ # handling time
 Above parameters take default values, but you can specify custom values.
 For instance if you want to set the handling time (`h`) to 0.1 instead of 1, you can do:
 
-```@repl befwm2
+```@example befwm2
 f = ClassicResponse(foodweb, hₜ=0.1);
 f.hₜ # custom handling time
 ```
@@ -216,7 +216,7 @@ corresponding to the classic functional response.
 To do so, you just need to provide the species biomass vector (`B`),
 where `B[i]` is the biomass of species ``i``.
 
-```@repl befwm2
+```@example befwm2
 f = ClassicResponse(foodweb);
 B = [1, 1, 1]; # defining species biomass
 f(B) # F matrix, F[i,j] = Fᵢⱼ
diff --git a/docs/src/man/modelparameters.md b/docs/src/man/modelparameters.md
index 76b0e1a8d..58c8a6acf 100644
--- a/docs/src/man/modelparameters.md
+++ b/docs/src/man/modelparameters.md
@@ -21,7 +21,7 @@ By default, the model parameters can be simply generated by providing only the f
 using BEFWM2
 ```
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]);
 params = ModelParameters(foodweb)
 ```
@@ -34,7 +34,7 @@ We explain how to do so for each of the 3 containers in the following sections.
 [`BioRates`](@ref) contains the species metabolic demand (`x`)
 and the species intrinsic growth rates (`r`).
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0], Z=10);
 biorates = BioRates(foodweb)
 biorates.x # metabolic demand
@@ -43,7 +43,7 @@ biorates.r # intrinsic growth rate
 
 If you want to change the default rate values you can specify your own rate vector.
 
-```@repl befwm2
+```@example befwm2
 x_custom = [1,2,3];
 r_custom = [4,5,6];
 biorates = BioRates(foodweb, x=x_custom, r=r_custom);
@@ -55,7 +55,7 @@ In case you want to attribute the same rate value for every species,
 instead of giving a vector filled with a constant,
 you can directly pass the constant.
 
-```@repl befwm2
+```@example befwm2
 biorates = BioRates(foodweb, x=1, r=2);
 biorates.x
 biorates.r
@@ -75,14 +75,14 @@ are stored in [`AllometricParams`](@ref) objects.
 You can access allometric parameters values taken from literature
 for the intrinsic growth rate and the metabolic demand.
 
-```@repl befwm2
+```@example befwm2
 DefaultGrowthParams() # default a, b for intrinsic growth rate
 DefaultMetabolismParams() # default a, b for intrinsic metabolic demand
 ```
 
 But you can also define your own allometric parameters.
 
-```@repl befwm2
+```@example befwm2
 my_allometric_params = AllometricParams(1,2,3,4,5,6)
 ```
 
@@ -92,7 +92,7 @@ which will compute the vector of the species rates:
 !!! note
     The foodweb is provided to get the species trophic levels.
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0], Z=10)
 rate = allometric_rate(foodweb, my_allometric_params) # 1 value per species
 ```
@@ -113,7 +113,7 @@ rate = allometric_rate(foodweb, my_allometric_params) # 1 value per species
 Then the `rate` vector computed with your custom allometric parameters
 can be given to [`BioRates`](@ref) as seen previously.
 
-```@repl befwm2
+```@example befwm2
 biorates = BioRates(foodweb, x=rate) # custom metabolic demand
 ```
 
@@ -125,7 +125,7 @@ As for the biological rates,
 environmental variables can be generated by only providing the foodweb,
 in which case the parameters take default values.
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 0 0 0; 0 1 0], Z=10)
 environment = Environment(foodweb)
 ```
@@ -134,13 +134,13 @@ By default carrying capacities (K) are set to 1 for producers and 0 for consumer
 as their growth term is assumed to be null.
 
 
-```@repl befwm2
+```@example befwm2
 environment.K
 ```
 
 And the temperature is set to 293.15 K.
 
-```@repl befwm2
+```@example befwm2
 environment.K
 ```
 
@@ -149,7 +149,7 @@ by providing custom values to the [`Environment`](@ref) method.
 For the temperature
 you just have to provide the scalar corresponding the the wanted temperature.
 
-```@repl befwm2
+```@example befwm2
 environment = Environment(foodweb, T=273.15);
 environment.T
 ```
@@ -159,7 +159,7 @@ And to change the carrying capacities you can either provide:
 - a scalar if you want to set all the producer carrying capacities to the
 same value
 
-```@repl befwm2
+```@example befwm2
 environment = Environment(foodweb, K=[1,2,0]); # using a vector
 environment.K
 
@@ -173,7 +173,7 @@ You can choose between 3 functional responses:
 [Linear response](@ref), [Bioenergetic response](@ref) and [Classic response](@ref).
 By default the [Bioenergetic response](@ref) response is selected.
 
-```@repl befwm2
+```@example befwm2
 params = ModelParameters(foodweb);
 params.functional_response
 ```
@@ -181,14 +181,14 @@ params.functional_response
 However, if you want you can choose another functional response by giving it as an argument.
 For the classic response, do:
 
-```@repl befwm2
+```@example befwm2
 classic_response = ClassicResponse(foodweb); # define your response
 params = ModelParameters(foodweb, functional_response = classic_response);
 params.functional_response
 ```
 
 And for the linear response, do:
-```@repl befwm2
+```@example befwm2
 linear_response = LinearResponse(foodweb); # define your response
 params = ModelParameters(foodweb, functional_response = linear_response);
 params.functional_response
diff --git a/docs/src/man/multiplexnetworks.md b/docs/src/man/multiplexnetworks.md
index 7b03e5974..9ac223d36 100644
--- a/docs/src/man/multiplexnetworks.md
+++ b/docs/src/man/multiplexnetworks.md
@@ -55,7 +55,7 @@ to generate a [`FoodWeb`](@ref),
 which is the backbone of the [`MultiplexNetwork`](@ref)
 (for more details on food web generation, see [How to generate foodwebs?](@ref)).
 
-```@repl befwm2
+```@example befwm2
 A = [0 0 0; 1 0 0; 0 1 0]; # 1 <- 2 <- 3
 foodweb = FoodWeb(A); # build food web from adjacency matrix
 ```
@@ -63,7 +63,7 @@ foodweb = FoodWeb(A); # build food web from adjacency matrix
 Now that your [`FoodWeb`](@ref) is created,
 you can directly create a [`MultiplexNetwork`](@ref) as follows:
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(foodweb);
 ```
 
@@ -71,7 +71,7 @@ As you only gave the food web to the [`MultiplexNetwork`](@ref) method
 without any additional arguments,
 the number of non-trophic links is set to zero.
 
-```@repl befwm2
+```@example befwm2
 n_links(multi_net)
 ```
 
@@ -84,19 +84,19 @@ because its values can be accessed either with the full key
 or a alias of the key.
 For instance if you want to access the trophic layer you can either use the full key...
 
-```@repl befwm2
+```@example befwm2
 multi_net.layers[:trophic]
 ```
 
 ... or an alias of the `:trophic` key (e.g. `:t`)
 
-```@repl befwm2
+```@example befwm2
 multi_net.layers[:t]
 ```
 
 Hopefully there is a cheat-sheet to know what are the aliases of each interaction.
 
-```@repl befwm2
+```@example befwm2
 interaction_names()
 ```
 
@@ -108,7 +108,7 @@ Now coming back to the trophic [`Layer`](@ref), you can see that it has three fi
 - `f`: the functional form of the non-trophic effect on the corresponding parameter
     (for more details see [Specifying non-trophic functional forms](@ref))
 
-```@repl befwm2
+```@example befwm2
 multi_net.layers[:facilitation].A # empty
 multi_net.layers[:facilitation].intensity # 1.0 by default
 multi_net.layers[:facilitation].f # (r,B_f) -> r(1+B_f) ++ of growth rate
@@ -118,7 +118,7 @@ Before explaining how to fill the non-trophic layers,
 note that the [`MultiplexNetwork`](@ref) contains, like the [`FoodWeb`](@ref),
 information about the species identities, the metabolic classes, and the body-masses:
 
-```@repl befwm2
+```@example befwm2
 multi_net.species; # species identities
 multi_net.metabolic_class; # metabolic classes
 multi_net.M; # individual body mass
@@ -139,7 +139,7 @@ and `<interaction_name>` is the full name or an alias of an interaction
 (e.g. `facilitation` or `f`).
 Thus if you want to set the connectance of the facilitation layer to `1.0` you can do:
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]); # 2 and 3 consumes 1
 net1 = MultiplexNetwork(foodweb, connectance_facilitation=1.0); # full names
 net2 = MultiplexNetwork(foodweb, C_f=1.0); # aliases
@@ -150,14 +150,14 @@ n_links(net1)[:f] == n_links(net2)[:f] == n_links(net3)[:f]
 As for the interactions, there is a cheat-sheet to retrieve the aliases
 of the [`MultiplexNetwork`](@ref) parameters.
 
-```@repl befwm2
+```@example befwm2
 multiplex_network_parameters_names()
 ```
 
 Moreover if you want to specify the same parameters
 for two or more interaction you can group them as follow
 
-```@repl befwm2
+```@example befwm2
 net1 = MultiplexNetwork(foodweb, C=(facilitation=0.5, interference=1.0));
 ```
 
@@ -165,7 +165,7 @@ Here we have set the connectance of the facilitation and interference layer
 to `0.5` and `1.0` respectively.
 This is equivalent to
 
-```@repl befwm2
+```@example befwm2
 net2 = MultiplexNetwork(foodweb, C_facilitation=0.5, C_interference=1.0);
 n_links(net1) == n_links(net2) # both ways are equivalent
 ```
@@ -173,7 +173,7 @@ n_links(net1) == n_links(net2) # both ways are equivalent
 Reversely if you want to specify two or more parameters for the same interaction
 you can group them as follow
 
-```@repl befwm2
+```@example befwm2
 net1 = MultiplexNetwork(foodweb, facilitation=(C=0.5, intensity=0.1));
 ```
 
@@ -181,7 +181,7 @@ Here we have set the connectance and the intensity of the facilitation layer
 to `0.5` and `0.1` respectively.
 This is equivalent to
 
-```@repl befwm2
+```@example befwm2
 net2 = MultiplexNetwork(foodweb, C_facilitation=0.5, intensity_facilitation=0.1);
 n_links(net1) == n_links(net2) # both ways are equivalent
 ```
@@ -210,8 +210,9 @@ For instance it does not make sense to
 define the connectance *and* the number of links of the same layer.
 Thus if you don't respect this rule an error will be thrown.
 
-```@repl befwm2
-MultiplexNetwork(foodweb, facilitation=(C=0.5, L=2))
+```jldoctest befwm2; setup = :(using BEFWM2; foodweb=FoodWeb([0 0; 0 1]))
+julia> MultiplexNetwork(foodweb, facilitation=(C=0.5, L=2))
+ERROR: ArgumentError: Ambiguous specifications for facilitation matrix adjacency: both connectance ('C' within a 'facilitation' argument) and number_of_links ('L' within a 'facilitation' argument) have been specified. Consider removing one.
 ```
 
 !!! note "Definition of connectance"
@@ -235,7 +236,7 @@ MultiplexNetwork(foodweb, facilitation=(C=0.5, L=2))
 To change the intensity `value` simply specify `intensity_<interaction_name>=value`.
 For instance if you want to set the intensity of refuge interactions to `2.0` you can do
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(foodweb, intensity_refuge=2.0);
 multi_net.layers[:refuge].intensity
 ```
@@ -243,7 +244,7 @@ multi_net.layers[:refuge].intensity
 But you can also use aliases either for the interaction (`refuge`)
 or the parameter (`intensity`), for instance
 
-```@repl befwm2
+```@example befwm2
 multi_net1 = MultiplexNetwork(foodweb, intensity_r=2.0);
 multi_net2 = MultiplexNetwork(foodweb, I_r=2.0);
 multi_net1.layers[:refuge].intensity == multi_net2.layers[:refuge].intensity == 2.0
@@ -305,7 +306,7 @@ With:
 Let's create a small [`MultiplexNetwork`](@ref) that contains competition interactions.
 To illustrate let's use the apparent competition module (1 consumer feeding on 2 plants).
 
-```@repl befwm2
+```@example befwm2
 comp_module = FoodWeb([0 0 0; 0 0 0; 1 1 0]);
 ```
 
@@ -313,7 +314,7 @@ The possible competition interactions can occur between the two producers
 i.e. species 1 and 2.
 These potential links can be accessed with:
 
-```@repl befwm2
+```@example befwm2
 A_competition_full(comp_module)
 ```
 
@@ -325,15 +326,18 @@ Now you can create a [`MultiplexNetwork`](@ref) including competition interactio
 For instance, if you want to add competition interactions as much as possible,
 you can do
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(comp_module, C_competition=1.0);
 multi_net.layers[:competition].A == A_competition_full(comp_module)
 ```
 
 Now let's say that you want to only add 1 competition link
 
-```@repl befwm2
-multi_net = MultiplexNetwork(comp_module, L_competition=1)
+```jldoctest befwm2; setup = :(using BEFWM2; comp_module = FoodWeb([0 0 0; 0 0 0; 1 1 0]))
+julia> multi_net = MultiplexNetwork(comp_module, L_competition=1)
+ERROR: ArgumentError: L should be even.
+  Evaluated: L = 1
+  Expected: L % 2 = 0
 ```
 
 You have an error because because
@@ -343,7 +347,7 @@ If you want to change that assumption because
 to add an odd number of competition links
 you can simply specify `sym_competition=false`.
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(comp_module, L_competition=1, sym_competition=false);
 multi_net.layers[:competition].A # only one link
 ```
@@ -380,13 +384,13 @@ With:
 Let's create a small [`MultiplexNetwork`](@ref) that contains competition interactions.
 We can consider the food chain module (of length 3).
 
-```@repl befwm2
+```@example befwm2
 food_chain = FoodWeb([0 0 0; 1 0 0; 0 1 0]); # 1 <- 2 <- 3
 ```
 
 Let's have look where possible facilitation interactions can occur
 
-```@repl befwm2
+```@example befwm2
 A_facilitation_full(food_chain)
 ```
 
@@ -394,7 +398,7 @@ Now you can create a [`MultiplexNetwork`](@ref) that includes facilitation links
 This time let's specify the links with and adjacency matrix.
 We want to have only one link which occurs from species 2 to species 1.
 
-```@repl befwm2
+```@example befwm2
 A_facilitation = [0 0 0; 1 0 0; 0 0 0];
 multi_net = MultiplexNetwork(food_chain, A_facilitation=A_facilitation);
 multi_net.layers[:facilitation].A == A_facilitation
@@ -454,20 +458,20 @@ Let's create a small [`MultiplexNetwork`](@ref) that contains interference inter
 We can consider the exploitative competition module
 (2 consumers feeding on the same resource).
 
-```@repl befwm2
+```@example befwm2
 exp_module = FoodWeb([0 0 0; 1 0 0; 1 0 0]);
 ```
 
 The interference links can occur between the two consumers.
 
-```@repl befwm2
+```@example befwm2
 A_interference_full(exp_module)
 ```
 
 Now you can create a [`MultiplexNetwork`](@ref)
 which includes interspecific interference links.
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(exp_module, L_i=2);
 n_links(multi_net)[:interference]
 ```
@@ -475,7 +479,7 @@ n_links(multi_net)[:interference]
 As for competition, interference is assumed by default to be symmetric
 but this can be modified.
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(exp_module, i=(sym=false, L=1));
 multi_net.layers[:interference].A
 ```
@@ -508,7 +512,7 @@ With:
 Let's create a small [`MultiplexNetwork`](@ref) that contains refuge interactions.
 To illustrate, we consider the intraguild predation module.
 
-```@repl befwm2
+```@example befwm2
 intraguild_module = FoodWeb([0 0 0; 1 0 0; 1 1 0]);
 ```
 
@@ -516,14 +520,14 @@ In this module, the producer (species 1) can possibly provide a refuge
 to the intermediate predator (species 2)
 who is eaten by the top predator (species 3).
 
-```@repl befwm2
+```@example befwm2
 A_refuge_full(intraguild_module)
 ```
 
 You can create a [`MultiplexNetwork`](@ref) that includes this refuge link.
 Moreover let's say that you also want to set the intensity of refuge interaction to `3.0`.
 
-```@repl befwm2
+```@example befwm2
 multi_net = MultiplexNetwork(intraguild_module, r=(L=1, intensity=3.0))
 ```
 
@@ -541,7 +545,7 @@ For instance if you want that growth rate function becomes quadratic i.e.
 ``r (1 + (f_0 \sum_{k \in \{\text{fac}\}} (A_\text{fac})_{ik} B_k)^2)``
 you can do
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0; 1 0]); # define a simple food web to illustrate
 custom_f(x,δx) = x*(1+δx^2) # default is x*(1+δx);
 multi_net = MultiplexNetwork(foodweb, L_f=1, functional_form_facilitation=custom_f);
diff --git a/docs/src/man/simulations.md b/docs/src/man/simulations.md
index 3387abb37..cb9cd3296 100644
--- a/docs/src/man/simulations.md
+++ b/docs/src/man/simulations.md
@@ -11,7 +11,7 @@ By default, this can be done in one line
 and with no other arguments than the [`ModelParameters`](@ref)
 and the initial biomass (`B0`).
 
-```@repl befwm2
+```@example befwm2
 foodweb = FoodWeb([0 0; 1 0]); # step 1: create the foodweb
 params = ModelParameters(foodweb); # step 2: generate model parameters
 B0 = [0.5, 0.5]; # set initial biomass
@@ -31,20 +31,20 @@ However these values depends on the simulated system.
 If your system has a fast dynamic you can decrease `tmax` value,
 by contrary if your system has slow dynamic you can increase `tmax` value.
 
-```@repl befwm2
+```@example befwm2
 solution = simulate(params, B0, tmax=50); # fast dynamic => decrease 'tmax'
 show(IOContext(stdout, :limit=>true, :displaysize=>(10, 10)), "text/plain", solution)
 ```
 
 Moreover we can note that trajectories are saved every `δt=0.25`.
 
-```@repl befwm2
+```@example befwm2
 solution.t == collect(0:0.25:50)
 ```
 
 But the timestep can be changed if you want to have lower or higher time resolution.
 
-```@repl befwm2
+```@example befwm2
 solution = simulate(params, B0, tmax=50, δt=0.5); # lower time resolution
 solution.t == collect(0:0.5:50)
 ```
diff --git a/src/macros.jl b/src/macros.jl
index a228b4ca5..fe9f89cd0 100644
--- a/src/macros.jl
+++ b/src/macros.jl
@@ -1,5 +1,5 @@
 #### Macros ####
-"Check that `var` is lower or equal than `max`."
+# Check that `var` is lower or equal than `max`.
 macro check_lower_than(var, max)
     :(
         if $(esc(var)) > $(esc(max))
@@ -10,7 +10,7 @@ macro check_lower_than(var, max)
     )
 end
 
-"Check that `var` is greater or equal than `min`."
+# Check that `var` is greater or equal than `min`.
 macro check_greater_than(var, min)
     :(
         if $(esc(var)) < $(esc(min))
@@ -21,7 +21,7 @@ macro check_greater_than(var, min)
     )
 end
 
-"Check that `var` is between `min` and `max` (bounds included)."
+# Check that `var` is between `min` and `max` (bounds included).
 macro check_is_between(var, min, max)
     :(
         if !($(esc(min)) <= $(esc(var)) <= $(esc(max)))
@@ -44,7 +44,7 @@ macro check_is_even(x)
     )
 end
 
-"Check that `var` takes one value of the vector `values`."
+# Check that `var` takes one value of the vector `values`.
 macro check_in(var, values)
     :(
         if $(esc(var)) ∉ $values
@@ -87,7 +87,7 @@ macro check_size_is_richness²(mat, S)
     )
 end
 
-"Check that `mat` has a size `size`."
+# Check that `mat` has a size `size`.
 macro check_size(mat, size)
     :(
         if size($(esc(mat))) != $size