diff --git a/test-bench-flexible-joints.org b/test-bench-flexible-joints.org
index a7599c3..53f4131 100644
--- a/test-bench-flexible-joints.org
+++ b/test-bench-flexible-joints.org
@@ -702,9 +702,6 @@ Now we can estimate the stiffness with a linear fit.
 sprintf('Stiffness is %.3f [N/mm]', abs(1e-3*(d\F)))
 #+end_src
 
-#+RESULTS:
-: Stiffness is 0.039 [N/mm]
-
 This is very close to the 0.04 [N/mm] written in the [[file:doc/tmp3m0cvmue_7888038c-cdc8-48d8-a837-35de02760685.pdf][Millimar 1318 probe datasheet]].
 
 And compare the linear fit with the raw measurement data (Figure [[fig:mahr_stiffness_f_d_plot]]).
@@ -735,3 +732,359 @@ The Millimar 1318 probe has a stiffness of $\approx 0.04\,[N/mm]$.
 
 * Bending Stiffness Measurement
 <<sec:bending_stiffness_meas>>
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('./matlab/mat/');
+addpath('./matlab/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('./mat/');
+#+end_src
+
+** Results
+#+begin_src matlab
+load('meas_stiff_flex_12_x.mat', 't', 'd', 'F');
+#+end_src
+
+#+begin_src matlab
+d = d(t > 5.35 & t < 14.0);
+F = F(t > 5.35 & t < 14.0);
+t = t(t > 5.35 & t < 14.0);
+
+d = d - d(1);
+F = F - F(1);
+t = t - t(1);
+#+end_src
+
+#+begin_src matlab :exports none
+%% Time Domain plots
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(t, F);
+ylabel('Force [N]'); set(gca, 'XTickLabel',[]);
+
+ax2 = nexttile;
+plot(t, d);
+xlabel('Time [s]'); ylabel('Displacement [m]');
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(F, d, 'DisplayName', 'Raw data');
+% plot(F, F/(d\F), 'DisplayName', 'Linear fit');
+hold off;
+xlabel('Measured Force [N]');
+ylabel('Measured Displacement [m]');
+legend('location', 'southeast');
+#+end_src
+
+#+begin_src matlab
+ta = [0.05, 1.13];
+tb = [5.75, 7.2];
+
+d_l = [d(t > ta(1) & t < ta(2)); d(t > tb(1) & t < tb(2))];
+F_l = [F(t > ta(1) & t < ta(2)); F(t > tb(1) & t < tb(2))];
+
+d_l = d_l - d_l(1);
+F_l = F_l - F_l(1);
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(F_l, d_l, '.', 'DisplayName', 'Raw data');
+plot(F_l, F_l/(d_l\F_l), 'DisplayName', 'Linear fit');
+hold off;
+xlabel('Measured Force [N]');
+ylabel('Measured Displacement [m]');
+legend('location', 'southeast');
+#+end_src
+
+#+begin_src matlab :results value replace :exports results :tangle no
+h = 25e-3;
+sprintf('Stiffness is %.3f [Nm/rad]', abs(h^2*(d_l\F_l)))
+#+end_src
+
+#+RESULTS:
+: Stiffness is 5.579 [Nm/rad]
+
+** Results - Y
+#+begin_src matlab
+load('meas_stiff_flex_12_y.mat', 't', 'd', 'F');
+#+end_src
+
+#+begin_src matlab
+%% Automatic Zero of the force
+F = F - mean(F(t > 0.9 & t < 1.1));
+
+%% Start measurement at t = 1.0 s
+d = d(t > 1.0);
+F = F(t > 1.0);
+t = t(t > 1.0); t = t - t(1);
+#+end_src
+
+#+begin_src matlab :exports none
+%% Time Domain plots
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(t, F);
+ylabel('Force [N]'); set(gca, 'XTickLabel',[]);
+
+ax2 = nexttile;
+plot(t, d);
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab
+t_l = [5.58, 6.75]; % Time of flexible joint's linear region
+t_s = [6.9, 7.24]; % Time of stop's linear region
+#+end_src
+
+#+begin_src matlab :exports none
+%% Linear Fit
+h = 25e-3;
+
+F_l = F(t > t_l(1) & t < t_l(2));
+d_l = d(t > t_l(1) & t < t_l(2));
+
+F_s = F(t > t_s(1) & t < t_s(2));
+d_s = d(t > t_s(1) & t < t_s(2));
+
+fit_l = polyfit(F_l, d_l, 1);
+fit_s = polyfit(F_s, d_s, 1);
+
+%% Reset displacement based on fit
+d = d - fit_l(2);
+fit_s(2) = fit_s(2) - fit_l(2);
+fit_l(2) = 0;
+
+%% Estimated Stroke
+F_max = fit_s(2)/(fit_l(1) - fit_s(1));
+d_max = fit_l(1)*F_max;
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+% plot(F_l, d_l, 'k.', 'DisplayName', 'Raw Data');
+% plot(F_s, d_s, 'k.', 'HandleVisibility', 'off');
+plot(F, d, 'k.', 'DisplayName', 'Raw data');
+set(gca,'ColorOrderIndex',1)
+plot(F_l, fit_l(1)*F_l + fit_l(2), '--', 'DisplayName', sprintf('$k_{R_x} = %.1f [Nm/rad]$', (h)^2/fit_l(1)));
+plot(F_s, fit_s(1)*F_s + fit_s(2), '--', 'DisplayName', sprintf('$k_s = %.1f [Nm/rad]$', (h)^2/fit_s(1)));
+plot([0.8*F_max, 1.2*F_max], [d_max, d_max], '--', 'DisplayName', sprintf('$R_{x,max} = %.1f [mrad]$', 1e3*atan2(d_max,h)));
+hold off;
+xlabel('Measured Force [N]');
+ylabel('Measured Displacement [m]');
+legend('location', 'southeast');
+ylim([-1e-4,inf])
+#+end_src
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stiffness is %.1f [Nm/rad]', (h)^2/fit_l(1))
+#+end_src
+
+#+RESULTS:
+: Bending Stiffness is 5.5 [Nm/rad]
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stroke is %.1f [mrad]', 1e3*atan2(d_max,h))
+#+end_src
+
+#+RESULTS:
+: Bending Stroke is 17.9 [mrad]
+
+** Results - X
+#+begin_src matlab
+load('meas_stiff_flex_12_x.mat', 't', 'd', 'F');
+#+end_src
+
+#+begin_src matlab
+%% Automatic Zero of the force
+F = F - mean(F(t > 0.9 & t < 1.1));
+
+%% Start measurement at t = 1.0 s
+d = d(t > 1.0);
+F = F(t > 1.0);
+t = t(t > 1.0); t = t - t(1);
+#+end_src
+
+#+begin_src matlab :exports none
+%% Time Domain plots
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(t, F);
+ylabel('Force [N]'); set(gca, 'XTickLabel',[]);
+
+ax2 = nexttile;
+plot(t, d);
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab
+t_l = [4.365, 5.47]; % Time of flexible joint's linear region
+t_s = [5.65, 6.09]; % Time of stop's linear region
+#+end_src
+
+#+begin_src matlab :exports none
+%% Linear Fit
+h = 25e-3;
+
+F_l = F(t > t_l(1) & t < t_l(2));
+d_l = d(t > t_l(1) & t < t_l(2));
+
+F_s = F(t > t_s(1) & t < t_s(2));
+d_s = d(t > t_s(1) & t < t_s(2));
+
+fit_l = polyfit(F_l, d_l, 1);
+fit_s = polyfit(F_s, d_s, 1);
+
+%% Reset displacement based on fit
+d = d - fit_l(2);
+fit_s(2) = fit_s(2) - fit_l(2);
+fit_l(2) = 0;
+
+%% Estimated Stroke
+F_max = fit_s(2)/(fit_l(1) - fit_s(1));
+d_max = fit_l(1)*F_max;
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(F, d, 'k.', 'DisplayName', 'Raw data');
+set(gca,'ColorOrderIndex',1)
+plot(F_l, fit_l(1)*F_l + fit_l(2), '-', 'DisplayName', sprintf('$k_{R_x} = %.1f [Nm/rad]$', (h)^2/fit_l(1)));
+plot(F_s, fit_s(1)*F_s + fit_s(2), '-', 'DisplayName', sprintf('$k_s = %.1f [Nm/rad]$', (h)^2/fit_s(1)));
+plot([0.8*F_max, 1.2*F_max], [d_max, d_max], '-', 'DisplayName', sprintf('$R_{x,max} = %.1f [mrad]$', 1e3*atan2(d_max,h)));
+hold off;
+xlabel('Measured Force [N]');
+ylabel('Measured Displacement [m]');
+legend('location', 'southeast');
+ylim([-1e-4,inf])
+#+end_src
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stiffness is %.1f [Nm/rad]', (h)^2/fit_l(1))
+#+end_src
+
+#+RESULTS:
+: Bending Stiffness is 5.7 [Nm/rad]
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stroke is %.1f [mrad]', 1e3*atan2(d_max,h))
+#+end_src
+
+#+RESULTS:
+: Bending Stroke is 17.9 [mrad]
+
+** Results - XY
+#+begin_src matlab
+load('meas_stiff_flex_12_xy.mat', 't', 'd', 'F');
+#+end_src
+
+#+begin_src matlab
+%% Automatic Zero of the force
+F = F - mean(F(t > 0.9 & t < 1.1));
+
+%% Start measurement at t = 1.0 s
+d = d(t > 1.0);
+F = F(t > 1.0);
+t = t(t > 1.0); t = t - t(1);
+#+end_src
+
+#+begin_src matlab :exports none
+%% Time Domain plots
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+plot(t, F);
+ylabel('Force [N]'); set(gca, 'XTickLabel',[]);
+
+ax2 = nexttile;
+plot(t, d);
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab
+t_l = [4.99, 6.5]; % Time of flexible joint's linear region
+t_s = [6.8, 7.10]; % Time of stop's linear region
+#+end_src
+
+#+begin_src matlab :exports none
+%% Linear Fit
+h = 25e-3;
+
+F_l = F(t > t_l(1) & t < t_l(2));
+d_l = d(t > t_l(1) & t < t_l(2));
+
+F_s = F(t > t_s(1) & t < t_s(2));
+d_s = d(t > t_s(1) & t < t_s(2));
+
+fit_l = polyfit(F_l, d_l, 1);
+fit_s = polyfit(F_s, d_s, 1);
+
+%% Reset displacement based on fit
+d = d - fit_l(2);
+fit_s(2) = fit_s(2) - fit_l(2);
+fit_l(2) = 0;
+
+%% Estimated Stroke
+F_max = fit_s(2)/(fit_l(1) - fit_s(1));
+d_max = fit_l(1)*F_max;
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(F, d, 'k.', 'DisplayName', 'Raw data');
+set(gca,'ColorOrderIndex',1)
+plot(F_l, fit_l(1)*F_l + fit_l(2), '-', 'DisplayName', sprintf('$k_{R_x} = %.1f [Nm/rad]$', (h)^2/fit_l(1)));
+plot(F_s, fit_s(1)*F_s + fit_s(2), '-', 'DisplayName', sprintf('$k_s = %.1f [Nm/rad]$', (h)^2/fit_s(1)));
+plot([0.8*F_max, 1.2*F_max], [d_max, d_max], '-', 'DisplayName', sprintf('$R_{x,max} = %.1f [mrad]$', 1e3*atan2(d_max,h)));
+hold off;
+xlabel('Measured Force [N]');
+ylabel('Measured Displacement [m]');
+legend('location', 'southeast');
+ylim([-1e-4,inf])
+#+end_src
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stiffness is %.1f [Nm/rad]', (h)^2/fit_l(1))
+#+end_src
+
+#+RESULTS:
+: Bending Stiffness is 5.6 [Nm/rad]
+
+#+begin_src matlab :results value replace :exports results :tangle no
+sprintf('Bending Stroke is %.1f [mrad]', 1e3*atan2(d_max,h))
+#+end_src
+
+#+RESULTS:
+: Bending Stroke is 23.1 [mrad]