Try this:
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(real(Transfer_func), [0 5E+4*pi], '--')
hold on
fplot(imag(Transfer_func), [0 5E+4*pi], '--')
fplot(abs(Transfer_func), [0 5E+4*pi])
hold off
legend('\Re', '\Im', '|H( j\omega )|')
title('Amplitude')
grid
subplot(2,1,2)
fplot(angle(Transfer_func), [0 5E+4*pi])
grid
title('Phase')
xlabel('Frequency (radians)')
It then uses fplot to evaluate the second version of the function, producing:
.
