Modeling a linear system with Python
As @Matt said, I know this is old. But this came up as my first google hit, so I wanted to edit it.
You can use scipy.signal.lti to model linear, time invariant systems. That gives you lti.bode.
For an impulse response in the form of H(s) = (As^2 + Bs + C)/(Ds^2 + Es + F), you would enter h = scipy.signal.lti([A,B,C],[D,E,F]). To get the bode plot, you would do plot(*h.bode()[:2]).
I would like to simulate/model a closed-loop, linear, time-invariant system (specifically a locked PLL approximation) with python. 
Each sub-block within the model has a known transfer function which is given in terms of complex frequency H(s) = K / ( s * tau + 1 ).  Using the model, I would like to see how the system response as well as the noise response is affected as parameters (e.g. the VCO gain) are changed.  This would involve using Bode plots and root-locus plots.  
What Python modules should I seek out to get the job done?