diff --git a/.gitignore b/.gitignore
index 3f0cc4e6e9a0002812891f6f95c867bdbe76f272..82c4583d8c90e9ab70dec1fe990abb955bfd8ad6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,4 @@
 /src/IQCs.ipynb
 /src/NeuralNetworks.ipynb
 /examples/random_systems.ipynb
+/examples/ToolboxExamples.ipynb
diff --git a/examples/ToolboxExamples.jl b/examples/ToolboxExamples.jl
new file mode 100644
index 0000000000000000000000000000000000000000..da87fd5641e07bd4a3541bf3f8129a0d016d14a0
--- /dev/null
+++ b/examples/ToolboxExamples.jl
@@ -0,0 +1,732 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .jl
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.3
+#   kernelspec:
+#     display_name: Julia 1.10.2
+#     language: julia
+#     name: julia-1.10
+# ---
+
+# %% [markdown]
+# # Prelude
+
+# %%
+import Pkg; Pkg.activate("."); Pkg.develop(path=".."); Pkg.resolve(); Pkg.instantiate(); 
+
+# %%
+using Revise
+
+# %%
+using IQC
+using Plots
+using ProgressMeter
+
+# %%
+import MosekTools
+
+# %%
+import ControlSystems
+CS = ControlSystems
+s = CS.tf("s")
+
+# %% [markdown]
+# **TODO:** Add simulations for all cases
+
+# %% [markdown]
+# **TODO:** Perform everything in discrete time as well. Check if disturbance-free convergece rate and robustness are different things there. In continuous time, I expect them to be the same and be high-gain control.
+
+# %% [markdown]
+# **TODO:** Run longer trainings on both convergence rate and gain bound, plot controller parameter values versus gradient steps, check if path followed is similar 
+
+# %% [markdown]
+# **TODO:** Plot both gain bound and convergence rate versus gradient steps in all training examples. Use missing for gain bound if optimization fails. This is also intended to check if they are related.
+
+# %% [markdown]
+# **TODO:** Add checks using other methods for P and PID control cases
+
+# %% [markdown]
+# **TODO:** Introduce uncertainty as an unknown parameter in model estimation. As a simple example, take $x = 1/(s+\beta)$ and pretend that I only know $\beta \in [0.1, 0.3]$. Rewrite as $sx+\beta x = u$ and then $sx = -(\beta_{\mbox{nominal}} + \delta \beta ) x + u$ with the new signal $(\delta \beta \; x)^2 \le 0.1^2 x^2$. Do not use the uncertainty in $u$ that I had before.
+
+# %% [markdown]
+# **TODO:** Ha en gain bound från mätbrus till tillstånd. Gör detta för alla processer.
+
+# %% [markdown]
+# # P-controller
+
+# %% [markdown]
+# ## Continuous Times
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    w, 1
+end
+
+@parameters begin
+    k = 0
+end
+
+@variables sys begin
+    γ², Scalar(), Min
+end
+
+@assume sys begin
+    x = (1/s)*(u+w)
+    u = -k*x
+end
+
+@require sys x^2 <= γ²*w^2
+
+# %%
+for k_val in [ 0.0, 1.0 ]
+    set!(k, k_val)
+    if solve!(sys, silent=true)
+        gain_bound = round(sqrt(IQC.value(γ²)), sigdigits = 2)
+        println("Succesfully verified a gain bound γ=$(gain_bound) for k=$(k_val)")
+    else
+        println("Failed to verify a finite gain bound for k=$(k_val)")
+    end
+end
+
+# %%
+k_vals = 0:.01:1
+γ²_vals = []
+@showprogress for k_val in k_vals
+    set!(k, k_val)
+    if solve!(sys, silent=true)
+        γ²_val = IQC.value(γ²)
+    else
+        γ²_val = missing    
+    end
+    push!(γ²_vals, γ²_val)
+end
+
+# %%
+plot(k_vals, sqrt.(γ²_vals), xlabel = "Control Gain", ylabel = "Gain Bound", yscale = :log10, label = "")
+vline!([0.0], color = :red, label = "")
+
+# %% [markdown]
+# # P controller with uncertainty
+
+# %% [markdown]
+# **TODO:** It would be good if any scalar inequality assumption not containing a variable gets a nonnegative variable.
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    δu, 1
+    w, 1
+end
+
+@parameters begin
+    k = 0
+    δk² = 0.1^2
+end
+
+@variables sys begin
+    τ >= 0.0
+    γ², Scalar(), Min
+end
+
+@assume sys begin
+    x = 1/s*(u+w)
+    u = -k*x + δu
+    τ*δu^2 <= τ*δk²*x^2
+end
+
+@require sys x^2 <= γ²*w^2
+
+# %%
+for k_val in [ 0.0, 1.0 ]
+    set!(k, k_val)
+    if solve!(sys, silent=true)
+        gain_bound = round(sqrt(IQC.value(γ²)), sigdigits = 2)
+        println("Succesfully verified a gain bound γ=$(gain_bound) for k=$(k_val)")
+    else
+        println("Failed to verify a finite gain bound for k=$(k_val)")
+    end
+end
+
+# %%
+k_vals = 0:.001:1
+γ²_vals = []
+@showprogress for k_val in k_vals
+    set!(k, k_val)
+    if solve!(sys, silent=true)
+        γ²_val = IQC.value(γ²)
+    else
+        γ²_val = missing    
+    end
+    push!(γ²_vals, γ²_val)
+end
+
+# %%
+plot(k_vals, sqrt.(γ²_vals), xlabel = "Control Gain", ylabel = "Gain Bound", yscale = :log10, label = "")
+vline!([0.1], color = :red, label = "0.1")
+
+# %% [markdown]
+# **TODO:** Add an NN controller
+#
+# If I can use flux with parameters, then that would be great.
+
+# %% [markdown]
+# # PID
+
+# %% [markdown]
+# PID controller of a high order system. Tore says that having a few zeros makes it a lot trickier, so I should have stable poles and som eright-half-plane zeros. Find a good gain bound using gradients.
+#
+# Placing poles at $-1$, $-1/2\pm i$, $1/2\pm i$ and zeros at $1$ and $2$
+#
+# Take small steps
+
+# %% [markdown]
+# ## Vanilla case
+
+# %%
+k_val = zeros(3)
+
+G = 1/(s-1)^2
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+function CS_gain_bound(k_val)
+    G = 1/(s-1)^2
+    C = k_val'*vcat(1, 1/s, s)
+    
+    sys = CS.feedback(G, C)
+    (g, p) = CS.bode(sys)
+    maximum(g)
+end
+
+function CS_conv_rate(k_val)
+    G = 1/(s-1)^2
+    C = k_val'*vcat(1, 1/s, s)
+    
+    sys = CS.feedback(G, C)
+    -maximum(real.(CS.pole(sys)))
+end
+
+# %%
+CS_conv_rate([1,.1,.1])
+
+# %%
+CS_gain_bound([1, .1, .1])
+
+# %% [markdown]
+# **TODO:** I should plot the P, I and D, in  both cases, probably.
+
+# %%
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1], method = :zoh)
+plot()
+for ix in 1:2
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %% [markdown]
+# **TODO:** Add a macro with a nice syntax for delayed equalities. (Can I detect automatically?)
+
+# %%
+G_pid = vcat(1, 1/s, s)
+
+# %% [markdown]
+# **TODO:** Log PID params and comment on them tracing a curve if they do. Do this in all cases.
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    x_pid, 3
+end
+
+@parameters begin
+    k = zeros(3, 1)
+end
+
+@assume sys begin
+    x = G*u
+    x_pid = G_pid*x
+end
+
+IQC.add_delayed_linear_equality!(sys, u, -k'*x_pid)
+
+y = x_pid
+K = [k]
+
+conv_rate_sys = (; sys, y, k )
+
+# %%
+IQC.convergence_rate(sys, y)
+
+# %%
+IQC.convergence_rate_gradient(sys, y, K)
+
+# %%
+rates = []
+ks = []
+λ = 1e-1
+
+@showprogress for _ in 1:100
+    (; conv_rate, grads ) = IQC.convergence_rate_gradient(sys, y, K)
+    push!(rates, conv_rate)
+    push!(ks, copy(value(k)))
+    set!(k, value(k) + λ * grads[k])
+end
+
+IQC.convergence_rate(sys, y)
+
+# %%
+plot(rates, label = "")
+hline!([0.0], label = "", c = :red)
+plot!(title = "Convergence Rates", xaxis = "Gradient Steps", yaxis = "Convergence Rate")|> display
+
+# %%
+rates = CS_conv_rate.(ks)
+plot(rates, label = "")
+plot!(title = "Convergence Rates from ControlSystems.jl", xaxis = "Gradient Step", yaxis = "Convergence Rate")
+
+# %%
+plt = plot()
+
+for ix in eachindex(ks[1])
+    plot!(getindex.(ks, ix), label = "k[$ix]")
+end
+plot!()
+
+# %%
+k_val = value(k)
+
+G = 1/(s-1)^2
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+# %%
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1], method = :zoh)
+plot()
+for ix in 1:3
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %%
+plt = CS.bodeplot(sys, label = "Closed Loop")
+hline!(plt[1], [1.0], label = "1.0")
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    w, 1
+    x_pid, 3
+end
+
+@parameters begin
+    k = value(k)
+end
+
+@assume sys begin
+    x = G*(u+w)
+    x_pid = G_pid*x
+end
+
+IQC.add_delayed_linear_equality!(sys, u, -k'*x_pid)
+
+@variables sys begin
+    γ², Scalar(), Min
+end
+
+@require sys x_pid^2 <= γ²*w^2
+
+y = x_pid
+K = [k]
+
+# %%
+gain_bounds = []
+for k_val in ks
+    set!(k, k_val)
+    if solve!(sys, silent = true)
+        push!(gain_bounds, sqrt(value(γ²)))
+    else
+        push!(gain_bounds, missing)
+    end
+end
+
+plot(gain_bounds, yscale = :log10)
+
+# %% [markdown]
+# **TODO:** This plot shows gain bound from w to x but I think the above plot shows gain bound to x_pid. I should build yet another system to obtain gain bounds to x only.
+
+# %%
+gain_bounds = []
+for k_val in ks
+    push!(gain_bounds, CS_gain_bound(k_val))
+end
+plot(gain_bounds)
+plot!(yaxis = :log10)
+    
+
+# %%
+solve!(sys)
+
+# %%
+value(γ²)
+
+# %%
+value(k)
+
+# %%
+sys′ = IQC.apply_KYP(sys)
+IQC.gradient(sys′, K)
+
+# %%
+ks = []
+γ²s = []
+λ = 1e-2
+
+@showprogress for _ in 1:100
+    solve!(sys′, silent = true)
+    push!(γ²s, value(γ²))
+    push!(ks, copy(value(k)))
+    grads = IQC.gradient(sys′, K)
+    set!(k, value(k) - λ * grads[k])
+end
+
+plot(sqrt.(γ²s), xaxis = "Gradient Steps", yscale = :log10, yaxis = "Gain Bound", label = "") |> display
+
+solve!(sys′)
+value(γ²)
+
+# %% [markdown]
+# **TODO:** Use IQC.jl to get gain bound from w to x only and compare to that. Can I get gain bound to all as max of gain bounds? If gain bounds are correlated for a specific w, then maximal gain bound of components might be lower than true gain bound for MIMO system. Bodemag can get the gain bound to the euclidian norm of the magnitude, and it has an implementatin in `ControlSystems.jl`.
+
+# %%
+plot(CS_gain_bound.(ks))
+
+# %%
+plot()
+for ix in eachindex(ks[begin])
+    plot!(getindex.(ks, ix), label = "k[ix]")
+end
+plot!()
+
+# %%
+rates = []
+
+for k_val in ks
+    set!(conv_rate_sys.k, k_val)
+    push!(rates, IQC.convergence_rate(conv_rate_sys.sys, conv_rate_sys.y))
+end
+
+plot(rates, label = "")
+
+# %%
+plot(CS_conv_rate.(ks))
+
+# %%
+value(k)
+
+# %%
+k_val = value(k)
+
+G = 1/(s-1)^2
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+# %%
+(gain, phase) = CS.bode(sys)
+gain_bound = maximum(gain)
+
+# %%
+plt = CS.bodeplot(sys, label = "Closed Loop")
+hline!(plt[1], [gain_bound], color = :green, label = "$(round(gain_bound, sigdigits = 2))")
+hline!(plt[1], [1.0], color = :red, label = "1.0")
+
+# %%
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1], method = :zoh)
+plot()
+for ix in 1:3
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %% [markdown]
+# ## Adding Zeros
+
+# %% [markdown]
+# **TODO:** Plot both convergence rate and gain_bound in all cases. Create systems for both for all cases, and check if convergence rate improved gain bound.
+
+# %% [markdown]
+# The state is stable, but the integrator signal is not. Talk to Tore about what I should want.
+
+# %%
+k_val = zeros(3)
+G = (s-1)/((s+3)*(s+2)*(s+1))
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1], method = :zoh)
+plot()
+for ix in 1:3
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %%
+k_val = 1e-3*ones(3)
+G = ((s-1)/((s+3)*(s+2)*(s+1)))
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1, 1], method = :zoh)
+plot()
+for ix in 1:4
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    x_pid, 3
+end
+
+@parameters begin
+    k = zeros(3, 1)
+end
+
+G_pid = vcat(1, 1/s, s)
+@assume sys begin
+    x = ((s-1)/((s+2)*(s+1)*(s+3)))*u
+    x_pid = G_pid*x
+end
+
+IQC.add_delayed_linear_equality!(sys, u, -k'*x_pid)
+
+y = x_pid
+K = [k]
+
+@require sys x_pid^2 <= 0
+
+# %%
+(; x, u) = IQC.SimplifyAndSolve.get_KYP_signals(sys)
+
+# %%
+sys.lin_eqs[y]
+
+# %%
+IQC.convergence_rate(sys, x, margin=1e-5)
+
+# %%
+IQC.convergence_rate_gradient(sys, x, K)
+
+# %%
+rates = []
+λ = 1e-1
+
+@showprogress for _ in 1:100
+    (; conv_rate, grads ) = IQC.convergence_rate_gradient(sys, x, K)
+    push!(rates, conv_rate)
+    set!(k, value(k) + λ * grads[k])
+end
+
+plot(rates, xaxis = "Gradient Steps", yaxis = "Convergence Rate", label = "")
+hline!([0.0], label = "", c = :red) |> display
+
+IQC.convergence_rate(sys, x)
+
+# %%
+k_val = value(k)
+G = ((s-1)/((s+3)*(s+2)*(s+1)))
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1, 1], method = :zoh)
+plot()
+for ix in 1:4
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    w, 1
+    x_pid, 3
+end
+
+@parameters begin
+    k = value(k)
+end
+
+G_pid = vcat(1, 1/s, s)
+@assume sys begin
+    x = ((s-1)/((s+1)*(s+2)*(s+3)))*(u+w)
+    x_pid = G_pid*x
+end
+
+IQC.add_delayed_linear_equality!(sys, u, -k'*x_pid)
+
+@variables sys begin
+    γ², Scalar(), Min
+end
+
+@require sys x_pid^2 <= γ²*w^2
+
+y = x_pid
+K = [k]
+
+# %%
+solve!(sys)
+
+# %%
+value(γ²)
+
+# %%
+value(k)
+
+# %%
+sys′ = IQC.apply_KYP(sys)
+IQC.gradient(sys′, K)
+
+# %%
+γ²s = []
+λ = 1e-2
+
+@showprogress for _ in 1:1000
+    solve!(sys′, silent = true)
+    push!(γ²s, value(γ²))
+    grads = IQC.gradient(sys′, K)
+    set!(k, value(k) - λ * grads[k])
+end
+
+plot(γ²s, xaxis = "Gradient Steps", yscale = :log10, yaxis = "Gain Bound", label = "") |> display
+
+solve!(sys′)
+value(γ²)
+
+# %%
+value(k)
+
+# %%
+k_val = value(k)
+G = ((s-1)/((s+3)*(s+2)*(s+1)))
+C = k_val'*vcat(1, 1/s, s)
+
+sys = CS.feedback(G, C)
+
+ts = 0:.01:10
+y, t, x, u = CS.lsim(sys, zeros(1, length(ts)), ts; x0 = [1, 1, 1, 1], method = :zoh)
+plot()
+for ix in 1:4
+    plot!(t, x[ix, :], xaxis = "t", yaxis = "State", label = "x[$ix]")
+end
+plot!()
+
+# %% [markdown]
+# ## More PID?
+#
+# Consider adding oscillations, making the base system unstable, and adding uncertainty.
+
+# %% [markdown]
+# # NN Controller
+
+# %% [markdown]
+# I really want to have bias terms in my NN, and I do not know what my state of that is currently.
+
+# %%
+n_relus = [ 8, 8 ]
+
+# %%
+sys = ContinuousTimeSystem(opt = MosekTools.Optimizer)
+
+@signals sys begin
+    x, 1
+    u, 1
+    δu, 1
+    w, 1
+    η_1, 8
+    η_2, 8
+    η_3, 1
+    ξ_1, 8
+    ξ_2, 8
+    ξ_3, 1
+end
+
+@parameters begin
+    k = 0
+    δk² = 0.1^2
+    W1 = randn(8, 1)
+    b1 = randn(8,1)
+    W2 = randn(8, 8)
+    b2 = randn(8,1)
+    W3 = randn(1, 8)
+    b3 = randn(1,1)
+end
+
+@variables sys begin
+    τ >= 0.0
+    γ², Scalar(), Min
+end
+
+unit = sys.one
+
+@assume sys begin
+    x = 1/s*(u+w)
+    τ*δu^2 <= τ*δk²*x^2
+    η_1 = W1*x+b1*unit
+    η_2 = W2*ξ_1+b2*unit
+    η_3 = W3*ξ_2+b3*unit
+    u = ξ_3 + δu
+end
+
+@IQC sys Absolute relu [η_1, η_2, η_3] [ξ_1, ξ_2, ξ_3] begin
+    QuadEq
+    CombLin
+    RateLimit
+    ZamesFalb
+end
+
+@require sys x^2 <= γ²*w^2
+
+# %%
+solve!(sys)
+
+# %%
+sys′ = IQC.apply_KYP(sys)
+IQC.gradient(sys′, W1)
+
+# %%