diff --git a/transient-response/examples/first-order-dynamic-systems/impulse.py b/transient-response/examples/first-order-dynamic-systems/impulse.py
new file mode 100644
index 0000000..947a9b7
--- /dev/null
+++ b/transient-response/examples/first-order-dynamic-systems/impulse.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+System Dynamics and Control                                 August 25, 2021
+Prof. M. Diaz-Maldonado
+
+
+Synopsis:
+Solves for the transient response of a first-order system when subjected
+to a step and impulse responses:
+
+                    y' + k * y = b * u(t),
+
+where k is the rate constant, b is the ``forcing'' constant, and u(t)
+is either the unit step or the unit impulse input functions.
+
+
+Copyright (c) 2021 Misael Diaz-Maldonado
+This file is released under the GNU General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+
+References:
+[0] A Gilat and V Subramanian, Numerical Methods for Engineers and
+    Scientists: An Introduction with Applications using MATLAB
+[1] R Johansson, Numerical Python: Scientific Computing and Data
+    Science Applications with NumPy, SciPy, and Matplotlib, 2nd edition
+[2] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+"""
+
+
+import numpy as np
+from numpy import exp, linspace
+from scipy.integrate import solve_ivp
+import matplotlib as mpl
+mpl.use("qt5agg")
+from matplotlib import pyplot as plt
+
+
+# initial value, and rate and forcing constants, respectively
+yi, k, b = (0.0, 1.0, 1.0)
+# defines lambdas for the analytic step and impulse-responses
+step    = lambda t: (yi - b / k) * exp(-k * t) + b / k 
+impulse = lambda t: -k * (yi - b / k) * exp(-k * t)
+# defines the right-hand side RHS of the ODE: dy/dt = f(t, y) as a lambda
+odefun = lambda t, y: (b - k * y)
+
+
+""" solves the transient response via the 4th-order Runge-Kutta Method """
+ti, tf = (0.0, 12.0)            # initial time and final integration times
+tspan  = (ti, tf)               # time span
+t = linspace(ti, tf, 256)
+odesol = solve_ivp(odefun, tspan, [yi], method="RK45")
+# unpacks the numerical solution
+t_scipy, y_scipy = (odesol.t, odesol.y[0, :])
+#y_scipy = y_scipy[0, :]
+
+
+# plots step response
+plt.close("all")
+plt.ion()
+fig, ax = plt.subplots()
+
+ax.plot(t, step(t), color="black", linewidth=2.0,
+        label="analytic solution")
+ax.plot(t_scipy, y_scipy, color="blue", marker="v", linestyle="",
+        label="scipy's 4th-order Runge-Kutta method")
+
+ax.grid()
+ax.legend()
+ax.set_xlabel("time, t")
+ax.set_ylabel("dynamic response, y(t)")
+ax.set_title("Step Response of a First Order Dynamic System")
+
+
+
+
+# plots impulse response
+fig, ax = plt.subplots()
+
+ax.plot(t, impulse(t), color="black", linewidth=2.0,
+        label="analytic solution")
+ax.plot(t_scipy, odefun(t_scipy, y_scipy), color="blue", marker="v",
+        linestyle="", label="scipy's 4th-order Runge-Kutta method")
+
+ax.grid()
+ax.legend()
+ax.set_xlabel("time, t")
+ax.set_ylabel("dynamic response, y(t)")
+ax.set_title("Impulse Response of a First Order Dynamic System")
diff --git a/transient-response/examples/higher-order-systems/step-response-underdamped-third-order-system-sympy.py b/transient-response/examples/higher-order-systems/step-response-underdamped-third-order-system-sympy.py
new file mode 100644
index 0000000..9006180
--- /dev/null
+++ b/transient-response/examples/higher-order-systems/step-response-underdamped-third-order-system-sympy.py
@@ -0,0 +1,254 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+System Dynamics and Control                                 August 29, 2021
+Prof. M Diaz-Maldonado
+
+
+Synopsis:
+Solves the Initial Value Problem IVP that models the underdamped response,
+y(t), of an underdamped third-order dynamic system:
+
+y''' + a2 * y'' + a1 * y' + a0 * y = b0 * u(t), y'' = y' = y = 0 at t = 0,
+
+where a0, a1, a2, and b0 are the ``standard'' coefficients of the Ordinary
+Differential Equation ODE, and u(t) is the step response of magnitude h.
+
+
+Copyright (c) 2021 Misael Diaz-Maldonado
+This file is released under the GNU General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+
+References:
+[0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+[1] R Johansson, Numerical Python: Scientific Computing and Data
+    Science Applications with NumPy, SciPy, and Matplotlib, 2nd edition
+[2] (linsolve: docs.sympy.org/latest/modules/solvers/solveset.html)
+[3] (sequence unpacking: stackoverflow.com/questions/38418205/
+     get-a-value-from-solution-set-returned-as-finiteset-by-sympy)
+"""
+
+
+import sympy
+from sympy import diff as D
+from numpy import pi, sqrt
+from numpy import empty, linspace
+from scipy.optimize import bisect
+from scipy.integrate import solve_ivp
+import matplotlib as mpl
+mpl.use("qt5agg")
+from matplotlib import pyplot as plt
+
+
+def fstep(R, Alpha, Beta, B0, H):
+    """
+    Synopsis:
+    Returns lambdas to obtain the step and impulse responses of the
+    third-order system given the values of the real, complex roots, and
+    input step parameters.
+
+    inputs:
+    R               real root
+    Alpha, Beta     real and imaginary components of the complex roots
+    B0, H           forcing constant and step magnitude
+
+    outputs:
+    step, impulse   lambdas
+    """
+
+    # time
+    t = sympy.Symbol('t')
+    # characteristic roots for the underdamped dynamic system
+    r = sympy.Symbol('r')                           # real
+    alpha, beta = sympy.symbols( ('alpha', 'beta') )# complex conjugates
+
+
+    # expresses the standard coefficients of the ODE in terms of the roots
+    a0 = -r * (alpha**2 + beta**2)
+    a1 = 2 * r * alpha  + (alpha**2 + beta**2)
+    a2 = -(r + 2 * alpha)
+    b0 = sympy.Symbol('b0')
+
+
+    # defines the homogeneous solution, y
+    y1 = sympy.exp(r * t)
+    y2 = sympy.exp(alpha * t) * sympy.sin(beta * t)
+    y3 = sympy.exp(alpha * t) * sympy.cos(beta * t)
+
+    y = y1 + y2 + y3
+
+    # defines the particular solution, yp
+    h = sympy.Symbol('h')
+    u = h                   # u(t), step
+    yp = b0 * h / a0
+
+
+    """ solves the Initial Value Problem IVP """
+    # matrix
+    Mat = sympy.Matrix([
+        [y1.subs({t: 0}), y2.subs({t: 0}), y3.subs({t: 0})],
+        [D(y1,t).subs({t: 0}), D(y2,t).subs({t: 0}), D(y3,t).subs({t: 0})],
+        [D(y1,t,t).subs({t: 0}), D(y2,t,t).subs({t: 0}),
+            D(y3,t,t).subs({t: 0})],
+    ])
+    # coefficient vector
+    vec = sympy.Matrix( [-yp.subs({t: 0}), 0, 0] )
+
+    # obtains the undetermined coefficient of the general solution
+    A, B, C = sympy.symbols( ('A', 'B', 'C') )
+    """Note: trailing comma enables sequence unpacking of the Finite Set"""
+    sol, = sympy.linsolve( (Mat, vec), A, B, C )
+
+    # defines the general solution (transient response `y(t)')
+    A, B, C = sol
+    y = A * y1 + B * y2 + C * y3 + yp
+
+    step = sympy.lambdify(
+        t, y.subs({r: R, alpha: Alpha, beta: Beta, b0: B0, h: H}), 'numpy'
+    )
+
+    impulse = sympy.lambdify(
+        t, D(y, t).subs({r: R, alpha: Alpha, beta: Beta, b0: B0, h: H}), 'numpy'
+    )
+
+    return (step, impulse)
+
+
+def perf(r, alpha, beta, step, impulse):
+    """
+    Synopsis:
+    Obtains the performance characteristics of the third-order underdamped
+    dynamic system.
+    """
+
+    omega = sqrt(alpha**2 + beta**2)                    # natural frequency
+    zeta  = -alpha / omega                              # damping ratio
+    
+    ts = max(-4 / r,  -4 / alpha)                       # settling time
+    P = Period = 2 * pi / ( omega * sqrt(1 - zeta**2) ) # period
+    N = Ncycles = ts / Period                           # number of cycles
+
+    # Note: See comments at the end of the source regarding the peak time
+    tp = bisect(impulse, 0.15 * P, 0.85 * P)            # peak time
+    y_max = step(tp)                                    # peak value
+    
+    A0 = ( -R * (Alpha**2 + Beta**2) )                  # A0 = a0
+    y_ss = B0 * H / A0                                  # steady-state
+
+
+    """ displays the performance characteristics of the dynamic system """
+    performance = (
+        f"\n\n"
+        f"Damping Ratio:                   {zeta}\n"
+        f"Natural Frequency:               {omega}\n"
+        f"peak time and value:             {tp}, {y_max}\n"
+        f"Settling time:                   {ts}\n"
+        f"steady-state value:              {y_ss}\n" 
+        f"Maximum Overshoot:               {y_max / y_ss - 1}\n"
+        f"Oscillation Period:              {Period}\n"
+        f"Number of Oscillation Cycles:    {Ncycles}\n\n"
+    )
+
+    print(performance)
+
+    return (omega, zeta, ts, tp, P, N, y_max)
+
+
+
+
+""" plots the analytic and numeric transient responses """
+time = linspace(0, 25, 256)
+
+R, Alpha, Beta, B0, H = (-1, -0.25, 1, 1, 1)
+
+step, impulse = fstep(R, Alpha, Beta, B0, H)
+
+
+def odefun(t, y):
+    """
+    Synopsis:
+    The state-variable equations for the third-order underdamped system.
+    Defines the third-order system as an equivalent system of first-order
+    ODEs to obtain the transient response y(t) numerically.
+    """
+
+    # defines the ``standard'' coefficients of the third-order ODE
+    A0 = -R * (Alpha**2 + Beta**2)
+    A1 = 2 * R * Alpha  + (Alpha**2 + Beta**2)
+    A2 = -(R + 2 * Alpha)
+
+    f = empty(3)
+
+    f[0] = y[1]
+    f[1] = y[2]
+    f[2] = -A2 * y[2] - A1 * y[1] - A0 * y[0] + B0 * H
+
+    return f
+
+
+# solves the third-order ODE numerically for verification
+tspan = [0, 25]
+""" initial values: y(0) = 0, y'(0) = 0, y''(0) = 0 """
+yi = [0, 0, 0]
+# applies the fourth-order Runge-Kutta method
+odesol = solve_ivp(odefun, tspan, yi, method="RK45")
+t, y = odesol.t, odesol.y
+y_step = y[0, :]    # step response
+y_impulse = y[1, :] # impulse response
+
+
+plt.close("all")
+plt.ion()
+
+fig, ax = plt.subplots()
+ax.plot(time, step(time), linewidth = 2, color = "black",
+        label = "analytic")
+ax.plot(t, y_step, linestyle = "", marker = "d", color = "orange",
+        label = "numeric")
+ax.set_xlabel("time, t")
+ax.set_ylabel("transient response, y(t)")
+ax.set_title("Step Response of an Underdamped Third-Order Dynamic System")
+ax.legend()
+ax.grid()
+
+
+# displays performance characteristics on the console
+perf(R, Alpha, Beta, step, impulse)
+
+
+""" plots the impulse response """
+fig, ax = plt.subplots()
+ax.plot(time, impulse(time), linewidth = 2, color = "black",
+        label = "analytic")
+ax.plot(t, y_impulse, linestyle = "", marker = "d", color = "orange",
+        label = "numeric")
+ax.set_xlabel("time, t")
+ax.set_ylabel("transient response, y(t)")
+ax.set_title(
+    "Impulse Response of an Underdamped Third-Order Dynamic System"
+)
+ax.legend()
+ax.grid()
+
+
+
+"""
+Comments on the Determination of the Peak Time:
+If the dynamics of the third-order system is dominated by the complex
+roots (rather than the real root) the peak time can be observed in the
+first oscillation cycle as for second-order systems. Only then one can
+solve for the peak time by applying a numerical technique. Note that
+solving for the peak time entails solving the nonlinear equation (in time)
+y'(t) = 0. Here we applied the bisection method, though the bracketing
+interval may need adjusting for other third-order underdamped systems
+(different set of parameters). You have the option of applying other
+numerical techniques such as Newton-Raphson.
+
+If the real root dominates, the system will exhibit weak oscillations
+and its transient response will resemble that of a first-order system.
+And its maximum value of y(t) is the steady-state value. This means
+that the bisection method will complain that there's no solution in
+the given interval. You may want to comment out that part of the code.
+"""
diff --git a/transient-response/examples/second-order-dynamic-systems/impulse.py b/transient-response/examples/second-order-dynamic-systems/impulse.py
new file mode 100644
index 0000000..05d3204
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/impulse.py
@@ -0,0 +1,172 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+System Dynamics and Control                                 August 29, 2021
+Prof. M Diaz-Maldonado
+
+
+Synopsis:
+Solves for the impulse response of an underdamped mechanical system
+described by the second-order Ordinary Differential Equation ODE:
+
+                m * y'' + b * y' + k * y = f * u(t),
+
+where m is the mass, b is the damper friction, k is the spring stiffness,
+f is the forcing constant, and u(t) is the impulse of magnitude P.
+
+The code produces a plot qualitatively similar to that shown in figure 7.22
+of Kluever's textbook.
+
+
+Copyright (c) 2021 Misael Diaz-Maldonado
+This file is released under the GNU General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+
+References:
+[0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+[1] R Johansson, Numerical Python: Scientific Computing and Data
+    Science Applications with NumPy, SciPy, and Matplotlib, 2nd edition
+"""
+
+
+from numpy import linspace
+from numpy import pi, cos, sin, exp, sqrt, arctan
+import matplotlib as mpl
+mpl.use("qt5agg")
+mpl.rcParams['text.usetex'] = True
+from matplotlib import pyplot as plt
+
+
+def impulse(prms, t):
+    """
+    Synopsis:
+    Computes the impulse-response y(t) of an underdamped second-order
+    dynamic system given the system mass, friction and stiffness constants,
+    and the forcing constant and magnitude P of the impulse u(t).
+    """
+
+    # unpacks mechanical system params (inertia, friction, stiffness, etc.)
+    m, b, k, f, P = prms
+    # defines the standard coefficients of the second-order ODE
+    a0, a1, b0 = (k / m, b / m, f / m)
+    # determines the natural frequency and damping ratio, respectively
+    omega, zeta = (  sqrt(a0),  a1 / ( 2 * sqrt(a0) )  )
+    # complains if the system is not underdamped
+    is_underdamped(zeta)
+    # obtains the real and imaginary components of the characteristic roots
+    alpha, beta = ( -zeta * omega, omega * sqrt(1 - zeta**2) )
+    
+    # computes the system response y(t) to a unit-impulse input
+    y = beta * (1 + (alpha / beta)**2) * exp(alpha * t) * sin(beta * t)
+    # rescales to obtain the system response to an impulse of magnitude P
+    y *= b0 * P / omega**2
+
+    return y
+
+
+def perf(prms):
+    """
+    Synopsis:
+    Obtains the performance characteristics of a second-order underdamped
+    dynamic system.
+    """
+
+    # calculates the underdamped mechanical system characteristics
+    a0, a1, b0, omega, zeta = underdamped(prms)
+    alpha, beta = ( -zeta * omega, omega * sqrt(1 - zeta**2) )
+  
+    """ peak time and value """
+    tp = arctan(-beta / alpha) / ( omega * sqrt(1 - zeta**2) )
+    y_max = impulse(prms, tp)
+
+    y_ss = 0.0                                  # steady-state value
+    ts = 4 / (zeta * omega)                     # settling time
+    T = 2 * pi / ( omega * sqrt(1 - zeta**2) )  # period of oscillation
+    N = 2 * sqrt(1 - zeta**2) / (pi * zeta)     # number of cycles
+
+    return (tp, ts, T, N, y_max, y_ss)
+
+
+def underdamped(prms):
+    """
+    Synopsis:
+    Returns the basic characteristics of an underdamped mechanical system
+    given the mass/moment of inertia and friction and stiffness constants.
+    """
+    # unpacks mechanical system parameters (inertia, friction, stiffness)
+    m, b, k, f, P = prms
+    # computes the standard coefficients of the second-order ODE
+    a0, a1, b0 = (k / m, b / m, f / m)
+    # determines the natural frequency and damping ratio, respectively
+    omega, zeta = (  sqrt(a0),  a1 / ( 2 * sqrt(a0) )  )
+    # complains if not underdamped
+    is_underdamped(zeta)
+    return (a0, a1, b0, omega, zeta)
+
+
+def is_underdamped(zeta):
+    # complains if the mechanical system is not underdamped
+
+    errMSG = (f"Not an underdamped mechanical system\n" +
+        f"              the damping ratio (zeta = {zeta}) is not " +
+        f"less than one\n"
+    )
+
+    if zeta >= 1.0:
+        raise RuntimeError(errMSG)
+    return
+
+
+""" mechanical system parameters """
+# moment of inertia, friction, and stiffness
+prms = m, b, k, f, P = (0.2, 1.6,  65.0, 1.0, 2.5)
+#prms = m, b, k, f, P = (1.0, 8.0, 325.0, 5.0, 2.5) # equivalent set
+
+
+""" simulation parameters """
+# time t, seconds
+t = linspace(0, 2, 512)
+
+
+""" displays performance characteristics """
+tp, ts, Period, Ncycles, y_max, y_ss = perf(prms)
+performance = (
+    f"\n\n"
+    f"peak time and value:             {tp}, {y_max}\n"
+    f"settling time:                   {ts}\n"
+    f"steady-state value:              {y_ss}\n"
+    f"Oscillation Period:              {Period}\n"
+    f"Number of Oscillation Cycles:    {Ncycles}\n\n"
+)
+print(performance)
+
+
+""" transient response visualization """
+plt.close("all")
+plt.ion()
+fig, ax = plt.subplots()
+ax.plot(t, impulse(prms, t), color = "black", linestyle = "-",
+        linewidth = 2)
+ax.set_xlabel(r'time, $t$, seconds')
+ax.set_ylabel(r'transient response, $y(t)$')
+ax.set_title(
+    'Impulse Response of an Underdamped Second Order Dynamic System'
+)
+ax.grid()
+
+fig.savefig("impulse-response-underdamped-2nd-order-system.png", dpi=300)
+
+
+
+"""
+Comments:
+The performance equations for the settling time, period of oscillation,
+and number of oscillation cycles are the same as those for the step
+response. The peak time, however, occurs at a different time which
+can be determined analytically.
+
+Obtaining the first-derivative of the impulse response, y'(t), and finding
+the time where y'(t) is equal to zero yields the sought expression.
+"""
diff --git a/transient-response/examples/second-order-dynamic-systems/step.m b/transient-response/examples/second-order-dynamic-systems/step.m
new file mode 100644
index 0000000..8634257
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/step.m
@@ -0,0 +1,61 @@
+% System Dynamics and Control                             December 27, 2019
+% ME 3030-21 WI19
+% Prof. M Diaz-Maldonado
+%
+% Synopsis:
+% Reproduces the step response of example 7.8 from Kluever's textbook.
+%
+%
+% Copyright (c) 2021 Misael Diaz-Maldonado
+% This file is released under the GNU General Public License as published
+% by the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+%
+% References:
+% [0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+% [1] A Gilat, MATLAB: An Introduction with Applications, 6th edition
+
+clear
+close all
+clc
+
+% mechanical system parameters
+J = 0.2;			% moment of inertia
+b = 1.6;			% friction coefficient
+k = 65.;			% stiffness
+
+% magnitude of step input
+A = 2.5;
+
+% standard coefficients of the second-order system
+a1 = b / J;	a0 = k / J;	b0 = 1 / J;
+
+
+% defines the second-order constants which specify a second-order system:
+
+Wn = sqrt(a0);			% natural frequency
+DR = a1 / (2 * Wn);		% damping ratio
+
+% real and imaginary parts of the roots of the characteristic equation
+alpha = -DR * Wn;		% real
+beta  = Wn * sqrt(1 - DR^2);	% imaginary
+
+
+% calculates the transient response, y(t), Equation (7.72)
+t = linspace(0, 2, 512);
+y = (  1 - exp(alpha * t) .* ( cos(beta*t) - alpha/beta * sin(beta*t) )  );
+y = (b0 * A / Wn^2) * y;
+plot(t, y, '-k', 'linewidth', 2)
+hold on
+
+
+
+% specifies the figure labels, the title, and the legend
+xlabel('time, t, seconds')
+ylabel('transient response, y(t), radians')
+title('Step Response of an Underdamped Mechanical System')
+grid on
+
+% exports figure to PNG format with 600 DPI
+print('step-response-mechanical-system.png','-r600','-dpng')
diff --git a/transient-response/examples/second-order-dynamic-systems/step.py b/transient-response/examples/second-order-dynamic-systems/step.py
new file mode 100644
index 0000000..054fbdc
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/step.py
@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+System Dynamics and Control                                 August 27, 2021
+Prof. M Diaz-Maldonado
+
+
+Synopsis:
+Solves for the step response of an underdamped mechanical system described
+by the second-order Ordinary Differential Equation ODE:
+
+                m * y'' + b * y' + k * y = f * u(t),
+
+where m is the mass, b is the damper friction, k is the spring stiffness,
+f is the forcing constant, and u(t) = H is the step-input of magnitude H.
+
+The code produces a plot similar to that shown in example 7.8 of Kluever's
+textbook.
+
+
+Copyright (c) 2021 Misael Diaz-Maldonado
+This file is released under the GNU General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+
+References:
+[0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+[1] R Johansson, Numerical Python: Scientific Computing and Data
+    Science Applications with NumPy, SciPy, and Matplotlib, 2nd edition
+"""
+
+
+from numpy import linspace
+from numpy import pi, cos, sin, exp, sqrt
+import matplotlib as mpl
+mpl.use("qt5agg")
+mpl.rcParams['text.usetex'] = True
+from matplotlib import pyplot as plt
+
+
+def step(prms, t):
+    """
+    Synopsis:
+    Computes the step-response y(t) of an underdamped second-order dynamic
+    system given the system mass, friction and stiffness constants, and the
+    forcing constant and magnitude H of the step-input u(t).
+    """
+
+    # unpacks mechanical system params (inertia, friction, stiffness, etc.)
+    m, b, k, f, H = prms
+    # defines the standard coefficients of the second-order ODE
+    a0, a1, b0 = (k / m, b / m, f / m)
+    # determines the natural frequency and damping ratio, respectively
+    omega, zeta = (  sqrt(a0),  a1 / ( 2 * sqrt(a0) )  )
+    # complains if the system is not underdamped
+    is_underdamped(zeta)
+    # obtains the real and imaginary components of the characteristic roots
+    alpha, beta = ( -zeta * omega, omega * sqrt(1 - zeta**2) )
+    
+    # computes the system response y(t) to a unit-step input
+    y = 1 - exp(alpha * t) * ( cos(beta * t) - alpha / beta *
+        sin(beta * t) )
+    # rescales to obtain the system response to a step of magnitude H
+    y *= b0 * H / omega**2
+
+    return y
+
+
+def perf(prms):
+    """
+    Synopsis:
+    Obtains the performance characteristics of a second-order underdamped
+    dynamic system.
+    """
+
+    # calculates the underdamped mechanical system characteristics
+    a0, a1, b0, omega, zeta = underdamped(prms)
+  
+    tp = pi / ( omega * sqrt(1 - zeta**2) )     # peak time
+    y_max = step(prms, tp)                      # peak value
+    y_ss = b0 * H / omega**2                    # steady-state value
+    ts = 4 / (zeta * omega)                     # settling time
+    P = 2 * pi / ( omega * sqrt(1 - zeta**2) )  # period of oscillation
+    N = 2 * sqrt(1 - zeta**2) / (pi * zeta)     # number of cycles
+    Mos = exp(-pi * zeta / sqrt(1 - zeta**2) )  # maximum overshoot
+
+    return (tp, ts, P, N, Mos, y_max, y_ss)
+
+
+def underdamped(prms):
+    """
+    Synopsis:
+    Returns the basic characteristics of an underdamped mechanical system
+    given the mass/moment of inertia and friction and stiffness constants.
+    """
+    # unpacks mechanical system parameters (inertia, friction, stiffness)
+    m, b, k, f, H = prms
+    # computes the standard coefficients of the second-order ODE
+    a0, a1, b0 = (k / m, b / m, f / m)
+    # determines the natural frequency and damping ratio, respectively
+    omega, zeta = (  sqrt(a0),  a1 / ( 2 * sqrt(a0) )  )
+    # complains if not underdamped
+    is_underdamped(zeta)
+    return (a0, a1, b0, omega, zeta)
+
+
+def is_underdamped(zeta):
+    # complains if the mechanical system is not underdamped
+
+    errMSG = (f"Not an underdamped mechanical system\n" +
+        f"              the damping ratio (zeta = {zeta}) is not " +
+        f"less than one\n"
+    )
+
+    if zeta >= 1.0:
+        raise RuntimeError(errMSG)
+    return
+
+
+""" mechanical system parameters """
+# moment of inertia, friction, and stiffness
+prms = m, b, k, f, H = (0.2, 1.6,  65.0, 1.0, 2.5)
+#prms = m, b, k, f, H = (1.0, 8.0, 325.0, 5.0, 2.5) # equivalent set
+
+
+""" simulation parameters """
+# time t, seconds
+t = linspace(0, 2, 256)
+
+
+""" displays performance characteristics """
+tp, ts, Period, Ncycles, Mos, y_max, y_ss = perf(prms)
+performance = (
+    f"\n\n"
+    f"peak time and value:             {tp}, {y_max}\n"
+    f"settling time:                   {ts}\n"
+    f"steady-state value:              {y_ss}\n"
+    f"Maximum Overshoot:               {Mos}\n"
+    f"Oscillation Period:              {Period}\n"
+    f"Number of Oscillation Cycles:    {Ncycles}\n\n"
+)
+print(performance)
+
+
+""" transient response visualization """
+plt.close("all")
+plt.ion()
+fig, ax = plt.subplots()
+ax.plot(t, step(prms, t), color = "black", linestyle = "-",
+        linewidth = 2)
+ax.set_xlabel(r'time, $t$, seconds')
+ax.set_ylabel(r'transient response, $y(t)$')
+ax.set_title('Step Response of an Underdamped Second Order Dynamic System')
+ax.grid()
+
+fig.savefig("figure-7.21.png", dpi=300)
+
+
+
+"""
+Comments:
+The ordinary differential equation that describes the dynamics of the
+mechanical system is:
+                    y'' + a1 * y' + a0 * y = b0 * u(t),
+where a1, a0, and b0 are the standard coefficients, u(t) is the input, and
+y = y(t) is the system response. The constants are defined for a mechanical
+system in the step() and underdamped() functions.
+"""
diff --git a/transient-response/examples/second-order-dynamic-systems/unit-step-response.py b/transient-response/examples/second-order-dynamic-systems/unit-step-response.py
new file mode 100644
index 0000000..c2e8d54
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/unit-step-response.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+System Dynamics and Control                                 August 26, 2021
+Prof. M Diaz-Maldonado
+
+
+Synopsis:
+Obtains the transient response of an underdamped second-order system to a
+unit-step input. (Reproduces Figure 7.19 of Kluever's textbook.)
+
+
+Copyright (c) 2021 Misael Diaz-Maldonado
+This file is released under the GNU General Public License as published
+by the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+
+References:
+[0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+[1] R Johansson, Numerical Python: Scientific Computing and Data
+    Science Applications with NumPy, SciPy, and Matplotlib, 2nd edition
+"""
+
+from numpy import linspace
+from numpy import cos, sin, exp, sqrt
+import matplotlib as mpl
+mpl.use("qt5agg")
+mpl.rcParams['text.usetex'] = True      # uses TeX for Math Symbols
+from matplotlib import pyplot as plt
+
+
+def step(zeta, omega, t):
+    """
+    Synopsis:
+    Computes the step response y(t) of an underdamped second-order system
+    given the damping ratio zeta, the natural frequency omega, and time t.
+    """
+    # obtains the real and imaginary components of the characteristic roots
+    alpha, beta = ( -zeta * omega, omega * sqrt(1 - zeta**2) )
+
+    y = 1 - exp(alpha * t) * ( cos(beta * t) - alpha / beta * 
+        sin(beta * t) )
+    
+    return y
+
+
+""" defines system parameters """
+zetas = [0.2, 0.4]          # damping ratios
+omega = 1.0                 # natural frequency, radians / second
+t = linspace(0, 25, 500)    # time, seconds
+
+
+""" plots the response for varying damping ratios on the same figure """
+# defines linestyles and labels
+colors = ["black", "red"]
+linestyles = ["-", "--"]
+labels = [r"$\zeta = 0.2$", r"$\zeta = 0.4$"]
+
+plt.close("all")
+plt.ion()
+fig, ax = plt.subplots()
+
+for n, zeta in enumerate(zetas):
+    ax.plot(t, step(zeta, omega, t), color = colors[n],
+            linestyle = linestyles[n], label = labels[n])
+
+ax.set_xlabel('time, t, seconds')
+ax.set_ylabel('transient response, y(t)')
+ax.set_title('Step Response of an Underdamped Second Order Dynamic System')
+ax.legend()
+
+fig.savefig("unit-step-response-second-order-dynamic-system.png", dpi=300)
+
+
+
+"""
+Final Remarks:
+The transient response obtained here applies when the dynamics of the
+mechanical system is described by equation (7.69):
+                y'' + a1 * y' + a0 * y = a0 * u(t).
+The equation describes the dynamics of a mass interconnected with a
+movable wall via a spring. The damper is interconnected with a fixed wall
+on the opposite side of the mass. The input step u(t) corresponds to the
+wall displacement, which can be verified by cheking the units of each
+term.
+"""
diff --git a/transient-response/examples/second-order-dynamic-systems/unit_impulse_response.m b/transient-response/examples/second-order-dynamic-systems/unit_impulse_response.m
new file mode 100644
index 0000000..c532b9f
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/unit_impulse_response.m
@@ -0,0 +1,51 @@
+% System Dynamics and Control                       	   January 22, 2020
+% ME 3030-21 WI19
+% Prof. M Diaz-Maldonado
+%
+% Synopsis:
+% Plots the response of an underdamped, second-order, dynamic system to a
+% unit-impulse (P = 1) input.
+%
+%
+% Copyright (c) 2021 Misael Diaz-Maldonado
+% This file is released under the GNU General Public License as published
+% by the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+%
+% References:
+% [0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+% [1] A Gilat, MATLAB: An Introduction with Applications, 6th edition
+
+clear
+close all
+clc
+
+% defines the second-order constants which specify a second-order system:
+DR = 0.15;		% damping ratio
+Wn = 1.00;		% natural frequency
+
+% real and imaginary parts of the roots of the characteristic equation
+alpha = -DR * Wn;
+beta  = Wn * sqrt(1 - DR^2);
+
+% calculates the transient response, y(t)
+% Note: b0 = Wn^2 and P=1, Kluever's Equation (7.66)
+t = linspace(0, 30, 1024);
+y = beta * (1 + (alpha / beta)^2) * exp(alpha * t) .* sin(beta * t);
+plot(t, y, '-k', 'linewidth', 2); hold on
+
+% plots the impulse response for a higher damping ratio on the same figure
+DR = 0.3;
+alpha = -DR * Wn;
+beta  = Wn * sqrt(1 - DR^2);
+y = beta * (1 + (alpha / beta)^2) * exp(alpha * t) .* sin(beta * t);
+plot(t, y, '-r', 'linewidth', 2)
+
+% specifies the figure labels, the title, and the legend
+xlabel('time, t, sec')
+ylabel('transient response, y(t)')
+title('Impulse Response of an Underdamped Second Order Dynamic System')
+grid on
+legend('\zeta = 0.15', '\zeta = 0.30')	% damping ratios
+print('impulse-response-second-order-system.png','-r600','-dpng')
diff --git a/transient-response/examples/second-order-dynamic-systems/unit_pulse_response.m b/transient-response/examples/second-order-dynamic-systems/unit_pulse_response.m
new file mode 100644
index 0000000..effc4ea
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/unit_pulse_response.m
@@ -0,0 +1,52 @@
+% System Dynamics and Control                             December 27, 2019
+% ME 3030-21 WI19
+% Prof. M Diaz-Maldonado
+%
+% Synopsis:
+% Obtains the pulse response of an underdamped second-order system by
+% applying the superposition property. The weight of the pulse A = 1.
+% It's seen that as T -> 0 the response converges to the impulse response.
+%
+%
+% Copyright (c) 2021 Misael Diaz-Maldonado
+% This file is released under the GNU General Public License as published
+% by the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+%
+% References:
+% [0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+% [1] A Gilat, MATLAB: An Introduction with Applications, 6th edition
+
+clear
+close all
+clc
+
+% defines the second-order constants which specify a second-order system:
+DR = 0.2;		% damping ratio
+Wn = 1.0;		% natural frequency
+T = 0.5;		% pulse duration (sec)
+P = 1 / T;		% pulse intensity (note the pulse weight is A=1)
+
+% real and imaginary parts of the roots of the characteristic equation
+alpha = -DR * Wn;
+beta  = Wn * sqrt(1 - DR^2);
+
+
+% calculates the transient response, y(t), Equation (7.72)
+t = linspace(0, 30, 1024);
+y = @(t) 1 - exp(alpha * t) .* ( cos(beta*t) - alpha/beta * sin(beta*t) );
+Y = P * ( y(t) - y(t - T) .* (t > T) );
+
+yim = beta * (1 + (alpha / beta)^2) * exp(alpha * t) .* sin(beta * t);
+plot(t, Y, '-k', 'linewidth', 2);    hold on	% pulse
+plot(t, yim, '--', 'linewidth', 2)		% impulse
+legend('pulse', 'impulse')
+
+ylim([-0.5, 1.0])
+xlabel('time, t, sec')
+ylabel('transient response, y(t)')
+title('Pulse Response of an Underdamped Second Order Dynamic System')
+grid on
+
+print('pulse-response-second-order-system.png','-r600','-dpng')
diff --git a/transient-response/examples/second-order-dynamic-systems/unit_step_response.m b/transient-response/examples/second-order-dynamic-systems/unit_step_response.m
new file mode 100644
index 0000000..05ba5e0
--- /dev/null
+++ b/transient-response/examples/second-order-dynamic-systems/unit_step_response.m
@@ -0,0 +1,62 @@
+% System Dynamics and Control		                  December 27, 2019
+% ME 3030-21 WI19
+% Prof. M Diaz-Maldonado
+%
+% Synopsis:
+% Obtains the unit-step response of an underdamped second-order system:
+%
+% 	y'' + a1 * y' + a0 * y = b0 * u(t),	y(0) = y'(0) = 0,
+%
+% where a1 = 2 * zeta * omega, a0 = b0 = omega^2, zeta is the damping
+% ratio, omega is the natural frequency, and u(t) is the unit step
+% function (u(t) = 1, for t > 0).
+%
+% (The code reproduces Figure 7.19 of Kluever's textbook.)
+%
+%
+% Copyright (c) 2021 Misael Diaz-Maldonado
+% This file is released under the GNU General Public License as published
+% by the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+%
+% References:
+% [0] CA Kluever, Dynamic Systems: Modeling, Simulation, and Control
+% [1] A Gilat, MATLAB: An Introduction with Applications, 6th edition
+
+clear
+close all
+clc
+
+% defines the second-order constants which specify a second-order system:
+DR = 0.2;			% damping ratio
+Wn = 1.0;			% natural frequency
+
+% real and imaginary parts of the roots of the characteristic equation
+alpha = -DR * Wn;		% real
+beta  = Wn * sqrt(1 - DR^2);	% imaginary
+
+% calculates the transient response, y(t), Equation (7.72)
+t = linspace(0, 25, 512);
+y = 1 - exp(alpha * t) .* ( cos(beta * t) - alpha / beta * sin(beta * t) );
+plot(t, y, '-k', 'linewidth', 2)
+hold on
+
+
+% plots the step response for a higher damping ratio on the same figure
+DR = 0.4;
+alpha = -DR * Wn;
+beta  = Wn * sqrt(1 - DR^2);
+y = 1 - exp(alpha * t) .* ( cos(beta * t) - alpha / beta * sin(beta * t) );
+plot(t, y, '-r', 'linewidth', 2)
+
+
+% specifies the figure labels, the title, and the legend
+xlabel('time, t, sec')
+ylabel('transient response, y(t)')
+title('Unit-step Response of an Underdamped Second Order Dynamic System')
+legend('\zeta = 0.2', '\zeta = 0.4')
+grid on
+
+% exports figure to PNG format with 600 DPI
+print('unit-step-response-second-order-system.png','-r600','-dpng')
diff --git a/transient-response/examples/second-order-response-transient-response/ch7_example_8_step_response_underdamped_second_order_system.m b/transient-response/examples/second-order-response-transient-response/ch7_example_8_step_response_underdamped_second_order_system.m
deleted file mode 100644
index 3a3f249..0000000
--- a/transient-response/examples/second-order-response-transient-response/ch7_example_8_step_response_underdamped_second_order_system.m
+++ /dev/null
@@ -1,55 +0,0 @@
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Course: System Dynamics and Control                       ME 3030-21 WI19
-% Author: Prof. M Diaz-Maldonado                          Date: Dec/27/2019
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Synopsis:
-% Reproduces the response of example 7.8 from Kluever's textbook.
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-clear;
-close all;
-clc;
-
-% mechanical system parameters
-% moment of inertia
-J = 0.2;
-% friction coefficient
-b = 1.6;
-% stiffness
-k = 65;
-
-% magnitude of step input
-A = 2.5; 
-
-% standard coefficients of the second-order system
-a1 = b/J; a0 = k/J; b0 = 1/J;
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% defines the second-order constants which specify a second-order system:
-
-% natural frequency
-Wn = sqrt(a0);
-% damping ratio
-DR = a1 / (2*Wn);
-
-% real and imaginary parts of the roots of the characteristic equation
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% calculates the transient response, y(t), Equation (7.72)
-t = linspace(0,2,500);
-y = b0 * A / Wn^2 *...
-    (1 - exp(alpha * t) .* (cos(beta*t) - alpha/beta * sin(beta*t)));
-plot(t,y,'-k'); hold on;
-
-
-
-% specifies the figure labels, the title, and the legend
-xlabel('time, t, sec');
-ylabel('underdamped unit-step response, y(t)');
-title('step response of an underdamped mechanical system');
-%title('step response of an underdamped second order system (\omega_n = 1 rad/s)');
-grid on;
-
-% exports figure to PNG format with 600 DPI
-print('step-response-mechanical-system.png','-r600','-dpng');
\ No newline at end of file
diff --git a/transient-response/examples/second-order-response-transient-response/impulse_response_underdamped_second_order_system.m b/transient-response/examples/second-order-response-transient-response/impulse_response_underdamped_second_order_system.m
deleted file mode 100644
index 7917b68..0000000
--- a/transient-response/examples/second-order-response-transient-response/impulse_response_underdamped_second_order_system.m
+++ /dev/null
@@ -1,49 +0,0 @@
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Course: System Dynamics and Control                       ME 3030-21 WI19
-% Author: Prof. M Diaz-Maldonado                          Date: Jan/22/2020
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Synopsis:
-% Plots the response of a second order system to a unit impulse input.
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-clear;
-close all;
-clc;
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% defines the second-order constants which specify a second-order system:
-% damping ratio
-DR = 0.15;
-% natural frequency
-Wn = 1;
-
-% real and imaginary parts of the roots of the characteristic equation
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% calculates the transient response, y(t)
-% Note: b0 = 1/Wn^2 and A=1.
-t = linspace(0,30,1e3);
-y = beta * (1+(alpha/beta)^2) * exp(alpha * t) .* sin(beta*t);
-plot(t,y,'-k'); hold on;
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% plots the impulse response for a higher damping ratio on the same figure
-DR = 0.3;
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-y = beta * (1+(alpha/beta)^2) * exp(alpha * t) .* sin(beta*t);
-plot(t,y,'-r');
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-
-% specifies the figure labels, the title, and the legend
-xlabel('time, t, sec');
-ylabel('underdamped impulse response, y(t)');
-title('impulse response of an underdamped second order system');
-%title('step response of an underdamped second order system (\omega_n = 1 rad/s)');
-grid on;
-
-legend('damping ratio \zeta = 0.15', 'damping ratio \zeta = 0.30');
-
-% exports figure to PNG format with 600 DPI
-print('impulse-response-second-order-system.png','-r600','-dpng');
diff --git a/transient-response/examples/second-order-response-transient-response/pulse_response_underdamped_second_order_system.m b/transient-response/examples/second-order-response-transient-response/pulse_response_underdamped_second_order_system.m
deleted file mode 100644
index 755aadb..0000000
--- a/transient-response/examples/second-order-response-transient-response/pulse_response_underdamped_second_order_system.m
+++ /dev/null
@@ -1,49 +0,0 @@
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Course: System Dynamics and Control                       ME 3030-21 WI19
-% Author: Prof. M Diaz-Maldonado                          Date: Dec/27/2019
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Synopsis:
-% Obtains the pulse response of an underdamped second-order system by
-% applying the superposition property. The weight of the pulse A = 1.
-% It's seen that as T -> 0 the response converges to the impulse response.
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-clear;
-close all;
-clc;
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% defines the second-order constants which specify a second-order system:
-% damping ratio
-DR = 0.2;
-% natural frequency
-Wn = 1;
-% pulse duration (sec)
-T = 0.5;
-% pulse intensity (note the pulse weight is A=1)
-P = 1/T;
-
-% real and imaginary parts of the roots of the characteristic equation
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% calculates the transient response, y(t), Equation (7.72)
-t = linspace(0,30,1e3);
-y = @(t) 1 - exp(alpha * t) .* (cos(beta*t) - alpha/beta * sin(beta*t));
-Y = P * ( y(t) - y(t-T) .* (t>T) );
-
-yim = beta * (1+(alpha/beta)^2) * exp(alpha * t) .* sin(beta*t);
-plot(t,Y,'-k'); hold on;
-plot(t,yim,'--');
-legend('pulse response', 'impulse response');
-
-% specifies the figure labels, the title, and the legend
-ylim([-0.5, 1.0])
-xlabel('time, t, sec');
-ylabel('transient response, y(t)');
-title('pulse response of an underdamped second order system');
-%title('step response of an underdamped second order system (\omega_n = 1 rad/s)');
-grid on;
-
-% exports figure to PNG format with 600 DPI
-print('pulse-response-second-order-system.png','-r600','-dpng');
\ No newline at end of file
diff --git a/transient-response/examples/second-order-response-transient-response/step_response_underdamped_second_order_system.m b/transient-response/examples/second-order-response-transient-response/step_response_underdamped_second_order_system.m
deleted file mode 100644
index 73e5a1a..0000000
--- a/transient-response/examples/second-order-response-transient-response/step_response_underdamped_second_order_system.m
+++ /dev/null
@@ -1,50 +0,0 @@
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Course: System Dynamics and Control                       ME 3030-21 WI19
-% Author: Prof. M Diaz-Maldonado                          Date: Dec/27/2019
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% Synopsis:
-% Obtains the unit-step response of an underdamped second-order system.
-% (Reproduces Figure 7.19 of Kluever's textbook.)
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-clear;
-close all;
-clc;
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% defines the second-order constants which specify a second-order system:
-% damping ratio
-DR = 0.2;
-% natural frequency
-Wn = 1;
-
-% real and imaginary parts of the roots of the characteristic equation
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% calculates the transient response, y(t), Equation (7.72)
-t = linspace(0,25,500);
-y = 1 - exp(alpha * t) .* (cos(beta*t) - alpha/beta * sin(beta*t));
-plot(t,y,'-k'); hold on;
-
-
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-% plots the step response for a higher damping ratio on the same figure
-DR = 0.4;
-alpha = -DR * Wn;
-beta  = Wn * sqrt(1-DR^2);
-y = 1 - exp(alpha * t) .* (cos(beta*t) -alpha/beta * sin(beta*t));
-plot(t,y,'-r');
-% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
-
-
-% specifies the figure labels, the title, and the legend
-xlabel('time, t, sec');
-ylabel('underdamped unit-step response, y(t)');
-title('step response of an underdamped second order system');
-%title('step response of an underdamped second order system (\omega_n = 1 rad/s)');
-grid on;
-legend('damping ratio \zeta = 0.2', 'damping ratio \zeta = 0.4');
-
-% exports figure to PNG format with 600 DPI
-print('unit-step-response-second-order-system.png','-r600','-dpng');
\ No newline at end of file