diff --git a/docs/src/api.md b/docs/src/api.md
index a0f386f3a2..f34cbfa8e3 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -272,6 +272,7 @@ hillar
 ```@docs
 Base.convert
 ModelingToolkit.structural_simplify
+set_default_noise_scaling
 ```
 
 ## Chemistry-related functionalities
diff --git a/docs/src/faqs.md b/docs/src/faqs.md
index d9faa38a73..45d69c9982 100644
--- a/docs/src/faqs.md
+++ b/docs/src/faqs.md
@@ -59,7 +59,7 @@ plot(sol; idxs = [A, B])
 ```
 
 ## How to disable rescaling of reaction rates in rate laws?
-As explained in the [Reaction rate laws used in simulations](@ref) section, for
+As explained in the [Reaction rate laws used in simulations](@ref introduction_to_catalyst_ratelaws) section, for
 a reaction such as `k, 2X --> 0`, the generated rate law will rescale the rate
 constant, giving `k*X^2/2` instead of `k*X^2` for ODEs and `k*X*(X-1)/2` instead
 of `k*X*(X-1)` for jumps. This can be disabled when directly `convert`ing a
@@ -105,7 +105,7 @@ parametric_stoichiometry) section.
 
 ## How to set default values for initial conditions and parameters?
 How to set defaults when using the `@reaction_network` macro is described in
-more detail [here](@ref dsl_description_defaults). There are several ways to do
+more detail [here](@ref dsl_advanced_options_default_vals). There are several ways to do
 this. Using the DSL, one can use the `@species` and `@parameters` options:
 ```@example faq3
 using Catalyst
diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index 452c9cd9c6..bf63cd4e60 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -1,7 +1,7 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
 In this tutorial we provide an introduction to using Catalyst to specify
 chemical reaction networks, and then to solve ODE, jump, and SDE models
-generated from them. At the end we show what mathematical rate laws and
+generated from them[1]. At the end we show what mathematical rate laws and
 transition rate functions (i.e. intensities or propensities) are generated by
 Catalyst for ODE, SDE and jump process models.
 
@@ -151,8 +151,7 @@ underlying problem.
     variable-based parameter mappings, `u₀symmap` and `psymmap`, while when
     directly passing `repressilator` we could use either those or the
     `Symbol`-based mappings, `u₀map` and `pmap`. `Symbol`-based mappings can
-    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref),
-    see the [Basic Syntax](@ref basic_examples) section for more details.
+    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref).
 
 
 !!! note
diff --git a/docs/src/inverse_problems/behaviour_optimisation.md b/docs/src/inverse_problems/behaviour_optimisation.md
index c3b879d33d..40da2254df 100644
--- a/docs/src/inverse_problems/behaviour_optimisation.md
+++ b/docs/src/inverse_problems/behaviour_optimisation.md
@@ -1,5 +1,5 @@
 # [Optimization for non-data fitting purposes](@id behaviour_optimisation)
-In previous tutorials we have described how to use [PEtab.jl](@ref petab_parameter_fitting) and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
+In previous tutorials we have described how to use PEtab.jl and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
 
 ## [Maximising the pulse amplitude of an incoherent feed forward loop](@id behaviour_optimisation_IFFL_example)
 Incoherent feedforward loops (network motifs where a single component both activates and deactivates a downstream component) are able to generate pulses in response to step inputs[^2]. In this tutorial we will consider such an incoherent feedforward loop, attempting to generate a system with as prominent a response pulse as possible.
diff --git a/docs/src/inverse_problems/global_sensitivity_analysis.md b/docs/src/inverse_problems/global_sensitivity_analysis.md
index 91019e2383..cd65631c2f 100644
--- a/docs/src/inverse_problems/global_sensitivity_analysis.md
+++ b/docs/src/inverse_problems/global_sensitivity_analysis.md
@@ -1,6 +1,6 @@
 # [Global Sensitivity Analysis](@id global_sensitivity_analysis)
 *Global sensitivity analysis* (GSA) is used to study the sensitivity of a function's outputs with respect to its input[^1]. Within the context of chemical reaction network modelling it is primarily used for two purposes:
