diff --git a/robust/mult_uncertainty.m b/robust/mult_uncertainty.m
index b5e565b..f0b9faa 100644
--- a/robust/mult_uncertainty.m
+++ b/robust/mult_uncertainty.m
@@ -30,17 +30,18 @@ I_m_analytic = tf(fitmagfrd(data, 2));
 % For MIMO systems, the division in the expression for I_m doesn't make
 % sense
 % This is another way to compute it:
-opt = wcOptions('MussvOptions','m3'); % for better lower bound
-[wcg, wcu, info_num] = wcgain(G-G_nom, w);
-[wcg, wcu, info_den] = wcgain(G_nom, w, opt);
 
-resp_im = info_num.Bounds(:, 2) ./ info_den.Bounds(:, 1);
+% generate some instances of the uss
+systems = usample(G, 20);
+data = frd(systems, w);
 
-data_2 = frd(resp_im, w);
-% constrain the uncertainty to stay above the magnitude data
-Constraint.LowerBound = data_2;
-Constraint.UpperBound = [];
-I_m_analytic_2 = tf(fitmagfrd(data_2, 2, [], [], Constraint));
+% cover the frequency response data using an order 2 multiplicative
+% uncertainty
+% ucover returns [G, Info] where G is the uncertain system given by
+% G(s) = [I+W1(s) Delta(s)] G_nom(s)
+% W1(s) can be recovered from Info.W1
+[G2, Info] = ucover(data, G_nom, 2);
+I_m_analytic_2 = tf(Info.W1);
 
 sigma(I_m, w); hold;
 semilogx(w, 20*log10(info.Bounds(:, 2)), 'green');