function varargout = gangof4(P, C, varargin)
global AM_bode_linewidth AM_ref_linewidth

% Compute the senstivity functions
L = C*P;
S = feedback(1, L);
T = feedback(L, 1);

% Determine what to do based on the number of input and output arguments
if (nargout == 0)
  % Plot the results in standard format
  location = {321, 322, 323, 324};
  system = {T, P*S, C*S, S};
  label = {'|T|', '|PS|', '|CS|', '|S|'};

  for i=1:4
    subplot(location{i});
    [mag, phase, omega] = bode(system{i}, varargin{:});
    loglog([min(omega), max(omega)], [1 1], 'k', 'Linewidth', AM_ref_linewidth);
    hold on;
    varargout{1} = loglog(omega, squeeze(mag), 'linewidth', AM_bode_linewidth); 

    % Reset the axes to something reasonable
    omega_min = min(omega);
    omega_max = max(omega);
    mag_min = 10^floor(log10(min(mag)));
    mag_max = 10^ceil(log10(max(mag)));

    amaxis([omega_min, omega_max, mag_min, mag_max]);
    xlabel('Frequency [rad/s]');
    ylabel(label{i});

    % Set the tick labels at reasonable powers of 10
    mag_incr = ceil((log10(mag_max) - log10(mag_min))/3);
    set(gca, 'YTick', 10.^(log10(mag_min):mag_incr:log10(mag_max)));

    % Set the frequency labels at reasonable powers of 10
    omega_incr = ceil((log10(omega_max) - log10(omega_min))/4);
    set(gca, 'XTick', 10.^(log10(omega_min):omega_incr:log10(omega_max)));
  end
end

if (nargout == 5)
  % Handle Karl's version
  W = logspace(varargin{1}, varargin{2}, varargin{3});
  varargout{1} = squeeze(freqresp(T, W));
  varargout{2} = squeeze(freqresp(P*S, W));
  varargout{3} = squeeze(freqresp(C*S, W));
  varargout{4} = squeeze(freqresp(S, W));
  varargout{5} = W;
end
