A0 = 10;
b = 1 / A0;
R1 = 10000;
R2 = R1 * (1/b - 1);
K = R1/(R1+R2);
C = [1:.2:3]*1e-12;
a0 = 1e5;
w1 = 1e4;
w2 = 1e6;
s = tf('s');
a = a0/(1+s/w1)/(1+s/w2);
A = feedback(a,b);

b_array = arrayfun (@(C) tf ([K*R2*C, K], [K*R2*C, 1]), C, 'uniformoutput', false);

A_array = cellfun (@feedback, {a}, b_array, 'uniformoutput', false);

L_array = cellfun (@mtimes, {a}, b_array, 'uniformoutput', false);

[Gm,Pm,Wcg,Wcp] = cellfun (@margin, L_array);

figure 1
step(A,'r',A_array{:});
figure 2
bode(A,A_array{:})
figure 3
plot(C,Pm)