-- [When fitting a model's parameters to data](@ref petab_parameter_fitting), it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
+- When fitting a model's parameters to data, it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
 - [When measuring some system behaviour or property](@ref behaviour_optimisation), it can help determine which parameters influence that property. E.g. for a model of a biofuel-producing circuit in a synthetic organism, GSA could determine which system parameters have the largest impact on the total rate of biofuel production.
 
 GSA can be carried out using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package. This tutorial contains a brief introduction of how to use it for GSA on Catalyst models, with [GlobalSensitivity providing a more complete documentation](https://docs.sciml.ai/GlobalSensitivity/stable/).
@@ -52,7 +52,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0
 
 !!! note
     We should make a couple of notes about the example above:
-    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_logarithmic_scale), this is advantageous in the context of inverse problems such as this one.
+    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_log_scale), this is advantageous in the context of inverse problems such as this one.
     - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set.
     - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`.
     - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.
diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md
index 4b82c0da70..ae729d2fc2 100644
--- a/docs/src/inverse_problems/optimization_ode_param_fitting.md
+++ b/docs/src/inverse_problems/optimization_ode_param_fitting.md
@@ -129,13 +129,13 @@ If we from previous knowledge know that $kD = 0.1$, and only want to fit the val
 fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2]))
 nothing # hide
 ```
-Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_remake)). Now we create our modified loss function:
+Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_problems_remake)). Now we create our modified loss function:
 ```@example diffeq_param_estim_1 
 loss_function_fixed_kD = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); prob_generator = fixed_p_prob_generator, maxiters=10000, verbose=false, save_idxs=4)
 nothing # hide
 ```
 
-We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$* and $kP$):
+We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$ and $kP$):
 ```@example diffeq_param_estim_1 
 optprob_fixed_kD = OptimizationProblem(loss_function_fixed_kD, [1.0, 1.0])
 optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead())
@@ -160,7 +160,7 @@ nothing # hide
 corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`).
 
 ## [Parameter fitting to multiple experiments](@id optimization_parameter_fitting_multiple_experiments)
-Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
+Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref ensemble_simulations). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
 
 ## [Optimisation solver options](@id optimization_parameter_fitting_solver_options)
 Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument:
