diff --git a/matlab/figs/bode_plot_controller.pdf b/matlab/figs/bode_plot_controller.pdf
index a91ba4e..a0606ba 100644
Binary files a/matlab/figs/bode_plot_controller.pdf and b/matlab/figs/bode_plot_controller.pdf differ
diff --git a/matlab/figs/bode_plot_controller.png b/matlab/figs/bode_plot_controller.png
index e774d97..65848b1 100644
Binary files a/matlab/figs/bode_plot_controller.png and b/matlab/figs/bode_plot_controller.png differ
diff --git a/matlab/figs/bode_plot_loop_gain.pdf b/matlab/figs/bode_plot_loop_gain.pdf
index 6f64dbf..bf5a4c7 100644
Binary files a/matlab/figs/bode_plot_loop_gain.pdf and b/matlab/figs/bode_plot_loop_gain.pdf differ
diff --git a/matlab/figs/bode_plot_loop_gain.png b/matlab/figs/bode_plot_loop_gain.png
index b92b960..36dd8d2 100644
Binary files a/matlab/figs/bode_plot_loop_gain.png and b/matlab/figs/bode_plot_loop_gain.png differ
diff --git a/matlab/figs/bode_plot_mech_sys.pdf b/matlab/figs/bode_plot_mech_sys.pdf
index 1e9a1ea..c926a52 100644
Binary files a/matlab/figs/bode_plot_mech_sys.pdf and b/matlab/figs/bode_plot_mech_sys.pdf differ
diff --git a/matlab/figs/bode_plot_mech_sys.png b/matlab/figs/bode_plot_mech_sys.png
index fb3ca67..e73bb3f 100644
Binary files a/matlab/figs/bode_plot_mech_sys.png and b/matlab/figs/bode_plot_mech_sys.png differ
diff --git a/matlab/figs/bode_requirements.pdf b/matlab/figs/bode_requirements.pdf
index 1de007b..98bd78b 100644
Binary files a/matlab/figs/bode_requirements.pdf and b/matlab/figs/bode_requirements.pdf differ
diff --git a/matlab/figs/bode_requirements.png b/matlab/figs/bode_requirements.png
index e45de81..a2f153a 100644
Binary files a/matlab/figs/bode_requirements.png and b/matlab/figs/bode_requirements.png differ
diff --git a/matlab/figs/bode_wi.pdf b/matlab/figs/bode_wi.pdf
index e499caa..8f2b74e 100644
Binary files a/matlab/figs/bode_wi.pdf and b/matlab/figs/bode_wi.pdf differ
diff --git a/matlab/figs/bode_wi.png b/matlab/figs/bode_wi.png
index 5b56254..95bfcf9 100644
Binary files a/matlab/figs/bode_wi.png and b/matlab/figs/bode_wi.png differ
diff --git a/matlab/figs/comp_filters_param_alpha.pdf b/matlab/figs/comp_filters_param_alpha.pdf
index 6c5dc4c..ba4a91f 100644
Binary files a/matlab/figs/comp_filters_param_alpha.pdf and b/matlab/figs/comp_filters_param_alpha.pdf differ
diff --git a/matlab/figs/comp_filters_param_alpha.png b/matlab/figs/comp_filters_param_alpha.png
index b637f4e..2081d1c 100644
Binary files a/matlab/figs/comp_filters_param_alpha.png and b/matlab/figs/comp_filters_param_alpha.png differ
diff --git a/matlab/figs/comp_hinf_analytical.pdf b/matlab/figs/comp_hinf_analytical.pdf
index 24f8070..d8cdaea 100644
Binary files a/matlab/figs/comp_hinf_analytical.pdf and b/matlab/figs/comp_hinf_analytical.pdf differ
diff --git a/matlab/figs/comp_hinf_analytical.png b/matlab/figs/comp_hinf_analytical.png
index 5ca712a..bbf0f64 100644
Binary files a/matlab/figs/comp_hinf_analytical.png and b/matlab/figs/comp_hinf_analytical.png differ
diff --git a/matlab/figs/compare_weights_upper_bounds_S_T.pdf b/matlab/figs/compare_weights_upper_bounds_S_T.pdf
index dfe10db..c6ec136 100644
Binary files a/matlab/figs/compare_weights_upper_bounds_S_T.pdf and b/matlab/figs/compare_weights_upper_bounds_S_T.pdf differ
diff --git a/matlab/figs/compare_weights_upper_bounds_S_T.png b/matlab/figs/compare_weights_upper_bounds_S_T.png
index 0305dc7..988b82c 100644
Binary files a/matlab/figs/compare_weights_upper_bounds_S_T.png and b/matlab/figs/compare_weights_upper_bounds_S_T.png differ
diff --git a/matlab/figs/complementary_filters_second_order.pdf b/matlab/figs/complementary_filters_second_order.pdf
index 090a8d5..b44aedb 100644
Binary files a/matlab/figs/complementary_filters_second_order.pdf and b/matlab/figs/complementary_filters_second_order.pdf differ
diff --git a/matlab/figs/complementary_filters_second_order.png b/matlab/figs/complementary_filters_second_order.png
index fa71777..b0700ff 100644
Binary files a/matlab/figs/complementary_filters_second_order.png and b/matlab/figs/complementary_filters_second_order.png differ
diff --git a/matlab/figs/complementary_filters_third_order.pdf b/matlab/figs/complementary_filters_third_order.pdf
index 1401431..25e8d87 100644
Binary files a/matlab/figs/complementary_filters_third_order.pdf and b/matlab/figs/complementary_filters_third_order.pdf differ
diff --git a/matlab/figs/complementary_filters_third_order.png b/matlab/figs/complementary_filters_third_order.png
index 12c867a..097dd08 100644
Binary files a/matlab/figs/complementary_filters_third_order.png and b/matlab/figs/complementary_filters_third_order.png differ
diff --git a/matlab/figs/hinf_filters_results.pdf b/matlab/figs/hinf_filters_results.pdf
index c76e67f..253f220 100644
Binary files a/matlab/figs/hinf_filters_results.pdf and b/matlab/figs/hinf_filters_results.pdf differ
diff --git a/matlab/figs/hinf_filters_results.png b/matlab/figs/hinf_filters_results.png
index efd7d47..2a7c2f1 100644
Binary files a/matlab/figs/hinf_filters_results.png and b/matlab/figs/hinf_filters_results.png differ
diff --git a/matlab/figs/loop_gain_robustness.pdf b/matlab/figs/loop_gain_robustness.pdf
index a10325f..87cc8f9 100644
Binary files a/matlab/figs/loop_gain_robustness.pdf and b/matlab/figs/loop_gain_robustness.pdf differ
diff --git a/matlab/figs/loop_gain_robustness.png b/matlab/figs/loop_gain_robustness.png
index d9ce8d7..f6c53fb 100644
Binary files a/matlab/figs/loop_gain_robustness.png and b/matlab/figs/loop_gain_robustness.png differ
diff --git a/matlab/figs/nyquist_robustness.pdf b/matlab/figs/nyquist_robustness.pdf
index ca3dec5..0f16d05 100644
Binary files a/matlab/figs/nyquist_robustness.pdf and b/matlab/figs/nyquist_robustness.pdf differ
diff --git a/matlab/figs/nyquist_robustness.png b/matlab/figs/nyquist_robustness.png
index b36bd7d..9f77149 100644
Binary files a/matlab/figs/nyquist_robustness.png and b/matlab/figs/nyquist_robustness.png differ
diff --git a/matlab/figs/plant_uncertainty_bode_plot.pdf b/matlab/figs/plant_uncertainty_bode_plot.pdf
index fe78089..2722d8f 100644
Binary files a/matlab/figs/plant_uncertainty_bode_plot.pdf and b/matlab/figs/plant_uncertainty_bode_plot.pdf differ
diff --git a/matlab/figs/plant_uncertainty_bode_plot.png b/matlab/figs/plant_uncertainty_bode_plot.png
index cf55d12..e800ea7 100644
Binary files a/matlab/figs/plant_uncertainty_bode_plot.png and b/matlab/figs/plant_uncertainty_bode_plot.png differ
diff --git a/matlab/figs/robust_performance_result.pdf b/matlab/figs/robust_performance_result.pdf
index 3ecbba8..cd7b60a 100644
Binary files a/matlab/figs/robust_performance_result.pdf and b/matlab/figs/robust_performance_result.pdf differ
diff --git a/matlab/figs/robust_performance_result.png b/matlab/figs/robust_performance_result.png
index 8f78d4f..f0455b9 100644
Binary files a/matlab/figs/robust_performance_result.png and b/matlab/figs/robust_performance_result.png differ
diff --git a/matlab/figs/robust_stability.pdf b/matlab/figs/robust_stability.pdf
index 2780f79..3693be0 100644
Binary files a/matlab/figs/robust_stability.pdf and b/matlab/figs/robust_stability.pdf differ
diff --git a/matlab/figs/robust_stability.png b/matlab/figs/robust_stability.png
index 340dde5..47c4d91 100644
Binary files a/matlab/figs/robust_stability.png and b/matlab/figs/robust_stability.png differ
diff --git a/matlab/figs/u_and_y_with_Kr.pdf b/matlab/figs/u_and_y_with_Kr.pdf
index 05fa3bf..3ddd69e 100644
Binary files a/matlab/figs/u_and_y_with_Kr.pdf and b/matlab/figs/u_and_y_with_Kr.pdf differ
diff --git a/matlab/figs/u_and_y_with_Kr.png b/matlab/figs/u_and_y_with_Kr.png
index 5cc4965..354cbb5 100644
Binary files a/matlab/figs/u_and_y_with_Kr.png and b/matlab/figs/u_and_y_with_Kr.png differ
diff --git a/matlab/figs/verification_NP.pdf b/matlab/figs/verification_NP.pdf
index 4d05c3c..172a2bc 100644
Binary files a/matlab/figs/verification_NP.pdf and b/matlab/figs/verification_NP.pdf differ
diff --git a/matlab/figs/verification_NP.png b/matlab/figs/verification_NP.png
index b1ffef6..b98d56b 100644
Binary files a/matlab/figs/verification_NP.png and b/matlab/figs/verification_NP.png differ
diff --git a/matlab/figs/weights_NP_RS_RP.pdf b/matlab/figs/weights_NP_RS_RP.pdf
index ac84737..6208d6d 100644
Binary files a/matlab/figs/weights_NP_RS_RP.pdf and b/matlab/figs/weights_NP_RS_RP.pdf differ
diff --git a/matlab/figs/weights_NP_RS_RP.png b/matlab/figs/weights_NP_RS_RP.png
index 2cb7bac..65d7f5b 100644
Binary files a/matlab/figs/weights_NP_RS_RP.png and b/matlab/figs/weights_NP_RS_RP.png differ
diff --git a/matlab/figs/weights_wl_wh.pdf b/matlab/figs/weights_wl_wh.pdf
index c1c4146..c40bc1e 100644
Binary files a/matlab/figs/weights_wl_wh.pdf and b/matlab/figs/weights_wl_wh.pdf differ
diff --git a/matlab/figs/weights_wl_wh.png b/matlab/figs/weights_wl_wh.png
index d134aa2..66e70fc 100644
Binary files a/matlab/figs/weights_wl_wh.png and b/matlab/figs/weights_wl_wh.png differ
diff --git a/matlab/index.org b/matlab/index.org
index c71fd75..c0d5e2b 100644
--- a/matlab/index.org
+++ b/matlab/index.org
@@ -33,19 +33,32 @@ The control architecture studied here is shown on figure [[fig:sf_arch_class_pre
 
 #+name: fig:sf_arch_class_prefilter
 #+caption: Control Architecture
-[[file:figs-tikz/sf_arch_class_prefilter.png]]
+[[file:figs-paper/sf_arch_class_prefilter.png]]
 
 Here is the outline of the =matlab= analysis for this control architecture:
-- Section [[sec:plant]]: the plant model $G$ is defined
-- Section [[sec:uncertainty]]: the plant uncertainty set $\Pi_I$ is defined using the multiplicative input uncertainty: $\Pi_I: \ G^\prime = G (1 + w_I \Delta)$. Thus the weight $w_I$ is defined such that the true system dynamics is included in the set $\Pi_I$
-- Section [[sec:specifications]]: From the specifications on performance that are expressed in terms of upper bounds of $S$ and $T$, performance weights $w_S$ and $w_T$ are derived such that the goal is to obtain $|S| < \frac{1}{|w_S|}$ and $|T| < \frac{1}{|w_T|}, \ \forall \omega$
-- Section [[sec:upper_bounds_filters]]: upper bounds on the magnitude of the complementary filters $|H_L|$ and $|H_H|$ are defined in order to ensure Nominal Performance (NP), Robust Stability (RS) and Robust Performance (RP)
-- Then, $H_L$ and $H_H$ are synthesize such that $|H_L|$ and $|H_H|$ are within the specified bounds and such that $H_L + H_H = 1$ (complementary property). This is done using two techniques, first $\mathcal{H}_\infty$ (section [[sec:h_infinity]]) and then analytical formulas (section [[sec:analytical_formula]]). Resulting complementary filters for both methods are compared in section [[sec:comp_filters]].
-- Section [[sec:controller_analysis]]: the obtain controller $K = G^{-1} H_H^{-1}$ is analyzed
-- Section [[sec:nominal_stability_performance]]: the Nominal Stability (NS) and Nominal Performance conditions are verified
-- Section [[sec:robustness_analysis]]: robust Stability and Robust Performance conditions are studied
-- Section [[sec:pre_filter]]: a pre-filter that is used to limit the input usage due to the change of the reference is added
-- Section [[sec:sisotool_controller]]: a controller is designed using =SISOTOOL= and then compared with the previously generated controller
+- Section [[sec:plant]]:
+  The plant model $G$ is defined
+- Section [[sec:uncertainty]]:
+  The plant uncertainty set $\Pi_I$ is defined using the multiplicative input uncertainty: $\Pi_I: \ G^\prime = G (1 + w_I \Delta)$.
+  Thus the weight $w_I$ is defined such that the true system dynamics is included in the set $\Pi_I$
+- Section [[sec:specifications]]:
+  From the specifications on performance that are expressed in terms of upper bounds of $S$ and $T$, performance weights $w_S$ and $w_T$ are derived such that the goal is to obtain $|S| < \frac{1}{|w_S|}$ and $|T| < \frac{1}{|w_T|}, \ \forall \omega$
+- Section [[sec:upper_bounds_filters]]:
+  Upper bounds on the magnitude of the complementary filters $|H_L|$ and $|H_H|$ are defined in order to ensure Nominal Performance (NP), Robust Stability (RS) and Robust Performance (RP)
+- Section [[sec:h_infinity]]:
+  $H_L$ and $H_H$ are synthesize using the $\mathcal{H}_\infty$ synthesis such that $|H_L|$ and $|H_H|$ are within the specified bounds and such that $H_L + H_H = 1$ (complementary property)
+- Section [[sec:analytical_formula]]:
+  Analytical formulas of complementary filters are used
+- Section [[sec:comp_filters]]:
+  The obtained complementary filters using the $\mathcal{H}_\infty$ synthesis and the analytical formula are compared
+- Section [[sec:controller_analysis]]:
+  The obtain controller $K = G^{-1} H_H^{-1}$ is analyzed
+- Section [[sec:nominal_stability_performance]]:
+  The Nominal Stability (NS) and Nominal Performance conditions are verified
+- Section [[sec:robustness_analysis]]:
+  Robust Stability and Robust Performance conditions are studied
+- Section [[sec:pre_filter]]:
+  A pre-filter that is used to limit the input usage due to the change of the reference is added
 
 * ZIP file containing the data and matlab files                      :ignore:
 #+begin_src bash :exports none :results none
@@ -71,86 +84,86 @@ Here is the outline of the =matlab= analysis for this control architecture:
 #+end_src
 
 #+begin_src matlab
-  freqs = logspace(-1, 3, 1000);
+  addpath('src');
+  freqs = 0.5*logspace(0, 3, 1000);
 #+end_src
 
 * Definition of the plant
-  <<sec:plant>>
+<<sec:plant>>
 
 The studied system consists of a solid positioned on top of a motorized uni-axial soft suspension.
 
 The absolute position $x$ of the solid is measured using an inertial sensor and a force $F$ can be applied to the mass using a voice coil actuator.
 
-The model of the system is represented on figure [[fig:mech_sys_alone]] where the mass of the solid is $m = 20\ [kg]$, the stiffness of the suspension is $k = 10^4\ [N/m]$ and the damping of the system is $c = 10^2\ [N/(m/s)]$.
+The model of the system is represented on figure [[fig:mech_sys_alone]] where the mass of the solid is $m = 10\ [kg]$, the stiffness of the suspension is $k = 10^4\ [N/m]$ and the damping of the system is $c = 10^2\ [N/(m/s)]$.
 
 #+name: fig:mech_sys_alone
 #+caption: One degree of freedom system
-[[file:figs-tikz/mech_sys_alone.png]]
+[[file:figs-paper/mech_sys_alone.png]]
 
 The plant $G$ is defined on matlab and its bode plot is shown on figure [[fig:bode_plot_mech_sys]].
 
 #+begin_src matlab
-  m = 20;  % [kg]
-  k = 1e4; % [N/m]
-  c = 1e2; % [N/(m/s)]
+  m = 10;  % mass [kg]
+  k = 1e4; % stiffness [N/m]
+  c = 1e2; % damping [N/(m/s)]
 
   G = 1/(m*s^2 + c*s + k);
 #+end_src
 
 #+begin_src matlab :exports none
   figure;
-  ax1 = subaxis(2,1,1);
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+  % Magnitude
+  ax1 = nexttile([2, 1]);
   hold on;
   plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
   hold off;
-  xlim([0.1, 100]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude [m/N]');
 
-  ax2 = subaxis(2,1,2);
+  % Phase
+  ax2 = nexttile;
   hold on;
   plot(freqs, 180/pi*angle(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
   hold off;
-  yticks(-180:90:180);
-  ylim([-180 180]);
-  xlim([0.1, 1000]);
+  yticks(-180:90:0);
+  ylim([-180 0]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
   linkaxes([ax1,ax2],'x');
+  xlim([freqs(1), freqs(end)])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/bode_plot_mech_sys.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/bode_plot_mech_sys.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:bode_plot_mech_sys
-#+CAPTION: Bode plot of $G$ ([[./figs/bode_plot_mech_sys.png][png]], [[./figs/bode_plot_mech_sys.pdf][pdf]])
+#+name: fig:bode_plot_mech_sys
+#+caption: Bode plot of $G$
+#+RESULTS:
 [[file:figs/bode_plot_mech_sys.png]]
 
 * Multiplicative input uncertainty
-  <<sec:uncertainty>>
+<<sec:uncertainty>>
 We choose to use the multiplicative input uncertainty to model the plant uncertainty:
-\[ \Pi_I: \ G^\prime(s) = G(s) (1 + w_I(s) \Delta(s)),\text{ with } |\Delta(j\omega)| < 1 \ \forall \omega \]
+\begin{equation}
+  \Pi_I: \ G^\prime(s) = G(s) (1 + w_I(s) \Delta(s)),\text{ with } |\Delta(j\omega)| < 1 \ \forall \omega
+\end{equation}
 
 
-The uncertainty weight $w_I$ has the following form:
-\[ w_I(s) = \frac{\tau s + r_0}{(\tau/r_\infty) s + 1} \]
-where $r_0=0.1$ is the relative uncertainty at steady-state, $1/\tau=80\text{Hz}$ is the frequency at which the relative uncertainty reaches 100%, and $r_\infty=10$ is the magnitude of the weight at high frequency.
-
-We defined the uncertainty weight on matlab. Its bode plot is shown on figure [[fig:bode_wi]].
-
+The uncertainty weight is shown in figure [[fig:bode_wi]].
 #+begin_src matlab
-  r0 = 0.1;
-  rinf = 10;
-  tau = 1/2/pi/80;
-
-  wI = (tau*s + r0)/((tau/rinf)*s+1);
+  wI = createWeight('n', 2, 'w0', 2*pi*80, 'G0', 0.1, 'G1', 10, 'Gc', 1);
 #+end_src
 
 #+begin_src matlab :exports none
   figure;
+  tiledlayout(1, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
   hold on;
   plot(freqs, abs(squeeze(freqresp(wI, freqs, 'Hz'))), 'k-');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -158,67 +171,66 @@ We defined the uncertainty weight on matlab. Its bode plot is shown on figure [[
   hold off;
   xlim([freqs(1), freqs(end)]);
   ylim([1e-1, 10]);
-  xticks([0.1, 1, 10, 100, 1000])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/bode_wi.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/bode_wi.pdf', 'width', 'half', 'height', 'small');
 #+end_src
 
-#+NAME: fig:bode_wi
-#+CAPTION: Bode plot of $w_I$ ([[./figs/bode_wi.png][png]], [[./figs/bode_wi.pdf][pdf]])
+#+name: fig:bode_wi
+#+CAPTION: Bode plot of $w_I$
+#+RESULTS:
 [[file:figs/bode_wi.png]]
 
-The uncertain model is created with the =ultidyn= function. Elements in the uncertainty set $\Pi_I$ are computed and their bode plot is shown on figure [[fig:plant_uncertainty_bode_plot]].
+Elements in the uncertainty set $\Pi_I$ are computed and their bode plot is shown on figure [[fig:plant_uncertainty_bode_plot]].
 
-#+begin_src matlab
-  Delta = ultidyn('Delta', [1 1]);
+#+begin_src matlab :exports none
+  % Delta = ultidyn('Delta', [1 1]);
 
-  Gd = G*(1+wI*Delta);
-  Gds = usample(Gd, 20);
+  % Gd = G*(1+wI*Delta);
+  % Gds = usample(Gd, 20);
 #+end_src
 
 #+begin_src matlab :exports none
   figure;
-  ax1 = subplot(2,1,1);
+  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+  % Magnitude
+  ax1 = nexttile;
   hold on;
-  for i=1:length(Gds)
-      plot(freqs, abs(squeeze(freqresp(Gds(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.1]);
-  end
-  plot(freqs, abs(squeeze(freqresp(Gd, freqs, 'Hz'))), 'k-');
+  plotMagUncertainty(wI, freqs, 'G', G, 'DisplayName', '$G$');
+  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude [m/N]');
+  ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]);
+  ylim([1e-8, 1e-3]);
   hold off;
+
   % Phase
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
-  for i=1:length(Gds)
-      plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gds(:, :, i), freqs, 'Hz')))), '-', 'color', [0, 0, 0, 0.1]);
-  end
+  plotPhaseUncertainty(wI, freqs, 'G', G);
   plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz')))), 'k-');
   set(gca,'xscale','log');
-  yticks(-360:90:180);
-  ylim([-360 0]);
+  yticks(-360:90:90);
+  ylim([-360 90]);
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
   hold off;
+
   linkaxes([ax1,ax2],'x');
   xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/plant_uncertainty_bode_plot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/plant_uncertainty_bode_plot.pdf', 'width', 'half', 'height', 'tall');
 #+end_src
 
-#+NAME: fig:plant_uncertainty_bode_plot
-#+CAPTION: Some elements in the uncertainty set $\Pi_I$ ([[./figs/plant_uncertainty_bode_plot.png][png]], [[./figs/plant_uncertainty_bode_plot.pdf][pdf]])
+#+name: fig:plant_uncertainty_bode_plot
+#+caption: Some elements in the uncertainty set $\Pi_I$
+#+RESULTS:
 [[file:figs/plant_uncertainty_bode_plot.png]]
 
 * Specifications and performance weights
-  <<sec:specifications>>
+<<sec:specifications>>
 
 The control objective is to isolate the displacement $x$ of the mass from the ground motion $w$.
 
@@ -230,6 +242,7 @@ These specifications can be represented as upper bounds on the closed loop trans
 
 #+begin_src matlab :exports none
   figure;
+
   hold on;
   set(gca,'ColorOrderIndex',1)
   plot([100, 1000], [0.1, 0.001], ':', 'DisplayName', '$|T|$ - Upper bound');
@@ -242,16 +255,16 @@ These specifications can be represented as upper bounds on the closed loop trans
   xlim([freqs(1), freqs(end)]);
   ylim([1e-3, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'northeast', 'FontSize', 8);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/bode_requirements.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/bode_requirements.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:bode_requirements
-#+CAPTION: Upper bounds on $S$ and $T$ ([[./figs/bode_requirements.png][png]], [[./figs/bode_requirements.pdf][pdf]])
+#+name: fig:bode_requirements
+#+caption: Upper bounds on $S$ and $T$
+#+RESULTS:
 [[file:figs/bode_requirements.png]]
 
 We now define two weights, $w_S(s)$ and $w_T(s)$ such that $1/|w_S|$ and $1/|w_T|$ are lower than the previously defined upper bounds.
@@ -281,20 +294,20 @@ The weights are defined as follow. They magnitude is compared with the upper bou
   xlim([freqs(1), freqs(end)]);
   ylim([1e-4, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/compare_weights_upper_bounds_S_T.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/compare_weights_upper_bounds_S_T.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:compare_weights_upper_bounds_S_T
-#+CAPTION: Weights $w_S$ and $w_T$ with the upper bounds on $S$ and $T$ obtained from the specifications ([[./figs/compare_weights_upper_bounds_S_T.png][png]], [[./figs/compare_weights_upper_bounds_S_T.pdf][pdf]])
+#+name: fig:compare_weights_upper_bounds_S_T
+#+caption: Weights $w_S$ and $w_T$ with the upper bounds on $S$ and $T$ obtained from the specifications
+#+RESULTS:
 [[file:figs/compare_weights_upper_bounds_S_T.png]]
 
 * Upper bounds on the norm of the complementary filters for NP, RS and RP
-  <<sec:upper_bounds_filters>>
+<<sec:upper_bounds_filters>>
 
 Now that we have defined $w_I$, $w_S$ and $w_T$, we can derive conditions for Nominal Performance, Robust Stability and Robust Performance ($j\omega$ is omitted here for readability):
 \begin{align*}
@@ -324,7 +337,7 @@ We plot these conditions on figure [[fig:weights_NP_RS_RP]].
   xlim([freqs(1), freqs(end)]);
   ylim([1e-3, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'northeast', 'FontSize', 8);
 #+end_src
 
 #+begin_src matlab :exports none
@@ -337,35 +350,35 @@ We plot these conditions on figure [[fig:weights_NP_RS_RP]].
   set(gca,'ColorOrderIndex',1)
   plot(freqs, 1./wT_resp, '--', 'DisplayName', 'NP - $H_L$');
   set(gca,'ColorOrderIndex',1)
-  plot(freqs, 1./(wT_resp .* (1 + wI_resp) + wI_resp), ':', 'DisplayName', 'RP for T - $H_L$');
+  plot(freqs, 1./(wT_resp .* (1 + wI_resp) + wI_resp), ':', 'DisplayName', 'RP - $H_L$ (T)');
   set(gca,'ColorOrderIndex',1)
   plot(freqs, 1./(wI_resp .* (2 + wI_resp)), '-.', 'DisplayName', 'RS - $H_L$');
   set(gca,'ColorOrderIndex',2)
   plot(freqs, 1./wS_resp, '--', 'DisplayName', 'NP - $H_H$');
   set(gca,'ColorOrderIndex',2)
-  plot(freqs, (1 + wI_resp)./(wS_resp .* (2 + wI_resp)), ':', 'DisplayName', 'RP for S - $H_H$');
+  plot(freqs, (1 + wI_resp)./(wS_resp .* (2 + wI_resp)), ':', 'DisplayName', 'RP - $H_H$ (S)');
   set(gca,'ColorOrderIndex',2)
-  plot(freqs, wT_resp.*(1 + wI_resp)./(wS_resp.*(wT_resp .* (1 + wI_resp) + wI_resp)), '-', 'DisplayName', 'RP for S - $H_H$');
+  plot(freqs, wT_resp.*(1 + wI_resp)./(wS_resp.*(wT_resp .* (1 + wI_resp) + wI_resp)), '-', 'DisplayName', 'RP - $H_H$ (S)');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   xlabel('Frequency [Hz]'); ylabel('Magnitude');
   hold off;
   xlim([freqs(1), freqs(end)]);
-  ylim([1e-3, 10]);
+  ylim([1e-4, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/weights_NP_RS_RP.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/weights_NP_RS_RP.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:weights_NP_RS_RP
-#+CAPTION: Upper bounds on the norm of the complementary filters for NP, RS and RP ([[./figs/weights_NP_RS_RP.png][png]], [[./figs/weights_NP_RS_RP.pdf][pdf]])
+#+name: fig:weights_NP_RS_RP
+#+caption: Upper bounds on the norm of the complementary filters for NP, RS and RP
+#+RESULTS:
 [[file:figs/weights_NP_RS_RP.png]]
 
 * H-Infinity synthesis of complementary filters
-  <<sec:h_infinity>>
+<<sec:h_infinity>>
 
 We here synthesize the complementary filters using the $\mathcal{H}_\infty$ synthesis.
 The goal is to specify upper bounds on the norms of $H_L$ and $H_H$ while ensuring their complementary property ($H_L + H_H = 1$).
@@ -374,7 +387,7 @@ In order to do so, we use the generalized plant shown on figure [[fig:sf_hinf_fi
 
 #+name: fig:sf_hinf_filters_plant_b
 #+caption: Generalized plant used for the $\mathcal{H}_\infty$ synthesis of the complementary filters
-[[file:figs-tikz/sf_hinf_filters_plant_b.png]]
+[[file:figs-paper/sf_hinf_filters_plant_b.png]]
 
 The $\mathcal{H}_\infty$ synthesis applied on this generalized plant will give a transfer function $H_L$ (figure [[fig:sf_hinf_filters_b]]) such that the $\mathcal{H}_\infty$ norm of the transfer function from $w$ to $[z_H,\ z_L]$ is less than one:
 \[ \left\| \begin{array}{c} H_L w_L \\ (1 - H_L) w_H \end{array} \right\|_\infty < 1 \]
@@ -391,7 +404,7 @@ We then see that $w_L$ and $w_H$ can be used to shape both $H_L$ and $H_H$ while
 
 #+name: fig:sf_hinf_filters_b
 #+caption: $\mathcal{H}_\infty$ synthesis of the complementary filters
-[[file:figs-tikz/sf_hinf_filters_b.png]]
+[[file:figs-paper/sf_hinf_filters_b.png]]
 
 
 Thus, if we choose $w_L$ and $w_H$ such that $1/|w_L|$ and $1/|w_H|$ lie below the upper bounds of figure [[fig:weights_NP_RS_RP]], we will ensure the NP, RS and RP of the controlled system.
@@ -405,42 +418,6 @@ Depending if we are interested only in NP, RS or RP, we can adjust the weights $
   wL = (s + omegab/(4.5)^(1/3))^3/(s*(1e-4)^(1/3) + omegab)^3;
 #+end_src
 
-#+begin_src matlab :exports none
-  wT_resp = abs(squeeze(freqresp(wT, freqs, 'Hz')));
-  wI_resp = abs(squeeze(freqresp(wI, freqs, 'Hz')));
-  wS_resp = abs(squeeze(freqresp(wS, freqs, 'Hz')));
-  Hh_resp = wT_resp.*(1 + wI_resp)./(wS_resp.*(wT_resp .* (1 + wI_resp) + wI_resp));
-  Hl_resp = 1./(wT_resp .* (1 + wI_resp) + wI_resp);
-
-  figure;
-  hold on;
-  plot(freqs, wS_resp .* Hh_resp + wI_resp .* Hl_resp, '--', 'DisplayName', 'NP - $H_L$');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]'); ylabel('Magnitude');
-  hold off;
-  xlim([freqs(1), freqs(end)]);
-  ylim([1e-3, 10]);
-  xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
-#+end_src
-
-#+begin_src matlab :exports none
-  wH_resp = abs(squeeze(freqresp(wH, freqs, 'Hz')));
-  wL_resp = abs(squeeze(freqresp(wL, freqs, 'Hz')));
-
-  figure;
-  hold on;
-  plot(freqs, wH_resp .* Hh_resp + wL_resp.*Hl_resp, '--', 'DisplayName', 'NP - $H_L$');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]'); ylabel('Magnitude');
-  hold off;
-  xlim([freqs(1), freqs(end)]);
-  ylim([1e-3, 10]);
-  xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
-#+end_src
-
-
 #+begin_src matlab :exports none
   figure;
   hold on;
@@ -464,17 +441,16 @@ Depending if we are interested only in NP, RS or RP, we can adjust the weights $
   hold off;
   xlim([freqs(1), freqs(end)]);
   ylim([1e-3, 10]);
-  xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/weights_wl_wh.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/weights_wl_wh.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:weights_wl_wh
-#+CAPTION: Weights on the complementary filters $w_L$ and $w_H$ and the associated performance weights ([[./figs/weights_wl_wh.png][png]], [[./figs/weights_wl_wh.pdf][pdf]])
+#+name: fig:weights_wl_wh
+#+caption: Weights on the complementary filters $w_L$ and $w_H$ and the associated performance weights
+#+RESULTS:
 [[file:figs/weights_wl_wh.png]]
 
 We define the generalized plant $P$ on matlab.
@@ -518,6 +494,10 @@ We then define the high pass filter $H_H = 1 - H_L$. The bode plot of both $H_L$
 
 #+begin_src matlab :exports none
   figure;
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+  % Magnitude
+  ax1 = nexttile([2, 1]);
   hold on;
   set(gca,'ColorOrderIndex',1)
   plot(freqs, 1./abs(squeeze(freqresp(wL, freqs, 'Hz'))), '--', 'DisplayName', '$w_L$');
@@ -528,26 +508,40 @@ We then define the high pass filter $H_H = 1 - H_L$. The bode plot of both $H_L$
   plot(freqs, abs(squeeze(freqresp(Hl_hinf, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$ - $\mathcal{H}_\infty$');
   set(gca,'ColorOrderIndex',2)
   plot(freqs, abs(squeeze(freqresp(Hh_hinf, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$ - $\mathcal{H}_\infty$');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]'); ylabel('Magnitude');
   hold off;
-  xlim([freqs(1), freqs(end)]);
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
   ylim([1e-3, 10]);
-  xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8);
+
+  % Phase
+  ax2 = nexttile;
+  hold on;
+  plot(freqs, 180/pi*phase(squeeze(freqresp(Hl_hinf, freqs, 'Hz'))), '-');
+  plot(freqs, 180/pi*phase(squeeze(freqresp(Hh_hinf, freqs, 'Hz'))), '-');
+  hold off;
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'XScale', 'log');
+  yticks([-360:90:360]);
+  axis padded
+
+  linkaxes([ax1,ax2],'x');
+  xlim([freqs(1), freqs(end)]);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/hinf_filters_results.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/hinf_filters_results.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:hinf_filters_results
-#+CAPTION: Obtained complementary filters using $\mathcal{H}_\infty$ synthesis ([[./figs/hinf_filters_results.png][png]], [[./figs/hinf_filters_results.pdf][pdf]])
+#+name: fig:hinf_filters_results
+#+caption: Obtained complementary filters using $\mathcal{H}_\infty$ synthesis
+#+RESULTS:
 [[file:figs/hinf_filters_results.png]]
 
-* Complementary filters using analytical formula
-  <<sec:analytical_formula>>
+* TODO Complementary filters using analytical formula
+<<sec:analytical_formula>>
+
+- [ ] The problem here is that the High pass filter gain goes to zero at high frequency, and when inverted, it will go to infinity.
 
 We here use analytical formula for the complementary filters $H_L$ and $H_H$.
 
@@ -584,8 +578,7 @@ The slope of those filters at high and low frequencies is $-2$ and $2$ respectiv
     plot(freqs_study, abs(squeeze(freqresp(Hl2, freqs_study, 'Hz'))));
   end
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude');
+  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
   hold off;
   ylim([1e-3, 20]);
   % Phase
@@ -604,22 +597,21 @@ The slope of those filters at high and low frequencies is $-2$ and $2$ respectiv
   yticks(-180:90:180);
   ylim([-180 180]);
   xlabel('Relative Frequency $\frac{\omega}{\omega_0}$'); ylabel('Phase [deg]');
-  legend('Location', 'northeast');
+  legend('Location', 'northeast', 'FontSize', 8);
   hold off;
   linkaxes([ax1,ax2],'x');
   xlim([freqs_study(1), freqs_study(end)]);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/comp_filters_param_alpha.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/comp_filters_param_alpha.pdf', 'width', 'half', 'height', 'tall');
 #+end_src
 
-#+NAME: fig:comp_filters_param_alpha
-#+CAPTION: Effect of the parameter $\alpha$ on the shape of the generated second order complementary filters ([[./figs/comp_filters_param_alpha.png][png]], [[./figs/comp_filters_param_alpha.pdf][pdf]])
+#+name: fig:comp_filters_param_alpha
+#+caption: Effect of the parameter $\alpha$ on the shape of the generated second order complementary filters
+#+RESULTS:
 [[file:figs/comp_filters_param_alpha.png]]
 
-
 The parameters $\alpha$ and $\omega_0$ are chosen in order to have that the complementary filters stay below the defined upper bounds.
 
 The obtained complementary filters are shown on figure [[fig:complementary_filters_second_order]].
@@ -657,19 +649,18 @@ The Robust Performance is not fulfilled for $T$, and we see that the RP conditio
   xlim([freqs(1), freqs(end)]);
   ylim([1e-3, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/complementary_filters_second_order.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/complementary_filters_second_order.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:complementary_filters_second_order
-#+CAPTION: Second order complementary filters using the analytical formula ([[./figs/complementary_filters_second_order.png][png]], [[./figs/complementary_filters_second_order.pdf][pdf]])
+#+name: fig:complementary_filters_second_order
+#+caption: Second order complementary filters using the analytical formula
+#+RESULTS:
 [[file:figs/complementary_filters_second_order.png]]
 
-
 The following formula gives complementary filters with slopes of $-3$ and $3$:
 \begin{align*}
   H_L(s) &= \frac{\left(1+(\alpha+1)(\beta+1)\right) (\frac{s}{\omega_0})^2 + (1+\alpha+\beta)(\frac{s}{\omega_0}) + 1}{\left(\frac{s}{\omega_0} + 1\right) \left( (\frac{s}{\omega_0})^2 + \alpha (\frac{s}{\omega_0}) + 1 \right) \left( (\frac{s}{\omega_0})^2 + \beta (\frac{s}{\omega_0}) + 1 \right)}\\
@@ -715,20 +706,20 @@ The filters are defined below and the result is shown on figure [[fig:complement
   xlim([freqs(1), freqs(end)]);
   ylim([1e-3, 10]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  legend('location', 'southeast', 'FontSize', 8);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/complementary_filters_third_order.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/complementary_filters_third_order.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:complementary_filters_third_order
-#+CAPTION: Third order complementary filters using the analytical formula ([[./figs/complementary_filters_third_order.png][png]], [[./figs/complementary_filters_third_order.pdf][pdf]])
+#+name: fig:complementary_filters_third_order
+#+caption: Third order complementary filters using the analytical formula
+#+RESULTS:
 [[file:figs/complementary_filters_third_order.png]]
 
 * Comparison of complementary filters
-  <<sec:comp_filters>>
+<<sec:comp_filters>>
 The generated complementary filters using $\mathcal{H}_\infty$ and the analytical formulas are compared on figure [[fig:comp_hinf_analytical]].
 
 Although they are very close to each other, there is some difference to note here:
@@ -739,8 +730,10 @@ The complementary filters obtained with the $\mathcal{H}_\infty$ will be used fo
 
 #+begin_src matlab :exports none
   figure;
+  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax1 = subplot(2,1,1);
+  % Magnitude
+  ax1 = nexttile;
   hold on;
   set(gca,'ColorOrderIndex',1)
   plot(freqs, abs(squeeze(freqresp(Hl_hinf, freqs, 'Hz'))), '--');
@@ -757,11 +750,12 @@ The complementary filters obtained with the $\mathcal{H}_\infty$ will be used fo
   set(gca,'ColorOrderIndex',2)
   plot(freqs, abs(squeeze(freqresp(Hh3_ana, freqs, 'Hz'))), ':');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Magnitude');
+  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
   hold off;
   ylim([1e-4, 10]);
 
-  ax2 = subplot(2,1,2);
+  % Phase
+  ax2 = nexttile;
   hold on;
   set(gca,'ColorOrderIndex',1)
   plot(freqs, 180/pi*phase(squeeze(freqresp(Hl_hinf, freqs, 'Hz'))), '--', 'DisplayName', '$H_L$ - $\mathcal{H}_\infty$');
@@ -781,24 +775,24 @@ The complementary filters obtained with the $\mathcal{H}_\infty$ will be used fo
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
   hold off;
   yticks([-360:90:360]);
-  legend('location', 'northeast');
+  legend('location', 'northeast', 'FontSize', 8);
 
   linkaxes([ax1,ax2],'x');
   xlim([freqs(1), freqs(end)]);
   xticks([0.1, 1, 10, 100, 1000]);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/comp_hinf_analytical.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/comp_hinf_analytical.pdf', 'width', 'half', 'height', 'tall');
 #+end_src
 
-#+NAME: fig:comp_hinf_analytical
-#+CAPTION: Comparison of the complementary filters obtained with $\mathcal{H}_\infty$ synthesis and with the analytical formula ([[./figs/comp_hinf_analytical.png][png]], [[./figs/comp_hinf_analytical.pdf][pdf]])
+#+name: fig:comp_hinf_analytical
+#+caption: Comparison of the complementary filters obtained with $\mathcal{H}_\infty$ synthesis and with the analytical formula
+#+RESULTS:
 [[file:figs/comp_hinf_analytical.png]]
 
 * Controller Analysis
-  <<sec:controller_analysis>>
+<<sec:controller_analysis>>
 
 The controller $K$ is computed from the plant model $G$ and the low pass filter $H_H$:
 \[ K = G^{-1} H_H^{-1} \]
@@ -806,7 +800,7 @@ The controller $K$ is computed from the plant model $G$ and the low pass filter
 As this is not proper and thus realizable, a second order low pass filter is added with a crossover frequency much larger than the control bandwidth.
 
 #+begin_src matlab
-  omega = 2*pi*1000;
+  omega = 2*pi*500;
   K = 1/(Hh_hinf*G) * 1/((1+s/omega)*(1+s/omega+(s/omega)^2));
 #+end_src
 
@@ -838,18 +832,19 @@ The bode plot of the controller is shown on figure [[fig:bode_plot_controller]]:
 
 #+begin_src matlab :exports none
   figure;
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
   % Magnitude
-  ax1 = subplot(2,1,1);
+  ax1 = nexttile([2, 1]);
   hold on;
   plot(freqs, abs(squeeze(freqresp(K, freqs, 'Hz'))), 'k-');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude [N/m]');
-  % ylim([1e3, 1e8])
+  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
+  ylim([8e3, 1e8])
   hold off;
 
   % Phase
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
   plot(freqs, 180/pi*angle(squeeze(freqresp(K, freqs, 'Hz'))), 'k-');
   set(gca,'xscale','log');
@@ -857,22 +852,22 @@ The bode plot of the controller is shown on figure [[fig:bode_plot_controller]]:
   ylim([-180 180]);
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
   hold off;
+
   linkaxes([ax1,ax2],'x');
   xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/bode_plot_controller.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/bode_plot_controller.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:bode_plot_controller
-#+CAPTION: Bode plot of the obtained controller $K$ ([[./figs/bode_plot_controller.png][png]], [[./figs/bode_plot_controller.pdf][pdf]])
+#+name: fig:bode_plot_controller
+#+caption: Bode plot of the obtained controller $K$
+#+RESULTS:
 [[file:figs/bode_plot_controller.png]]
 
 * Nominal Stability and Nominal Performance
-  <<sec:nominal_stability_performance>>
+<<sec:nominal_stability_performance>>
 
 The nominal stability of the system is first checked with the =allmargin= matlab command.
 
@@ -886,12 +881,12 @@ allmargin(K*G*Hl_hinf)
 ans =
   struct with fields:
 
-     GainMargin: 4.46426896164391
-    GMFrequency: 243.854595348016
-    PhaseMargin: 35.7045152899792
-    PMFrequency: 88.3664383511655
-    DelayMargin: 0.00705201387841809
-    DMFrequency: 88.3664383511655
+     GainMargin: 3.86196691340257
+    GMFrequency: 225.992831144248
+    PhaseMargin: 34.0890006269983
+    PMFrequency: 88.3632737129296
+    DelayMargin: 0.0067331740287082
+    DMFrequency: 88.3632737129296
          Stable: 1
 #+end_example
 
@@ -901,17 +896,20 @@ The bode plot of the loop gain $L = K*G*H_L$ is shown on figure [[fig:bode_plot_
 
 #+begin_src matlab :exports none
   figure;
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
   % Magnitude
-  ax1 = subplot(2,1,1);
+  ax1 = nexttile([2, 1]);
   hold on;
   plot(freqs, abs(squeeze(freqresp(K*G*Hl_hinf, freqs, 'Hz'))), 'k-');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude [m/N]');
+  ylabel('Loop Gain');
   hold off;
+  ylim([1e-3, 1e2])
 
   % Phase
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
   plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K*G*Hl_hinf, freqs, 'Hz')))), 'k-');
   set(gca,'xscale','log');
@@ -919,18 +917,18 @@ The bode plot of the loop gain $L = K*G*H_L$ is shown on figure [[fig:bode_plot_
   ylim([-270 0]);
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
   hold off;
+
   linkaxes([ax1,ax2],'x');
   xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/bode_plot_loop_gain.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/bode_plot_loop_gain.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:bode_plot_loop_gain
-#+CAPTION: Bode Plot of the Loop Gain $L = K G H_L$ ([[./figs/bode_plot_loop_gain.png][png]], [[./figs/bode_plot_loop_gain.pdf][pdf]])
+#+name: fig:bode_plot_loop_gain
+#+caption: Bode Plot of the Loop Gain $L = K G H_L$
+#+RESULTS:
 [[file:figs/bode_plot_loop_gain.png]]
 
 In order to check the Nominal Performance of the system, we compute the sensibility and the complementary sensibility transfer functions.
@@ -959,22 +957,21 @@ As expected, we guarantee the Nominal Performance of the system.
   xlabel('Frequency [Hz]'); ylabel('Magnitude');
   hold off;
   xlim([freqs(1), freqs(end)]);
-  ylim([1e-4, 10]);
-  xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeast');
+  ylim([1e-4, 5]);
+  legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/verification_NP.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/verification_NP.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:verification_NP
-#+CAPTION: Bode plot of $S$ and $T$ in order to verify the nominal performance of the system ([[./figs/verification_NP.png][png]], [[./figs/verification_NP.pdf][pdf]])
+#+name: fig:verification_NP
+#+caption: Bode plot of $S$ and $T$ in order to verify the nominal performance of the system
+#+RESULTS:
 [[file:figs/verification_NP.png]]
 
 * Robust Stability and Robust Performance
-  <<sec:robustness_analysis>>
+<<sec:robustness_analysis>>
 In order to verify the Robust stability of the system, we can use the following equivalence:
 \[ \text{RS} \Leftrightarrow \left| w_I T \right| < 1 \quad \forall \omega  \]
 
@@ -992,61 +989,64 @@ This is shown on figure [[fig:robust_stability]].
   plot([freqs(1) freqs(end)], [1 1], 'k--');
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]');
-  ylim([0.02, 2])
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
+  xlim([freqs(1), freqs(end)]);
+  ylim([0.02, 2]);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/robust_stability.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/robust_stability.pdf', 'width', 'half', 'height', 'small');
 #+end_src
 
-#+NAME: fig:robust_stability
-#+CAPTION: Robust Stability Check: $|w_I T| < 1, \quad \forall \omega$ ([[./figs/robust_stability.png][png]], [[./figs/robust_stability.pdf][pdf]])
-[[file:figs-tikz/robust_stability.png]]
+#+name: fig:robust_stability
+#+caption: Robust Stability Check: $|w_I T| < 1, \quad \forall \omega$
+#+RESULTS:
+[[file:figs/robust_stability.png]]
 
 To check Robust Stability, we can also look at the loop gain of the uncertain system (figure [[fig:loop_gain_robustness]]) or the Nyquist plot (figure [[fig:nyquist_robustness]]).
 
-#+begin_src matlab :results silent :exports none
+#+begin_src matlab :exports none
   figure;
-  ax2 = subplot(2,1,1);
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+  % Magnitude
+  ax1 = nexttile([2, 1]);
   hold on;
   for i=1:length(Gds)
       plot(freqs, abs(squeeze(freqresp(Gds(:, :, i)*K*Hl_hinf, freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.1]);
   end
   plot(freqs, abs(squeeze(freqresp(G*K*Hl_hinf, freqs, 'Hz'))), 'k-');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude [m/N]');
-  ylim([1e-4 1e4]);
   hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Loop Gain'); set(gca, 'XTickLabel',[]);
+  ylim([1e-3, 1e2])
+
   % Phase
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
   for i=1:length(Gds)
       plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gds(:, :, i)*K*Hl_hinf, freqs, 'Hz')))), '-', 'color', [0, 0, 0, 0.1]);
   end
   plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G*K*Hl_hinf, freqs, 'Hz')))), 'k-');
-  set(gca,'xscale','log');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   yticks(-360:90:180);
   ylim([-270 0]);
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
   hold off;
+
   linkaxes([ax1,ax2],'x');
   xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000])
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/loop_gain_robustness.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/loop_gain_robustness.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:loop_gain_robustness
-#+CAPTION: Loop Gain of the uncertain system ([[./figs/loop_gain_robustness.png][png]], [[./figs/loop_gain_robustness.pdf][pdf]])
+#+name: fig:loop_gain_robustness
+#+caption: Loop Gain of the uncertain system
+#+RESULTS:
 [[file:figs/loop_gain_robustness.png]]
 
-
 #+begin_src matlab :exports none
   freqs_nyquist = logspace(0, 4, 100);
 
@@ -1062,13 +1062,13 @@ To check Robust Stability, we can also look at the loop gain of the uncertain sy
   xlabel('Real Part'); ylabel('Imaginary Part');
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/nyquist_robustness.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/nyquist_robustness.pdf', 'width', 'half', 'height', 'small');
 #+end_src
 
-#+NAME: fig:nyquist_robustness
-#+CAPTION: Nyquist plot of the uncertain system ([[./figs/nyquist_robustness.png][png]], [[./figs/nyquist_robustness.pdf][pdf]])
+#+name: fig:nyquist_robustness
+#+caption: Nyquist plot of the uncertain system
+#+RESULTS:
 [[file:figs/nyquist_robustness.png]]
 
 The Robust Performance is verified by plotting $|S|$ and $|T|$ for the uncertain system along side the upper bounds defined for performance.
@@ -1077,44 +1077,41 @@ This is shown on figure [[fig:robust_performance_result]] and we can indeed conf
 #+begin_src matlab :exports none
   figure;
   hold on;
-  plot(freqs, abs(squeeze(freqresp(Ts(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'DisplayName', '$|T|$ - Uncertain');
-  plot(freqs, abs(squeeze(freqresp(Ss(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'DisplayName', '$|S|$ - Uncertain');
-
-  for i=2:length(Gds)
+  for i=1:length(Gds)
       plot(freqs, abs(squeeze(freqresp(Ts(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'HandleVisibility', 'off');
       plot(freqs, abs(squeeze(freqresp(Ss(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'HandleVisibility', 'off');
   end
 
   set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(G*K*Hl_hinf/(1+G*K*Hl_hinf), freqs, 'Hz'))), 'DisplayName', '$|T|$ - Nominal');
+  plot(freqs, abs(squeeze(freqresp(G*K*Hl_hinf/(1+G*K*Hl_hinf), freqs, 'Hz'))), 'DisplayName', '$|T|$');
   set(gca,'ColorOrderIndex',2)
-  plot(freqs, abs(squeeze(freqresp(1/(1+G*K*Hl_hinf), freqs, 'Hz'))), 'DisplayName', '$|S|$ - Nominal');
+  plot(freqs, abs(squeeze(freqresp(1/(1+G*K*Hl_hinf), freqs, 'Hz'))), 'DisplayName', '$|S|$');
 
   set(gca,'ColorOrderIndex',1)
-  plot([100, 1000], [0.1, 0.001], ':', 'DisplayName', '$|T|$ - Upper bound');
+  plot([100, 1000], [0.1, 0.001], ':', 'DisplayName', '$|T|$ - Spec.');
   set(gca,'ColorOrderIndex',2)
-  plot([0.1, 0.2, 2], [0.001, 0.001, 0.1], ':', 'DisplayName', '$|S|$ - Upper bound');
+  plot([0.1, 0.2, 2], [0.001, 0.001, 0.1], ':', 'DisplayName', '$|S|$ - Spec.');
 
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   hold off;
-  xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
   xlim([freqs(1), freqs(end)]);
   ylim([1e-4, 5]);
   xticks([0.1, 1, 10, 100, 1000]);
-  legend('location', 'northeastoutside');
+  legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/robust_performance_result.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/robust_performance_result.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:robust_performance_result
-#+CAPTION: Verification of the Robust Performance ([[./figs/robust_performance_result.png][png]], [[./figs/robust_performance_result.pdf][pdf]])
+#+name: fig:robust_performance_result
+#+caption: Verification of the Robust Performance
+#+RESULTS:
 [[file:figs/robust_performance_result.png]]
 
 * Pre-filter
-  <<sec:pre_filter>>
+<<sec:pre_filter>>
 
 For now, we have not specified any performance requirements on the input usage due to change of reference.
 Do limit the input usage due to change of reference, we can use a pre-filter $K_r$ as shown on figure [[fig:sf_arch_class_prefilter]].
@@ -1129,7 +1126,7 @@ We then run a simulation for a step of $1\mu m$ for the system without and with
 This confirms that a pre-filter can be used to limit the input usage due to change of reference using this architecture.
 
 #+begin_src matlab
-  t  = linspace(0, 0.02, 1000);
+  t  = linspace(0, 0.01, 500);
 
   opts = stepDataOptions;
   opts.StepAmplitude = 1e-6;
@@ -1142,204 +1139,187 @@ This confirms that a pre-filter can be used to limit the input usage due to chan
 
 #+begin_src matlab :exports none
   figure;
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(t, u, 'k--', 'DisplayName', 'Without Pre-filter');
-  plot(t, uf, 'k-', 'DisplayName', 'With Pre-Filter');
-  hold off;
-  ylabel('Command Input [N]');
-  set(gca, 'XTickLabel',[]);
-  legend('location', 'northeast');
+  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax2 = subplot(2,1,2);
+  ax1 = nexttile;
   hold on;
-  plot(t, y,  'k--');
-  plot(t, yf, 'k-' );
+  plot(t, u,  'DisplayName', 'Without Pre-filter');
+  plot(t, uf, 'DisplayName', 'With Pre-Filter');
   hold off;
-  xlabel('Time [s]');
-  ylabel('Output [m]');
+  ylabel('Input $u$ [N]'); set(gca, 'XTickLabel',[]);
+  legend('location', 'northeast', 'FontSize', 8);
+
+  ax2 = nexttile;
+  hold on;
+  plot(t, y);
+  plot(t, yf);
+  hold off;
+  xlabel('Time [s]'); ylabel('Output [m]');
 #+end_src
 
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/u_and_y_with_Kr.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/u_and_y_with_Kr.pdf', 'width', 'half', 'height', 'normal');
 #+end_src
 
-#+NAME: fig:u_and_y_with_Kr
-#+CAPTION: Input usage and response due to a step change of reference when using a pre-filter $K_r$ ([[./figs/u_and_y_with_Kr.png][png]], [[./figs/u_and_y_with_Kr.pdf][pdf]])
+#+name: fig:u_and_y_with_Kr
+#+caption: Input usage and response due to a step change of reference when using a pre-filter $K_r$
+#+RESULTS:
 [[file:figs/u_and_y_with_Kr.png]]
 
-* Controller using classical techniques
-  <<sec:sisotool_controller>>
-A controller is designed using =SISOTOOL= with a bandwidth of approximately $20\ \text{Hz}$ and with two integrator.
+* Matlab Functions
+<<sec:matlab_functions>>
+:PROPERTIES:
+:header-args:matlab+: :comments none :mkdirp yes :eval no
+:END:
 
-The obtained controller is shown below.
-#+begin_src matlab
-  Kf = 1.1814e12*(s+10.15)*(s+9.036)*(s+53.8)/(s^2*(s+216.1)*(s+1200)*(s+1864));
+** =createWeight=
+<<sec:createWeight>>
+
+This Matlab function is accessible [[file:src/createWeight.m][here]].
+
+#+begin_src matlab :tangle src/createWeight.m
+    function [W] = createWeight(args)
+    % createWeight -
+    %
+    % Syntax: [in_data] = createWeight(in_data)
+    %
+    % Inputs:
+    %    - n  - Weight Order
+    %    - G0 - Low frequency Gain
+    %    - G1 - High frequency Gain
+    %    - Gc - Gain of W at frequency w0
+    %    - w0 - Frequency at which |W(j w0)| = Gc
+    %
+    % Outputs:
+    %    - W - Generated Weight
+
+        arguments
+            args.n  (1,1) double {mustBeInteger, mustBePositive} = 1
+            args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1
+            args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10
+            args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1
+            args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1
+        end
+
+      mustBeBetween(args.G0, args.Gc, args.G1);
+
+      s = tf('s');
+
+      W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (args.G0/args.Gc)^(1/args.n))/((1/args.G1)^(1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (1/args.Gc)^(1/args.n)))^args.n;
+
+      end
+
+      % Custom validation function
+      function mustBeBetween(a,b,c)
+          if ~((a > b && b > c) || (c > b && b > a))
+              eid = 'createWeight:inputError';
+              msg = 'Gc should be between G0 and G1.';
+              throwAsCaller(MException(eid,msg))
+          end
+      end
 #+end_src
 
-#+begin_src matlab :results output :exports results replace :wrap example
-  zpk(Kf)
-#+end_src
+** =plotMagUncertainty=
+<<sec:plotMagUncertainty>>
 
-#+RESULTS:
-#+begin_example
-zpk(Kf)
+This Matlab function is accessible [[file:src/plotMagUncertainty.m][here]].
 
-ans =
+#+begin_src matlab :tangle src/plotMagUncertainty.m
+  function [p] = plotMagUncertainty(W, freqs, args)
+  % plotMagUncertainty -
+  %
+  % Syntax: [p] = plotMagUncertainty(W, freqs, args)
+  %
+  % Inputs:
+  %    - W     - Multiplicative Uncertainty Weight
+  %    - freqs - Frequency Vector [Hz]
+  %    - args  - Optional Arguments:
+  %      - G
+  %      - color_i
+  %      - opacity
+  %
+  % Outputs:
+  %    - p - Plot Handle
 
-  1.1814e12 (s+10.15) (s+9.036) (s+53.8)
-  --------------------------------------
-     s^2 (s+216.1) (s+1200) (s+1864)
-
-Continuous-time zero/pole/gain model.
-#+end_example
-
-The loop gain for both cases are compared on figure [[fig:loop_gain_compare]].
-
-#+begin_src matlab :exports none
-  figure;
-  % Magnitude
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(K*G*Hl_hinf, freqs, 'Hz'))),  'k--');
-  plot(freqs, abs(squeeze(freqresp(Kf*G, freqs, 'Hz'))), 'k-');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  set(gca, 'XTickLabel',[]);
-  ylabel('Magnitude [N/m]');
-  % ylim([1e3, 1e8])
-  hold off;
-
-  % Phase
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*angle(squeeze(freqresp(K*G*Hl_hinf, freqs, 'Hz'))), 'k--', 'DisplayName', '$K G H_L$ - $\mathcal{H}_\infty$');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(Kf*G, freqs, 'Hz'))), 'k-', 'DisplayName', '$K G$ - SISOTOOL');
-  set(gca,'xscale','log');
-  yticks(-180:90:180);
-  ylim([-180 180]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  legend('location', 'northwest');
-
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000])
-#+end_src
-
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/loop_gain_compare.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
-#+end_src
-
-#+NAME: fig:loop_gain_compare
-#+CAPTION: Comparison of the Loop Gains ([[./figs/loop_gain_compare.png][png]], [[./figs/loop_gain_compare.pdf][pdf]])
-[[file:figs/loop_gain_compare.png]]
-
-The Robust Stability of the system is verified using the Nyquist plot on figure [[fig:nyquist_plot_sisotool_controller]].
-
-#+begin_src matlab :exports none
-  freqs_nyquist = logspace(0, 4, 100);
-
-  figure;
-  hold on;
-  for i=1:length(Gds)
-      plot(real(squeeze(freqresp(Gds(:, :, i)*Kf, freqs_nyquist, 'Hz'))), imag(squeeze(freqresp(Gds(:, :, i)*Kf, freqs_nyquist, 'Hz'))), 'color', [0, 0, 0, 0.1]);
-  end
-  plot(real(squeeze(freqresp(G*Kf, freqs_nyquist, 'Hz'))), imag(squeeze(freqresp(G*Kf, freqs_nyquist, 'Hz'))), 'k');
-  hold off;
-  xlim([-1.4, 0.2]); ylim([-1.4, 0.2]);
-  xticks(-1.4:0.2:0.2); yticks(-1.4:0.2:0.2);
-  xlabel('Real Part'); ylabel('Imaginary Part');
-#+end_src
-
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/nyquist_plot_sisotool_controller.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
-#+end_src
-
-#+NAME: fig:nyquist_plot_sisotool_controller
-#+CAPTION: Nyquist Plot of the uncertain system ([[./figs/nyquist_plot_sisotool_controller.png][png]], [[./figs/nyquist_plot_sisotool_controller.pdf][pdf]])
-[[file:figs/nyquist_plot_sisotool_controller.png]]
-
-The closed loop sensitivity and complementary sensitivity transfer functions are computed.
-And finally, the Robust Performance of both systems are compared on figure [[fig:robust_performance_compare]].
-
-#+begin_src matlab :exports none
-  Sf = 1/(Kf*G + 1);
-  Tf = Kf*G/(Kf*G + 1);
-
-  Tfs = Gds*Kf/(Gds*Kf + 1);
-  Sfs = 1/(Gds*Kf + 1);
-#+end_src
-
-
-#+begin_src matlab :exports none
-  figure;
-  ax1 = subplot(1, 2, 1);
-  title('$K$ - SISOTOOL');
-  hold on;
-
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(Tf, freqs, 'Hz'))), 'DisplayName', '$|T|$ - Nominal');
-  set(gca,'ColorOrderIndex',2)
-  plot(freqs, abs(squeeze(freqresp(Sf, freqs, 'Hz'))), 'DisplayName', '$|S|$ - Nominal');
-
-  plot(freqs, abs(squeeze(freqresp(Tfs(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'DisplayName', '$|T|$ - Uncertain');
-  plot(freqs, abs(squeeze(freqresp(Sfs(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'DisplayName', '$|S|$ - Uncertain');
-
-  for i=2:length(Gds)
-      plot(freqs, abs(squeeze(freqresp(Tfs(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'HandleVisibility', 'off');
-      plot(freqs, abs(squeeze(freqresp(Sfs(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'HandleVisibility', 'off');
+      arguments
+      W
+      freqs double {mustBeNumeric, mustBeNonnegative}
+      args.G = tf(1)
+      args.color_i (1,1) double {mustBeInteger, mustBeNonnegative} = 0
+      args.opacity (1,1) double {mustBeNumeric, mustBeNonnegative} = 0.3
+      args.DisplayName char = ''
   end
 
-  set(gca,'ColorOrderIndex',1)
-  plot([100, 1000], [0.1, 0.001], ':', 'DisplayName', '$|T|$ - Upper bound');
-  set(gca,'ColorOrderIndex',2)
-  plot([0.1, 0.2, 2], [0.001, 0.001, 0.1], ':', 'DisplayName', '$|S|$ - Upper bound');
+  % Get defaults colors
+  colors = get(groot, 'defaultAxesColorOrder');
 
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  hold off;
-  xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
-  xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000]);
+  p = patch([freqs flip(freqs)], ...
+            [abs(squeeze(freqresp(args.G, freqs, 'Hz'))).*(1 + abs(squeeze(freqresp(W, freqs, 'Hz')))); ...
+             flip(abs(squeeze(freqresp(args.G, freqs, 'Hz'))).*max(1 - abs(squeeze(freqresp(W, freqs, 'Hz'))), 1e-6))], 'w', ...
+            'DisplayName', args.DisplayName);
 
-  ax2 = subplot(1, 2, 2);
-  title('$K$ - complementary filters');
-  hold on;
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(T, freqs, 'Hz'))), 'DisplayName', '$|T|$ - Nominal');
-  set(gca,'ColorOrderIndex',2)
-  plot(freqs, abs(squeeze(freqresp(S, freqs, 'Hz'))), 'DisplayName', '$|S|$ - Nominal');
+  if args.color_i == 0
+      p.FaceColor = [0; 0; 0];
+  else
+      p.FaceColor = colors(args.color_i, :);
+  end
+  p.EdgeColor = 'none';
+  p.FaceAlpha = args.opacity;
 
-  plot(freqs, abs(squeeze(freqresp(Ts(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'DisplayName', '$|T|$ - Uncertain');
-  plot(freqs, abs(squeeze(freqresp(Ss(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'DisplayName', '$|S|$ - Uncertain');
+   end
+#+end_src
 
-  for i=2:length(Gds)
-      plot(freqs, abs(squeeze(freqresp(Ts(:, :, i), freqs, 'Hz'))), 'color', [0, 0.4470, 0.7410, 0.1]     , 'HandleVisibility', 'off');
-      plot(freqs, abs(squeeze(freqresp(Ss(:, :, i), freqs, 'Hz'))), 'color', [0.8500, 0.3250, 0.0980, 0.1], 'HandleVisibility', 'off');
+** =plotPhaseUncertainty=
+<<sec:plotPhaseUncertainty>>
+
+This Matlab function is accessible [[file:src/plotPhaseUncertainty.m][here]].
+
+#+begin_src matlab :tangle src/plotPhaseUncertainty.m
+  function [p] = plotPhaseUncertainty(W, freqs, args)
+  % plotPhaseUncertainty -
+  %
+  % Syntax: [p] = plotPhaseUncertainty(W, freqs, args)
+  %
+  % Inputs:
+  %    - W     - Multiplicative Uncertainty Weight
+  %    - freqs - Frequency Vector [Hz]
+  %    - args  - Optional Arguments:
+  %      - G
+  %      - color_i
+  %      - opacity
+  %
+  % Outputs:
+  %    - p - Plot Handle
+
+  arguments
+      W
+      freqs double {mustBeNumeric, mustBeNonnegative}
+      args.G = tf(1)
+      args.color_i (1,1) double {mustBeInteger, mustBeNonnegative} = 0
+      args.opacity (1,1) double {mustBeNumeric, mustBePositive} = 0.3
+      args.DisplayName char = ''
   end
 
-  set(gca,'ColorOrderIndex',1)
-  plot([100, 1000], [0.1, 0.001], ':', 'DisplayName', '$|T|$ - Upper bound');
-  set(gca,'ColorOrderIndex',2)
-  plot([0.1, 0.2, 2], [0.001, 0.001, 0.1], ':', 'DisplayName', '$|S|$ - Upper bound');
+  % Get defaults colors
+  colors = get(groot, 'defaultAxesColorOrder');
 
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  hold off;
-  xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
-  xlim([freqs(1), freqs(end)]);
-  xticks([0.1, 1, 10, 100, 1000]);
+  % Compute Phase Uncertainty
+  Dphi = 180/pi*asin(abs(squeeze(freqresp(W, freqs, 'Hz'))));
+  Dphi(abs(squeeze(freqresp(W, freqs, 'Hz'))) > 1) = 360;
 
-  linkaxes([ax1 ax2], 'y')
-  ylim([1e-4, 10]);
+  % Compute Plant Phase
+  G_ang = 180/pi*angle(squeeze(freqresp(args.G, freqs, 'Hz')));
+
+  p = patch([freqs flip(freqs)], [G_ang+Dphi; flip(G_ang-Dphi)], 'w', ...
+            'DisplayName', args.DisplayName);
+
+  if args.color_i == 0
+      p.FaceColor = [0; 0; 0];
+  else
+      p.FaceColor = colors(args.color_i, :);
+  end
+  p.EdgeColor = 'none';
+  p.FaceAlpha = args.opacity;
+
+  end
 #+end_src
-
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/robust_performance_compare.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
-#+end_src
-
-#+NAME: fig:robust_performance_compare
-#+CAPTION: Comparison of the Robust Performance for both controllers ([[./figs/robust_performance_compare.png][png]], [[./figs/robust_performance_compare.pdf][pdf]])
-[[file:figs/robust_performance_compare.png]]
diff --git a/paper/paper.org b/paper/paper.org
index 9da55e2..519bc79 100644
--- a/paper/paper.org
+++ b/paper/paper.org
@@ -427,7 +427,7 @@ With the $\mathcal{H}_\infty$ synthesis the condition eqref:eq:hinf_problem only
                                                              &\Rightarrow |w_L H_L| + |w_H H_H| \le \sqrt{2} \quad \forall\omega
 \end{align*}
 
-And thus we have almost robust stability.
+And thus we have almost robust performance.
 
 ** Choice of the weighting functions
    <<sec:hinf_weighting_func>>