MATLAB: How to manually replicate the bode() gain plot from a transfer function

bode plotcontrolControl System Toolboxtheory

Hello,
I have been trying to replicate a bode gain plot from a given transfer function (see code below).
I am not seeing where I went wrong, but my manual plots looks like this:
But the actual bode plots should look like this:
I thought that I was setting up my amplitude equation exactly the same as how the bode() function operates, but apparently not. Can you see where I made the mistake?
Here is my code:
syms C R_1 R_2 R_3 L s V_o i omega
%Equations
R_1 = 10;
R_2 = 10;
R_3 = 10;
L = 0.001; 
C = 2*10^-6;
Z_1 = R_1+L*s;
Z_2 = 1/(C*s)+R_2;
Z_3 = R_3;
Z_23 = 1/(1/(Z_2)+1/(Z_3));
Z_tot = Z_1+Z_23;
V_i = Z_1*i+V_o;
V_o = (V_i/(Z_1))/((1/((1/(C*s))+R_2))+1/R_3+1/(Z_1));
Transfer_func(s) = vpa(simplify(V_o/V_i),5);
amp = (sqrt((5.903*10^24)^2+(1.1806*10^20*omega)^2))/(sqrt((1.1806*10^25+(-2.3613*10^16*omega^2))^2+(9.4447*10^20*omega)^2))
amp_dB = 20*log(amp)
fplot(amp_dB,[0.001 1*10^6])
figure
subplot(2,1,1)
fplot(20*log(abs(Transfer_func)), [0.001 1*10^6])
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi()*angle(Transfer_func), [0.001 1*10^6])
grid
title('Phase (degrees)')
xlabel('Frequency (rad/s)')
legend('\phi( j\omega )')
H = tf([118059162071741125000 5902958103587056517120000],[23611832414348225 944473296573929026712 11805916207174113034240000]);
figure
bode(H,{0.001,1*10^6})

Best Answer

‘But the actual bode plots should look like this:’
It does.
Harkening back to yesterday: 
syms C R_1 R_2 R_3 L s V_o i omega
%Equations
R_1 = 10;
R_2 = 10;
R_3 = 10;
L = 0.001; 
C = 2*10^-6;
Z_1 = R_1+L*s;
Z_2 = 1/(C*s)+R_2;
Z_3 = R_3;
Z_23 = 1/(1/(Z_2)+1/(Z_3));
Z_tot = Z_1+Z_23;
V_i = Z_1*i+V_o;
V_o = (V_i/(Z_1))/((1/((1/(C*s))+R_2))+1/R_3+1/(Z_1));
Transfer_func(s) = vpa(simplify(V_o/V_i), 5);
Transfer_func(omega) = subs(Transfer_func, {s},{1j*omega});
figure
subplot(2,1,1)
fplot(20*log10(abs(Transfer_func)), [0.001 1E6])
set(gca, 'XScale','log')
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi*angle(Transfer_func), [0.001 1E6])
set(gca, 'XScale','log')
grid
title('Phase (degrees)')
xlabel('Frequency (rad/s)')
legend('\phi( j\omega )')
produces:
Note the slight changes between your code and mine.
Related Question