diff --git a/docs/src/model_creation/compositional_modeling.md b/docs/src/model_creation/compositional_modeling.md
index 82e9bfc9d3..1f0bfac0cd 100644
--- a/docs/src/model_creation/compositional_modeling.md
+++ b/docs/src/model_creation/compositional_modeling.md
@@ -139,7 +139,7 @@ nothing # hide
 Here we assume the user will pass in the repressor species as a ModelingToolkit
 variable, and specify a name for the network. We use Catalyst's interpolation
 ability to substitute the value of these variables into the DSL (see
-[Interpolation of Julia Variables](@ref dsl_description_interpolation_of_variables)). To make the repressilator we now make
+[Interpolation of Julia Variables](@ref dsl_advanced_options_symbolics_and_DSL_interpolation)). To make the repressilator we now make
 three genes, and then compose them together
 ```@example ex1
 t = default_t()
diff --git a/docs/src/model_creation/constraint_equations.md b/docs/src/model_creation/constraint_equations.md
index b8a08769cc..b8468dc62a 100644
--- a/docs/src/model_creation/constraint_equations.md
+++ b/docs/src/model_creation/constraint_equations.md
@@ -52,7 +52,7 @@ end
 Notice, here we interpolated the variable `V` with `$V` to ensure we use the
 same symbolic unknown variable in the `rn` as we used in building `osys`. See the
 doc section on [interpolation of variables](@ref
-dsl_description_interpolation_of_variables) for more information.
+dsl_advanced_options_symbolics_and_DSL_interpolation) for more information.
 
 We can now merge the two systems into one complete `ReactionSystem` model using
 [`ModelingToolkit.extend`](@ref):
diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md
index a7dfc694c6..ebfa8fadde 100644
--- a/docs/src/model_creation/dsl_advanced.md
+++ b/docs/src/model_creation/dsl_advanced.md
@@ -463,3 +463,77 @@ A reaction's metadata can be accessed using specific functions, e.g. `Catalyst.h
 rx = @reaction p, 0 --> X, [description="A production reaction"]
 Catalyst.getdescription(rx)
 ```
+
+## [Working with symbolic variables and the DSL](@id dsl_advanced_options_symbolics_and_DSL)
+We have previously described how Catalyst represents its models symbolically (enabling e.g. symbolic differentiation of expressions stored in models). While Catalyst utilises this for many internal operation, these symbolic representations can also be accessed and harnessed by the user. Primarily, doing so is much easier during programmatic (as opposed to DSL-based) modelling. Indeed, the section on [programmatic modelling](@ref programmatic_CRN_construction) goes into more details about symbolic representation in models, and how these can be used. It is, however, also ways to utilise these methods during DSL-based modelling. Below we briefly describe two methods for doing so.
+
+### [Using `@unpack` to extract symbolic variables from `ReactionSystem`s](@id dsl_advanced_options_symbolics_and_DSL_unpack)
+Let us consider a simple [birth-death process](@ref basic_CRN_library_bd) created using the DSL:
+```@example dsl_advanced_programmatic_unpack
+using Catalyst # hide
+bd_model = @reaction_network begin
+    (p,d), 0 <--> X
+end
+nothing # hide
+```
+Since we have not explicitly declared `p`, `d`, and `X` using `@parameters` and `@species`, we cannot represent these symbolically (only using `Symbol`s). If we wish to do so, however, we can fetch these into our current scope using the `@unpack` macro:
+```@example dsl_advanced_programmatic_unpack
+@unpack p, d, X = bd_model
+nothing # hide
+```
+This lists first the quantities we wish to fetch (does not need to be the model's full set of parameters and species), then `=`, followed by the model variable. `p`, `d` and `X` are now symbolic variables in the current scope, just as if they had been declared using `@parameters` or `@species`. We can confirm this:
+```@example dsl_advanced_programmatic_unpack
+X
+```
+Next, we can now use these to e.g. designate initial conditions and parameter values for model simulations:
+```@example dsl_advanced_programmatic_unpack
+using OrdinaryDiffEq, Plots # hide
+u0 = [X => 0.1]
+tspan = (0.0, 10.0)
+ps = [p => 1.0, d => 0.2]
+oprob = ODEProblem(bd_model, u0, tspan, ps)
+sol = solve(oprob)
+plot(sol)
+```
+
+!!! warn
+    Just like when using `@parameters` and `@species`, `@unpack` will overwrite any variables in the current scope which share name with the imported quantities.
+
+### [Interpolating variables into the DSL](@id dsl_advanced_options_symbolics_and_DSL_interpolation)
+Catalyst's DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species/stoichiometries within reactions. Using the lower-level symbolic interface we can then define symbolic variables and parameters outside of `@reaction_network`, which can then be used within expressions in the DSL. 
+
+Interpolation is carried out by pre-appending the interpolating variable with a `$`. For example, here we declare the parameters and species of a birth-death model, and interpolate these into the model:
+```@example dsl_advanced_programmatic_interpolation
+using Catalyst # hide
+t = default_t()
+@species X(t)
+@parameters p d
+bd_model = @reaction_network begin
+    ($p, $d), 0 <--> $X
+end
+```
+Additional information (such as default values or metadata) supplied to `p`, `d`, and `X` is carried through to the DSL. However, interpolation for this purpose is of limited value, as such information [can be declared within the DSL](@ref dsl_advanced_options_declaring_species_and_parameters). However, it is possible to interpolate larger algebraic expressions into the DSL, e.g. here
+```@example dsl_advanced_programmatic_interpolation
+@species X1(t) X2(t) X3(t) E(t)
+@parameters d
+d_rate = d/(1 + E)
+degradation_model = @reaction_network begin
+    $d_rate, X1 --> 0
+    $d_rate, X2 --> 0
+    $d_rate, X3 --> 0
+end
+```
+we declare an expression `d_rate`, which then can be inserted into the DSL via interpolation.
+
+It is also possible to use interpolation in combination with the `@reaction` macro. E.g. the reactions of the above network can be declared individually using
+```@example dsl_advanced_programmatic_interpolation
+rxs = [
+    @reaction $d_rate, $X1 --> 0
+    @reaction $d_rate, $X2 --> 0
+    @reaction $d_rate, $X3 --> 0
+]
+nothing # hide 
+```
+
+!!! note
+    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.
\ No newline at end of file
diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md
index ae4d204bf1..9d6249d5f3 100644
--- a/docs/src/model_simulation/ode_simulation_performance.md
+++ b/docs/src/model_simulation/ode_simulation_performance.md
@@ -8,7 +8,7 @@ Generally, this short checklist provides a quick guide for dealing with ODE perf
 1. If performance is not critical, use [the default solver choice](@ref ode_simulation_performance_solvers) and do not worry further about the issue.
 2. If improved performance would be useful, read about solver selection (both in [this tutorial](@ref ode_simulation_performance_solvers) and [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)) and then try a few different solvers to find one with good performance.
 3. If you have a large ODE (approximately 100 variables or more), try the [various options for efficient Jacobian computation](@ref ode_simulation_performance_jacobian) (noting that some are non-trivial to use, and should only be investigated if truly required).
-4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref investigating) (with preference for the former, which is easier to use).
+4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref ode_simulation_performance_parallelisation) (with preference for the former, which is easier to use).
 
 ## [Regarding stiff and non-stiff problems and solvers](@id ode_simulation_performance_stiffness)
 Generally, ODE problems can be categorised into [*stiff ODEs* and *non-stiff ODEs*](https://en.wikipedia.org/wiki/Stiff_equation). This categorisation is important due to stiff ODEs requiring specialised solvers. A common cause of failure to simulate an ODE is the use of a non-stiff solver for a stiff problem. There is no exact way to determine whether a given ODE is stiff or not, however, systems with several different time scales (e.g. a CRN with both slow and fast reactions) typically generate stiff ODEs.
diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md
index 54c84479b2..bf5fc7158e 100644
--- a/docs/src/model_simulation/simulation_introduction.md
+++ b/docs/src/model_simulation/simulation_introduction.md
@@ -275,7 +275,7 @@ nothing # hide
 ```
 If the `@default_noise_scaling` option is used, that term is only applied to reactions *without* `noise_scaling` metadata.
 
