function bode3(sys, range, max, N)
% BODE3(SYS, RANGE, MAX, N) plots a 3D magnitude response in the s-plane.
% RANGE is the range in [Hz] of sigma, omega axis
% MAX is the maximum magnitude (default: 3)
% N is the number of points (along each axis, default: 50)

if nargin < 4
    N = 50;
end
if nargin < 3
    max = 3;
end

% s-matrix (value of s in s-plane)
f = linspace(-range, range, N);
w = 2*pi*f;
s = zeros(N, N);                        % allocate matrix

for i = 1:N                             % populate
    s(i, :) = w;
end
ww = 1i * w';
for i = 1:N
    s(:, i) = s(:, i) + ww;
end

% magnitude response in s-plane, m
[ num, den ] = tfdata(sys, 'v');

n = polyval(num, s);
d = polyval(den, s);

m = abs(n ./ d);
m(m > max) = max;

mesh(f, f, abs(m));                     % plot

% frequency response on j-omega axis
hold on;
h = freqresp(sys, w);
h = abs(h(1, :));
h(h > max) = max;
plot3(0*f, f, h, 'r');

% axis labels
view(37.5, 30);
ylabel('Frequency [Hz]');
xlabel('Sigma [Hz]');
zlabel('Magnitude [linear]');
title('Magnitude Response (s-plane)');

% dome
hold off

