diff --git a/Project.toml b/Project.toml
index 40d25378b..5608b474c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "ControlSystems"
 uuid = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.13.0"
+version = "1.13.1"
 
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
diff --git a/docs/src/examples/example.md b/docs/src/examples/example.md
index e52c4802f..31a9fcb28 100644
--- a/docs/src/examples/example.md
+++ b/docs/src/examples/example.md
@@ -209,81 +209,6 @@ save_docs_plot("ppgofplot.svg"); # hide
 ![](../../plots/ppgofplot.svg)
 
 
-## Stability boundary for PID controllers
-The stability boundary, i.e., the surface of PID parameters where the transfer function ``P(s)C(s)`` equals -1, can be plotted with the command [`stabregionPID`](@ref). The process can be given in function form or as a regular LTIsystem.
-
-```jldoctest; output = false
-P1 = s -> exp(-sqrt(s))
-doplot = true
-form = :parallel
-kp, ki, f1 = stabregionPID(P1,exp10.(range(-5, stop=1, length=1000)); doplot, form); f1
-P2 = s -> 100*(s+6).^2. /(s.*(s+1).^2. *(s+50).^2)
-kp, ki, f2 = stabregionPID(P2,exp10.(range(-5, stop=2, length=1000)); doplot, form); f2
-P3 = tf(1,[1,1])^4
-kp, ki, f3 = stabregionPID(P3,exp10.(range(-5, stop=0, length=1000)); doplot, form); f3
-
-save_docs_plot(f1, "stab1.svg") # hide
-save_docs_plot(f2, "stab2.svg") # hide
-save_docs_plot(f3, "stab3.svg"); # hide
-
-# output
-
-```
-![](../../plots/stab1.svg)
-![](../../plots/stab2.svg)
-![](../../plots/stab3.svg)
-
-
-## PID plots
-This example utilizes the function [`pidplots`](@ref), which accepts vectors of PID-parameters and produces relevant plots. The task is to take a system with bandwidth 1 rad/s and produce a closed-loop system with bandwidth 0.1 rad/s. If one is not careful and proceed with pole placement, one easily get a system with very poor robustness.
-```jldoctest PIDPLOTS; output = false
-using ControlSystemsBase
-P = tf([1.], [1., 1])
-
-ζ = 0.5 # Desired damping
-ws = exp10.(range(-1, stop=2, length=8)) # A vector of closed-loop bandwidths
-kp = 2*ζ*ws .- 1 # Simple pole placement with PI given the closed-loop bandwidth, the poles are placed in a butterworth pattern
-ki = ws.^2
-
-ω = exp10.(range(-3, stop = 2, length = 500))
-pidplots(
-    P,
-    :nyquist;
-    params_p = kp,
-    params_i = ki,
-    ω = ω,
-    ylims = (-2, 2),
-    xlims = (-3, 3),
-    form = :parallel,
-)
-save_docs_plot("pidplotsnyquist1.svg") # hide
-pidplots(P, :gof; params_p = kp, params_i = ki, ω = ω, legend = false, form=:parallel, legendfontsize=6, size=(1000, 1000))
-# You can also request both Nyquist and Gang-of-four plots (more plots are available, see ?pidplots ):
-# pidplots(P,:nyquist,:gof;kps=kp,kis=ki,ω=ω);
-save_docs_plot("pidplotsgof1.svg"); # hide
-
-# output
-
-```
-![](../../plots/pidplotsnyquist1.svg)
-![](../../plots/pidplotsgof1.svg)
-
-
-Now try a different strategy, where we have specified a gain crossover frequency of 0.1 rad/s
-```jldoctest PIDPLOTS; output = false
-kp = range(-1, stop=1, length=8) #
-ki = sqrt.(1 .- kp.^2)/10
-
-pidplots(P,:nyquist,;params_p=kp,params_i=ki,ylims=(-1,1),xlims=(-1.5,1.5), form=:parallel)
-save_docs_plot("pidplotsnyquist2.svg") # hide
-pidplots(P,:gof,;params_p=kp,params_i=ki,legend=false,ylims=(0.08,8),xlims=(0.003,20), form=:parallel, legendfontsize=6, size=(1000, 1000))
-save_docs_plot("pidplotsgof2.svg"); # hide
-
-# output
-
-```
-![](../../plots/pidplotsnyquist2.svg)
-![](../../plots/pidplotsgof2.svg)
 
 ## Further examples
 - See the [examples folder](https://github.com/JuliaControl/ControlSystems.jl/tree/master/example) as well as the notebooks in [ControlExamples.jl](https://github.com/JuliaControl/ControlExamples.jl).
diff --git a/lib/ControlSystemsBase/src/pid_design.jl b/lib/ControlSystemsBase/src/pid_design.jl
index 1e223a73e..9e3ff7ef4 100644
--- a/lib/ControlSystemsBase/src/pid_design.jl
+++ b/lib/ControlSystemsBase/src/pid_design.jl
@@ -227,6 +227,40 @@ and should be supplied as additional arguments to the function.
 One can also supply a frequency vector `ω` to be used in Bode and Nyquist plots.
 
 See also `loopshapingPI`, `stabregionPID`
+
+## Example
+This example utilizes the function [`pidplots`](@ref), which accepts vectors of PID-parameters and produces relevant plots. The task is to take a system with bandwidth 1 rad/s and produce a closed-loop system with bandwidth 0.1 rad/s. If one is not careful and proceed with pole placement, one easily get a system with very poor robustness.
+```julia
+using ControlSystemsBase, Plots
+P = tf([1.], [1., 1])
+
+ζ = 0.5 # Desired damping
+ws = exp10.(range(-1, stop=2, length=8)) # A vector of closed-loop bandwidths
+kp = 2*ζ*ws .- 1 # Simple pole placement with PI given the closed-loop bandwidth, the poles are placed in a butterworth pattern
+ki = ws.^2
+
+ω = exp10.(range(-3, stop = 2, length = 500))
+pidplots(
+    P,
+    :nyquist;
+    params_p = kp,
+    params_i = ki,
+    ω = ω,
+    ylims = (-2, 2),
+    xlims = (-3, 3),
+    form = :parallel,
+)
+pidplots(P, :gof; params_p = kp, params_i = ki, ω = ω, legend = false, form=:parallel, legendfontsize=6, size=(1000, 1000))
+```
+Now try a different strategy, where we have specified a gain crossover frequency of 0.1 rad/s
+```julia
+kp = range(-1, stop=1, length=8) #
+ki = sqrt.(1 .- kp.^2)/10
+
+pidplots(P,:nyquist,;params_p=kp,params_i=ki,ylims=(-1,1),xlims=(-1.5,1.5), form=:parallel)
+pidplots(P,:gof,;params_p=kp,params_i=ki,legend=false,ylims=(0.08,8),xlims=(0.003,20), form=:parallel, legendfontsize=6, size=(1000, 1000))
+```
+
 """
 function pidplots(P::LTISystem, args...; 
     params_p, params_i, params_d=0, 
@@ -375,6 +409,20 @@ arg(P) + arg(C) = -π
 ```
 If `P` is a function (e.g. s -> exp(-sqrt(s)) ), the stability of feedback loops using PI-controllers can be analyzed for processes with models with arbitrary analytic functions
 See also [`loopshapingPI`](@ref), [`loopshapingPID`](@ref), [`pidplots`](@ref)
+
+## Example
+The stability boundary, i.e., the surface of PID parameters where the transfer function ``P(s)C(s)`` equals -1, can be plotted with the command [`stabregionPID`](@ref). The process can be given in function form or as a regular LTIsystem.
+
+```julia
+P1 = s -> exp(-sqrt(s))
+doplot = true
+form = :parallel
+kp, ki, f1 = stabregionPID(P1,exp10.(range(-5, stop=1, length=1000)); doplot, form); f1
+P2 = s -> 100*(s+6).^2. /(s.*(s+1).^2. *(s+50).^2)
+kp, ki, f2 = stabregionPID(P2,exp10.(range(-5, stop=2, length=1000)); doplot, form); f2
+P3 = tf(1,[1,1])^4
+kp, ki, f3 = stabregionPID(P3,exp10.(range(-5, stop=0, length=1000)); doplot, form); f3
+```
 """
 function stabregionPID(P, ω = _default_freq_vector(P,Val{:bode}()); kd=0, form=:standard, doplot=false)
     Pv  = freqrespv(P,ω)