-While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the [`remake_reactionsystem`](@ref) function can be used to achieve a similar effect.
+While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the `set_default_noise_scaling` function can be used to achieve a similar effect.
 
 ## [Performing jump simulations using stochastic chemical kinetics](@id simulation_intro_jumps)
 
diff --git a/docs/src/model_simulation/simulation_plotting.md b/docs/src/model_simulation/simulation_plotting.md
index 56471fafcc..f8c40d8684 100644
--- a/docs/src/model_simulation/simulation_plotting.md
+++ b/docs/src/model_simulation/simulation_plotting.md
@@ -58,7 +58,7 @@ can be used to plot `X` only. When only a single argument is given, the vector f
 plot(sol; idxs = brusselator.X + brusselator.Y)
 ```
 
-## [Multi-plot plots](@id simulation_plotting_options)
+## [Multi-plot plots](@id simulation_plotting_options_subplots)
 It is possible to save plots in variables. These can then be used as input to the `plot` command. Here, the plot command can be used to create plots containing multiple plots (by providing multiple inputs). E.g. here we plot the concentration of `X` and `Y` in separate subplots:
 ```@example simulation_plotting
 plt_X = plot(sol; idxs = [:X])
@@ -71,7 +71,7 @@ When working with subplots, the [`layout`](https://docs.juliaplots.org/latest/la
 plot(plt_X, plt_Y; layout = (2,1), size = (700,500))
 ```
 
-## [Saving plots](@id simulation_plotting_options)
+## [Saving plots](@id simulation_plotting_options_saving)
 Once a plot has been saved to a variable, the `savefig` function can be used to save it to a file. Here we save our Brusselator plot simulation (the first argument) to a file called "saved_plot.png" (the second argument):
 ```@example simulation_plotting
 plt = plot(sol)
@@ -80,7 +80,7 @@ rm("saved_plot.png") # hide
 ```
 The plot file type is automatically determined from the extension (if none is given, a .png file is created).
 
-## [Phase-space plots](@id simulation_plotting_options)
+## [Phase-space plots](@id simulation_plotting_options_phasespace)
 By default, simulations are plotted as species concentrations over time. However, [phase space](https://en.wikipedia.org/wiki/Phase_space#:~:text=In%20dynamical%20systems%20theory%20and,point%20in%20the%20phase%20space.) plots are also possible. This is done by designating the axis arguments using the `idxs` option, but providing them as a tuple. E.g. here we plot our simulation in $X-Y$ space:
 ```@example simulation_plotting
 plot(sol; idxs = (:X, :Y))
diff --git a/docs/src/model_simulation/simulation_structure_interfacing.md b/docs/src/model_simulation/simulation_structure_interfacing.md
index eae285f296..0b45057a41 100644
--- a/docs/src/model_simulation/simulation_structure_interfacing.md
+++ b/docs/src/model_simulation/simulation_structure_interfacing.md
@@ -162,7 +162,7 @@ get_S(oprob)
 ```
 
 ## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation)
-When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref rbasic_CRN_library_two_statesef) programmatically, and use its symbolic variables to check, and update, an initial condition:
+When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref basic_CRN_library_two_states) programmatically, and use its symbolic variables to check, and update, an initial condition:
 ```@example structure_indexing_symbolic_variables
 using Catalyst
 t = default_t()
diff --git a/docs/src/steady_state_functionality/dynamical_systems.md b/docs/src/steady_state_functionality/dynamical_systems.md
index 01b21fcc57..5b844e065d 100644
--- a/docs/src/steady_state_functionality/dynamical_systems.md
+++ b/docs/src/steady_state_functionality/dynamical_systems.md
@@ -4,7 +4,7 @@ The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSyst
 ## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
 Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will typically be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 
 
-DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create an `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` and `tspan` is unused, and their values are of little importance (the only exception is than `tspan`, for implementation reason, must provide a not too small interval, we recommend minimum `(0.0, 1.0)`). 
+DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation)). As a first step, we create an `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` and `tspan` is unused, and their values are of little importance (the only exception is than `tspan`, for implementation reason, must provide a not too small interval, we recommend minimum `(0.0, 1.0)`). 
 ```@example dynamical_systems_basins
 using Catalyst
 wilhelm_model = @reaction_network begin
diff --git a/docs/src/steady_state_functionality/steady_state_stability_computation.md b/docs/src/steady_state_functionality/steady_state_stability_computation.md
index 8e43dcfe28..b23a0b0d88 100644
--- a/docs/src/steady_state_functionality/steady_state_stability_computation.md
+++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md
@@ -1,5 +1,5 @@
 # Steady state stability computation
-After system steady states have been found using [HomotopyContinuation.jl](@ref homotopy_continuation), [NonlinearSolve.jl](@ref nonlinear_solve), or other means, their stability can be computed using Catalyst's `steady_state_stability` function. Systems with conservation laws will automatically have these removed, permitting stability computation on systems with singular Jacobian.
+After system steady states have been found using [HomotopyContinuation.jl](@ref homotopy_continuation), [NonlinearSolve.jl](@ref steady_state_solving), or other means, their stability can be computed using Catalyst's `steady_state_stability` function. Systems with conservation laws will automatically have these removed, permitting stability computation on systems with singular Jacobian.
 
 !!! warn 
     Catalyst currently computes steady state stabilities using the naive approach of checking whether a system's largest eigenvalue real part is negative. While more advanced stability computation methods exist (and would be a welcome addition to Catalyst), there is no direct plans to implement these. Furthermore, Catalyst uses a tolerance `tol = 10*sqrt(eps())` to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This threshold can be changed through the `tol` keyword argument.
diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/unpublished/petab_ode_param_fitting.md
similarity index 100%
rename from docs/src/inverse_problems/petab_ode_param_fitting.md
rename to docs/unpublished/petab_ode_param_fitting.md