diff --git a/docs/src/manual/nonlinmpc.md b/docs/src/manual/nonlinmpc.md
index 8fcf5bc32..2fc8a3bcd 100644
--- a/docs/src/manual/nonlinmpc.md
+++ b/docs/src/manual/nonlinmpc.md
@@ -64,10 +64,10 @@ model = setname!(NonLinModel(f, h, Ts, nu, nx, ny; p=p_model); u=vu, x=vx, y=vy)
 
 The output function ``\mathbf{h}`` converts the ``θ`` angle to degrees. Note that special
 characters like ``θ`` can be typed in the Julia REPL or VS Code by typing `\theta` and
-pressing the `<TAB>` key. Note that the parameter `p` can be of any type but use a mutable
-type like a vector of you want to modify it later. A 4th order [`RungeKutta`](@ref) method
-solves the differential equations by default. It is good practice to first simulate `model`
-using [`sim!`](@ref) as a quick sanity check:
+pressing the `<TAB>` key. The parameter `p` can be of any type but use a mutable type like a
+vector of you want to modify it later. A 4th order [`RungeKutta`](@ref) method solves the
+differential equations by default. It is good practice to first simulate `model` using
+[`sim!`](@ref) as a quick sanity check:
 
 ```@example 1
 using Plots
diff --git a/src/model/linearization.jl b/src/model/linearization.jl
index 9239919dc..b7ddc34eb 100644
--- a/src/model/linearization.jl
+++ b/src/model/linearization.jl
@@ -35,15 +35,15 @@ julia> linmodel.A
     With the nonlinear state-space model:
     ```math
     \begin{aligned}
-        \mathbf{x}(k+1) &= \mathbf{f}\Big(\mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k)\Big)  \\
-        \mathbf{y}(k)   &= \mathbf{h}\Big(\mathbf{x}(k), \mathbf{d}(k)\Big)
+        \mathbf{x}(k+1) &= \mathbf{f}\Big(\mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k), \mathbf{p}\Big)  \\
+        \mathbf{y}(k)   &= \mathbf{h}\Big(\mathbf{x}(k), \mathbf{d}(k), \mathbf{p}\Big)
     \end{aligned}
     ```
     its linearization at the operating point ``\mathbf{x_{op}, u_{op}, d_{op}}`` is:
     ```math
     \begin{aligned}
         \mathbf{x_0}(k+1) &≈ \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k)  
-                            + \mathbf{f(x_{op}, u_{op}, d_{op})} - \mathbf{x_{op}}         \\
+                            + \mathbf{f(x_{op}, u_{op}, d_{op}, p)} - \mathbf{x_{op}}         \\
         \mathbf{y_0}(k)   &≈ \mathbf{C x_0}(k) + \mathbf{D_d d_0}(k) 
     \end{aligned}
     ```
@@ -51,18 +51,18 @@ julia> linmodel.A
     documentation, and the Jacobian matrices:
     ```math
     \begin{aligned}
-        \mathbf{A}   &= \left. \frac{∂\mathbf{f(x, u, d)}}{∂\mathbf{x}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
-        \mathbf{B_u} &= \left. \frac{∂\mathbf{f(x, u, d)}}{∂\mathbf{u}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
-        \mathbf{B_d} &= \left. \frac{∂\mathbf{f(x, u, d)}}{∂\mathbf{d}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
-        \mathbf{C}   &= \left. \frac{∂\mathbf{h(x, d)}}{∂\mathbf{x}}    \right|_{\mathbf{x=x_{op},\, d=d_{op}}}             \\
-        \mathbf{D_d} &= \left. \frac{∂\mathbf{h(x, d)}}{∂\mathbf{d}}    \right|_{\mathbf{x=x_{op},\, d=d_{op}}}
+        \mathbf{A}   &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{x}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
+        \mathbf{B_u} &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{u}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
+        \mathbf{B_d} &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{d}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\
+        \mathbf{C}   &= \left. \frac{∂\mathbf{h(x, d, p)}}{∂\mathbf{x}}    \right|_{\mathbf{x=x_{op},\, d=d_{op}}}             \\
+        \mathbf{D_d} &= \left. \frac{∂\mathbf{h(x, d, p)}}{∂\mathbf{d}}    \right|_{\mathbf{x=x_{op},\, d=d_{op}}}
     \end{aligned}
     ```
     Following [`setop!`](@ref) notation, we find:
     ```math
     \begin{aligned}
-        \mathbf{f_{op}} &= \mathbf{f(x_{op}, u_{op}, d_{op})} \\
-        \mathbf{y_{op}} &= \mathbf{h(x_{op}, d_{op})}
+        \mathbf{f_{op}} &= \mathbf{f(x_{op}, u_{op}, d_{op}, p)} \\
+        \mathbf{y_{op}} &= \mathbf{h(x_{op}, d_{op}, p)}
     \end{aligned}
     ```
     Notice that ``\mathbf{f_{op} - x_{op} = 0}`` if the point is an equilibrium. The 
diff --git a/src/sim_model.jl b/src/sim_model.jl
index 4a49529ca..242bf5d37 100644
--- a/src/sim_model.jl
+++ b/src/sim_model.jl
@@ -64,16 +64,16 @@ The state-space description of [`LinModel`](@ref) around the operating points is
 ```math
 \begin{aligned}
     \mathbf{x_0}(k+1) &= \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k) 
-                         + \mathbf{f_{op}}   - \mathbf{x_{op}}                                 \\
+                         + \mathbf{f_{op}}   - \mathbf{x_{op}}                                              \\
     \mathbf{y_0}(k)   &= \mathbf{C x_0}(k) + \mathbf{D_d d_0}(k) 
 \end{aligned}
 ```
 and, for [`NonLinModel`](@ref):
 ```math
 \begin{aligned}
-    \mathbf{x_0}(k+1) &= \mathbf{f}\Big(\mathbf{x_0}(k), \mathbf{u_0}(k), \mathbf{d_0}(k)\Big) 
-                            + \mathbf{f_{op}} - \mathbf{x_{op}}                                \\
-    \mathbf{y_0}(k)   &= \mathbf{h}\Big(\mathbf{x_0}(k), \mathbf{d_0}(k)\Big)
+    \mathbf{x_0}(k+1) &= \mathbf{f}\Big(\mathbf{x_0}(k), \mathbf{u_0}(k), \mathbf{d_0}(k), \mathbf{p}\Big) 
+                            + \mathbf{f_{op}} - \mathbf{x_{op}}                                             \\
+    \mathbf{y_0}(k)   &= \mathbf{h}\Big(\mathbf{x_0}(k), \mathbf{d_0}(k), \mathbf{p}\Big)
 \end{aligned}
 ```
 The state `xop` and the additional `fop` operating points are frequently zero e.g.: when