function [gain,p_angle] = computeBode(G_jw,w)
    %input argument as a symfun
    %this function accepts transfer function G(jw) as input argument and
    %returns two vectors as a function of w (freq in rad/s):
    %1: |G(jw)| in dB
    %2: phase angle phi in degrees (neg if output lags behind input)
    
    G_vector = double(G_jw(1i.*w));

    magnitude = abs(G_vector); %finds the magnitude of G_jw
    gain = 20*log10(magnitude); %takes the log10 of magnitude and multiplies by 20
    offset = angle(G_vector(1));
    p_angle = (angle(G_vector./exp(1i*offset)) + offset).*(180/pi);
%      p_angle = (angle(G_vector).*180./pi);

    subplot(2,1,1)
    plot(gain)
    ylabel("Gain (dB)")
    xlabel("Frequency (rad/s)")
    grid

    subplot(2,1,2)
    plot(p_angle)
    ylabel("Phase angle (deg)")
    xlabel("Frequency (rad/s)")
    grid

    clc %clears command window
    
end