#   GNUPLOT v3.6 beta multiplot script file
#
# Second Order System Characteristics
#
#       D**2 + 2*zeta*wn*D + (wn**2)y = (wn**2)*x
#
#   x           input variable
#   y           output variable
#   w           frequency ratio (w/wn)
#   wn          natural frequency
#   wd          damped natural frequency
#   zeta        damping ratio 
#   mag(w)      amplitude response
#   phi(w)      phase response
#   wdwn        damped natural frequency ratio
#   wnt         normalized time
#
#   Plots:
#   Frequency domain    magnitude response
#                           phase response
#
#   Time domain         unit step response
#                       unit impulse response
#
#
# Created by: W. D. Kirby email: wdkirby@ix.netcom.com 
# Date: 1/18/96
# Released to the public domain with no warranty of any kind
#
reset
set style function lines
set size 1.0, 1.0
set origin 0.0, 0.0
set multiplot
set size 0.5,0.5
set origin 0.0,0.5
set grid
unset key
set angles radians
set samples 250
#  Plot Magnitude Response
set title "Second Order System Transfer Function - Magnitude"
mag(w) = -10*log10( (1-w**2)**2 + 4*(zeta*w)**2)
set dummy w
set logscale x
set xlabel "Frequency (w/wn)"
set ylabel "Magnitude (dB)" offset 1,0
set label 1 "Damping =.1,.2,.3,.4,.5,.707,1.0,2.0" at .14,17
set xrange [.1:10]
set yrange [-40:20]
plot \
  zeta=.1,mag(w), \
  zeta=.2,mag(w), \
  zeta=.3,mag(w), \
  zeta=.4,mag(w), \
  zeta=.5,mag(w), \
  zeta=.707,mag(w), \
  zeta=1.0,mag(w), \
  zeta=2.0,mag(w),-6
#  Plot Phase Response
set size 0.5,0.5
set origin 0.0,0.0
set title "Second Order System Transfer Function - Phase"
set label 1 ""
set ylabel "Phase (deg)" offset 1,0
set ytics -180, 30, 0 
set yrange [-180:0]
tmp(w) = (-180/pi)*atan( 2*zeta*w/(1-w**2) )
# Fix for atan function wrap problem
tmp1(w)= w<1?tmp(w):(tmp(w)-180)
phi(w)=zeta==1?(-2*(180/pi)*atan(w)):tmp1(w)
plot \
  zeta=.1,phi(w), \
  zeta=.2,phi(w), \
  zeta=.3,phi(w), \
  zeta=.4,phi(w), \
  zeta=.5,phi(w), \
  zeta=.707,phi(w), \
  zeta=1,phi(w), \
  zeta=2.0,phi(w), \
  -90
#  Plot Step Response
set size 0.5,0.5
set origin 0.5,0.5
set dummy wnt
unset logscale x
set title "Second Order System - Unit Step Response"
set ylabel "Amplitude y(wnt)" offset 1,0 
set xlabel "Normalized Time (wnt)"
set xrange [0:20]
set xtics 0,5,20
set yrange [0:2.0]
set ytics 0, .5, 2.0
set mytics 5
set mxtics 10
wdwn(zeta)=sqrt(1-zeta**2)
shift(zeta) = atan(wdwn(zeta)/zeta)
alpha(zeta)=zeta>1?sqrt(zeta**2-1.0):0
tau1(zeta)=1/(zeta-alpha(zeta))
tau2(zeta)=1/(zeta+alpha(zeta))
c1(zeta)=(zeta + alpha(zeta))/(2*alpha(zeta))
c2(zeta)=c1(zeta)-1
y1(wnt)=zeta==1?1 - exp(-wnt)*(wnt + 1):0
y2(wnt)=zeta<1?(1 - (exp(-zeta*wnt)/wdwn(zeta))*sin(wdwn(zeta)*wnt + shift(zeta))):y1(wnt)
y(wnt)=zeta>1?1-c1(zeta)*exp(-wnt/tau1(zeta))+c2(zeta)*exp(-wnt/tau2(zeta)):y2(wnt)
plot \
  zeta=.1,y(wnt), \
  zeta=.2,y(wnt), \
  zeta=.3,y(wnt), \
  zeta=.4,y(wnt), \
  zeta=.5,y(wnt), \
  zeta=.707,y(wnt), \
  zeta=1,y(wnt), \
  zeta=2,y(wnt)
#
#  Plot Impulse Response
set origin .5,0.
set title "Second Order System - Unit Impulse Response"
y(wnt)=exp(-zeta*wnt) * sin(wdwn(zeta)*wnt) / wdwn(zeta)
set yrange [-1. :1.]
set ytics -1,.5,1.
plot \
  zeta=.1,y(wnt), \
  zeta=.2,y(wnt), \
  zeta=.3,y(wnt), \
  zeta=.4,y(wnt), \
  zeta=.5,y(wnt), \
  zeta=.707,y(wnt), \
  zeta=1,y(wnt), \
  zeta=2,y(wnt)
unset multiplot

Click here for minimal script to generate this plot