clear, close, echo off, clc,
%
% + o-------/\/\/\/\-------o-------/\/\/\/\-------o-------o +
%              R1          |          R2          |
%                          |                      |
% Vout                  ======= C1             ======= C2   Vout
%                          |                      |
%                          |                      |
% - o----------------------o----------------------o-------o -
%
% Underdamped
R1 = 330e3; R2 = 100e3; C1 = .1e-6; C2 = .1e-6;

tau1 = R1*C1
tau2 = R2*C2

zeta = .5*(R1*C1 + R1*C2 + R2*C2)
wn = sqrt(1/(R1*R2*C1*C2))

tau = inv(zeta*wn)

A = [0 1; -wn^2, -2*zeta*wn];
B = [0,0]'; C = [1,0]; D = 0;
sys1 = ss(A,B,C,D);
y0 = [1,0];
[y1, t1, x1] = initial(sys1,y0,4*tau);;

% Critically Damped
% R1 = (2 - R2*C2)/(C1 + C2)
% tau1 = R1*C1
% tau2 = R2*C2
% zeta = .5*(R1*C1 + R1*C2 + R2*C2)
% wn = sqrt(1/(R1*R2*C1*C2))
% 
% tau = inv(zeta*wn)
% 
% sys2 = tf([wn^2],[1 2*zeta*wn wn^2]);
% [y2, t2] = step(sys2,3);

% Overdamped

figure(gcf)
%subplot(2,1,1)
plot(t1,y1,t1,zeros(1,length(t1)),'--')
%title('Underdamped System Response'); grid on
axis([0 round(max(t1)) -1 1]);
xlabel('time [s]')
ylabel('Amplitude')
% subplot(2,1,2)
% plot(t2,y2), title('Critically Damped'), axis([0,3,0,1.25]); grid on
% subplot(3,1,3)
% plot(t3,y2), title('Overdamped')