[Matlab] Bode plot without Control Toolbox

When it comes to Bode plot, it is easy to draw a Bode plot with control toolbox, but Not everybody can get this toolbox.
For those who don't have Control Toolbox, let's see how to draw a Bode plot with only basic Matlab functions.

Assuming that a transfer function is as below.
     1
  -------
   s + 1

Then, the magnitude and phase of the transfer function are neccsary for bode plot.
The magnitude and phase can be calculated by replacing s with jw.
where, j is the imaginary unit of complex number, and w is the frequency.

     1                                  1
  -------          -->         -----------
   s + 1                          j*w + 1

Possible matlab code for bode plot is as below.
--------------------------------------------------------
% Define frequency
freq = 0.01:0.01:1000;

% Transfer function
tf = 1 ./ (1 + i*freq);

% Magnitude
m = 20 * log10(abs(tf));

% Phase
phase = angle(tf);

% Plot
subplot(2,1,1)
semilogx(freq,m);
grid on
ylabel('Magnitude (dB)');

subplot(2,1,2)
semilogx(freq,phase);
grid on
ylabel('Phase (deg)');
xlabel('Frequency (rad/sec)');
--------------------------------------

As a result of the code above,

-----------------------------------------------------------------------

General method for Bode plot

% Define the numerator and denominator

num=[-0.1, -2.4, -181, -1950];

den = [1, 3.3, 990,2600];

tnum = [];

tden = [];

an = length(num);

bn = length(den);

% Numerator

for N = 1:1:an

    if num(N) > 0

        si ='+';

    elseif num(N) <0 p="">

        si = '';

    end

        

    temp = sprintf('%s%d%s%s%d ',si,num(N),'*','s.^',an-N);

    tnum = [tnum temp];

end

% Denominator

for N = 1:1:length(den)

    if den(N) > 0

        si ='+';

    elseif den(N) <0 p="">

        si = '-';

    end

        

    temp = sprintf('%s%d%s%s%d ',si,den(N),'*','s.^',bn-N);

    tden = [tden temp];

end

% Define a system transfer function

sys = @(s) eval(tnum)./ eval(tden);

% Define frequency as x 

w=0.01:0.01:1000;

% Magnitude (dB)

amp=20*log10(abs(sys(w*1i)));

% Phase

ang=angle(sys(w*1i))*180/pi;

% Plot using subplot

subplot(211)

semilogx(w,amp);

ylabel('Magnitude (dB)')

grid on

subplot(212)

semilogx(w,ang);

xlabel('frequency (rad/sec)')

ylabel('phase (deg)')

grid on