How to model transfer functions using Python (without using higher level packages) and plot the magnitude and phase response?
If the response is real you need not convert it to complex in order to compute the magnitude, but the process is correct for complex and real number since the abs function will use complex conjugate product prior to taking the square root.  If it was real, there would be a trivial phase response since the phase would always be 0 degrees. This would not be a causal system so not a realistic example to have a real frequency response. For complex cases the approach used is correct for magnitude and phase, just not the most efficient.
To confirm results (or to compute directly), use the freqs and freqz commands that are part of scipy.signal for continuous-domain (s) and discrete-domain (z) transfer functions.
For example the transfer function
$$H(s) = \frac{s+1}{s+4s+4}$$ could be plotted using:
w,h = scipy.signal.freqs([1, 1], [1, 4, 4])

Returning the frequency (w) and complex response (h)
Read the help for both of those functions (also available under same name in MATLAB and Octave) for further info.
For MIMO systems the control package has this all worked out as part of the bode function.