In system dynamics, time enters the model in one of two fundamental ways: continuously or discretely. A continuous-time (CT) signal is modeled as a function of a real-valued time variable \( t \), while a discrete-time (DT) signal is modeled as a sequence indexed by an integer \( k \). Formally, for an \( n \)-dimensional state:

\[ x_c : \mathbb{R}_{\ge 0} \rightarrow \mathbb{R}^n, \quad t \mapsto x_c(t), \qquad x_d : \mathbb{Z}_{\ge 0} \rightarrow \mathbb{R}^n, \quad k \mapsto x_d[k]. \]

In engineering, physical systems are inherently continuous in time, but measurements, simulations, and digital controllers operate on discrete samples. A central task in system dynamics is therefore to:

• Model physics with continuous-time differential equations.
• Derive discrete-time difference equations that approximate or exactly reproduce the continuous dynamics at sampling instants.
• Understand similarities and differences in behavior between the two models.

flowchart TD
  P["Physical system \n(continuous in time)"] --> M["Choose modeling viewpoint"]
  M --> C["Continuous-time model \n(differential equations)"]
  M --> D["Discrete-time model \n(difference equations)"]
  C --> AC["Analyze: calculus and \ndifferential equations"]
  D --> AD["Analyze: sequences \nand recursions"]
  C --> SC["Simulate: \nnumerical integration"]
  D --> SD["Simulate: \niterate x[k+1] = ..."]
        

In this lesson we stay within basic ordinary differential equation (ODE) theory and algebra. Transform-domain tools (Laplace, \( z \)-transform) and full state-space theory will be developed in later chapters.

A continuous-time dynamic model describes the evolution of a state vector \( x(t) \in \mathbb{R}^n \) driven by an input \( u(t) \in \mathbb{R}^m \) using an ODE:

\[ \frac{d}{dt} x(t) = f\big(x(t), u(t), t\big), \quad t \ge 0, \qquad x(0) = x_0. \]

For many engineering systems, a linear time-invariant (LTI) model is adequate around a given operating point. In that case the dynamics are:

\[ \frac{d}{dt} x(t) = A x(t) + B u(t), \qquad y(t) = C x(t) + D u(t), \]

where \( A \in \mathbb{R}^{n \times n} \), \( B \in \mathbb{R}^{n \times m} \), \( C \in \mathbb{R}^{p \times n} \), \( D \in \mathbb{R}^{p \times m} \), and \( y(t) \in \mathbb{R}^p \) is the output. The simplest nontrivial example is a first-order SISO (single-input single-output) system:

\[ \frac{dx(t)}{dt} = -\frac{1}{\tau}\,x(t) + \frac{K}{\tau}\,u(t), \quad \tau > 0. \]

Here, \( K \) is a static gain and \( \tau \) is a time constant. You have already seen how to solve such linear ODEs in mathematics: for constant input \( u(t) \equiv u_0 \), the solution is an exponential plus a constant.

A discrete-time signal is a sequence \( x_d[0], x_d[1], x_d[2], \dots \) defined for integer indices. In practice, discrete-time signals typically arise from sampling a continuous signal with sampling period \( T > 0 \):

\[ t_k = k T, \quad k \in \mathbb{Z}_{\ge 0}, \qquad x_d[k] = x_c(t_k). \]

A general nonlinear discrete-time model takes the form

\[ x[k+1] = g\big(x[k], u[k], k\big), \qquad x[0] = x_0, \]

where \( x[k] \in \mathbb{R}^n \) and \( u[k] \in \mathbb{R}^m \). For LTI systems, the dynamics are linear:

\[ x[k+1] = A_d x[k] + B_d u[k], \qquad y[k] = C_d x[k] + D_d u[k], \]

with constant matrices \( A_d, B_d, C_d, D_d \). A simple first-order SISO case is:

\[ x[k+1] = a\,x[k] + b\,u[k], \qquad y[k] = c\,x[k] + d\,u[k]. \]

The central question of this lesson: given a continuous-time model, how do we obtain meaningful discrete-time parameters \( a, b \) (or \( A_d, B_d \))? Conversely, when is it reasonable to regard a given discrete-time model as an approximation of continuous dynamics?

Consider a scalar continuous-time LTI system

\[ \frac{dx(t)}{dt} = a\,x(t) + b\,u(t), \qquad x(0) = x_0, \]

with constant coefficients \( a, b \in \mathbb{R} \). Suppose we sample every \( T > 0 \) seconds at times \( t_k = kT \), and assume a zero-order hold on the input:

\[ u(t) = u[k] \quad \text{for } t_k \le t < t_{k+1}, \qquad t_k = kT. \]

By standard ODE theory (variation of constants), for any interval \( [t_k, t_{k+1}] \) we have the exact solution

\[ x(t) = e^{a(t-t_k)} x(t_k) + \int_{t_k}^{t} e^{a(t-\sigma)} b\,u(\sigma)\,d\sigma. \]

Evaluate this at \( t = t_{k+1} = t_k + T \) and use the zero-order hold assumption \( u(\sigma) = u[k] \):

\[ \begin{aligned} x(t_{k+1}) &= e^{aT} x(t_k) + b\,u[k] \int_{t_k}^{t_k+T} e^{a(t_k+T-\sigma)} d\sigma \\ &= e^{aT} x(t_k) + b\,u[k] \int_{0}^{T} e^{a(T-\theta)} d\theta \quad (\theta = \sigma - t_k) \\ &= e^{aT} x(t_k) + b\,u[k] \, e^{aT} \int_{0}^{T} e^{-a\theta} d\theta. \end{aligned} \]

If \( a \neq 0 \), we integrate the exponential explicitly:

\[ \int_{0}^{T} e^{-a\theta}\,d\theta = \frac{1 - e^{-aT}}{a}, \]

so that

\[ x(t_{k+1}) = e^{aT} x(t_k) + b\,u[k]\,e^{aT}\,\frac{1 - e^{-aT}}{a} = e^{aT} x(t_k) + \frac{e^{aT}-1}{a}\,b\,u[k]. \]

Defining discrete-time variables \( x[k] := x(t_k) \), \( u[k] \) as above, we obtain the exact discrete-time model:

\[ x[k+1] = e^{aT} x[k] + \frac{e^{aT}-1}{a}\,b\,u[k], \qquad a \neq 0. \]

For the special case \( a = 0 \), the continuous model reduces to \( dx/dt = b\,u(t) \). With zero-order hold, the integral becomes

\[ x(t_{k+1}) = x(t_k) + b\,u[k] \int_{t_k}^{t_k+T} d\sigma = x(t_k) + b\,T\,u[k], \]

giving the discrete-time recursion

\[ x[k+1] = x[k] + b\,T\,u[k], \qquad a = 0. \]

For vector-valued systems with constant matrices \( A, B \), the same reasoning using the matrix exponential \( e^{AT} \) yields

\[ A_d = e^{AT}, \qquad B_d = \int_{0}^{T} e^{A\sigma} B\,d\sigma, \]

which defines an exact discrete-time model at sampling instants under zero-order hold.

flowchart TD
  CT["Continuous model: dx/dt = f(x,u,t)"] --> Tsel["Choose sampling period T"]
  Tsel --> Meth["Select discretization method"]
  Meth --> Exact["Exact: compute Phi = e^{AT}, Gamma = integral"]
  Meth --> Euler["Approximate: Euler or higher-order method"]
  Exact --> Rec["Recursion x[k+1] = Phi x[k] + Gamma u[k]"]
  Euler --> Rec
  Rec --> Sim["Iterate for k = 0..N to simulate"]
        

Exact formulas \( A_d = e^{AT} \), \( B_d = \int_0^T e^{A\sigma} B\,d\sigma \) may be computationally expensive or analytically intractable. A widely used alternative is to approximate the derivative by a finite difference:

\[ \dot{x}(t_k) \approx \frac{x(t_{k+1}) - x(t_k)}{T}, \qquad t_k = kT. \]

Inserting this into the ODE \( \dot{x} = f(x,u,t) \) leads to the forward Euler scheme:

\[ x[k+1] = x[k] + T\,f\big(x[k], u[k], t_k\big). \]

To understand the accuracy, use Taylor expansion of the exact solution around \( t_k \):

\[ x(t_{k+1}) = x(t_k) + T\,\dot{x}(t_k) + \frac{T^2}{2}\,\ddot{x}(\xi_k), \quad \text{for some } \xi_k \in (t_k, t_{k+1}). \]

With \( \dot{x}(t_k) = f(x(t_k),u(t_k),t_k) \), the local truncation error (one-step error assuming exact \( x(t_k) \)) is

\[ x(t_{k+1}) - \Bigl[x(t_k) + T f\big(x(t_k),u(t_k),t_k\big)\Bigr] = \mathcal{O}(T^2), \]

i.e., forward Euler has first-order accuracy (local error proportional to \( T^2 \), global error roughly proportional to \( T \)).

For the scalar linear system \( \dot{x}(t) = a\,x(t) \) (with \( u \equiv 0 \)), forward Euler yields the discrete-time recursion

\[ x[k+1] = x[k] + T a x[k] = (1 + aT)\,x[k]. \]

By induction,

\[ x[k] = (1 + aT)^k x[0], \]

which approximates the exact solution \( x(t) = e^{at} x(0) \) at \( t = kT \) as \( (1 + aT)^k \approx e^{a kT} \) for small \( T \).

Even before a full treatment of stability, it is instructive to see that discretization can change stability properties. Consider again the scalar ODE \( \dot{x}(t) = a\,x(t) \) with \( a < 0 \). The continuous-time solution \( x(t) = e^{at} x(0) \) decays to zero as \( t \rightarrow \infty \).

Under forward Euler, the discrete-time recursion is

\[ x[k+1] = (1 + aT)\,x[k]. \]

The discrete-time solution is \( x[k] = (1 + aT)^k x[0] \). For this to decay to zero as \( k \rightarrow \infty \), we require

\[ |1 + aT| < 1. \]

This inequality is equivalent to

\[ -1 < 1 + aT < 1 \quad \Longleftrightarrow \quad -2 < aT < 0 \quad \Longleftrightarrow \quad 0 < T < -\frac{2}{a}. \]

Because \( a < 0 \), the upper bound \( -2/a \) is positive. Thus a continuous-time stable system can become unstable under forward Euler if the time step \( T \) is chosen too large. This illustrates that discrete-time models are not merely algebraic rewritings of continuous models: the discretization method and step size matter.

We now implement a first-order system with time constant \( \tau \) and unit gain:

\[ \frac{dx(t)}{dt} = -\frac{1}{\tau}\,x(t) + \frac{1}{\tau}\,u(t), \quad u(t) \equiv 1,\; x(0) = 0. \]

The exact continuous-time solution is

\[ x(t) = 1 - e^{-t/\tau}. \]

The exact discrete-time model at sampling period \( T \) with zero-order hold is

\[ x[k+1] = e^{-T/\tau}\,x[k] + \bigl(1 - e^{-T/\tau}\bigr) u[k], \quad u[k] \equiv 1. \]

Forward Euler gives instead

\[ x[k+1] = \Bigl(1 - \frac{T}{\tau}\Bigr)x[k] + \frac{T}{\tau} u[k]. \]

The following Python script compares continuous-time analytical values with both exact discrete-time and Euler approximations. We use numpy and matplotlib; in more advanced work one would also use scipy.integrate or the python-control package.

import numpy as np
import matplotlib.pyplot as plt

tau = 1.0      # time constant
T = 0.1        # sampling period
t_final = 5.0  # simulation horizon

# time grid for samples
k_max = int(t_final / T)
t_samples = np.arange(0.0, (k_max + 1) * T, T)

# continuous-time analytical solution
x_cont = 1.0 - np.exp(-t_samples / tau)

# exact discrete-time recursion coefficients
a_exact = np.exp(-T / tau)
b_exact = 1.0 - np.exp(-T / tau)

# Euler discrete-time recursion coefficients
a_euler = 1.0 - T / tau
b_euler = T / tau

# simulate exact discrete-time model
x_exact = np.zeros_like(t_samples)
# simulate Euler approximation
x_euler = np.zeros_like(t_samples)

u = 1.0  # unit step input

for k in range(k_max):
    x_exact[k + 1] = a_exact * x_exact[k] + b_exact * u
    x_euler[k + 1] = a_euler * x_euler[k] + b_euler * u

# plot
plt.figure()
plt.plot(t_samples, x_cont, label="continuous analytic")
plt.stem(t_samples, x_exact, linefmt="C1-", markerfmt="C1o", basefmt=" ", label="exact discrete", use_line_collection=True)
plt.stem(t_samples, x_euler, linefmt="C2-", markerfmt="C2s", basefmt=" ", label="Euler discrete", use_line_collection=True)
plt.xlabel("time t [s]")
plt.ylabel("state x")
plt.legend()
plt.grid(True)
plt.show()
      

By decreasing \( T \), the Euler approximation approaches both the continuous-time solution and the exact discrete-time model. However, if \( T \) violates the stability bound in Section 6, the Euler solution diverges even though the continuous-time system is stable.

In C++, we can implement the same discrete-time recursion using basic standard library facilities. For more advanced projects, linear algebra libraries such as Eigen and ODE solvers such as boost::numeric::odeint are commonly used in system dynamics.

#include <iostream>
#include <vector>
#include <cmath>

int main() {
    double tau = 1.0;
    double T = 0.1;
    double t_final = 5.0;
    int k_max = static_cast<int>(t_final / T);

    // Euler coefficients for dx/dt = -(1/tau) x + (1/tau) u, u = 1
    double a_euler = 1.0 - T / tau;
    double b_euler = T / tau;

    std::vector<double> t(k_max + 1);
    std::vector<double> x(k_max + 1);
    std::vector<double> x_cont(k_max + 1);

    x[0] = 0.0;
    for (int k = 0; k <= k_max; ++k) {
        t[k] = k * T;
        x_cont[k] = 1.0 - std::exp(-t[k] / tau);  // continuous-time reference
        if (k < k_max) {
            double u = 1.0;
            x[k + 1] = a_euler * x[k] + b_euler * u;
        }
    }

    std::cout << "k\t t\t x_euler\t x_cont\n";
    for (int k = 0; k <= k_max; ++k) {
        std::cout << k << "\t" << t[k]
                  << "\t" << x[k]
                  << "\t" << x_cont[k] << "\n";
    }

    return 0;
}
      

Compiling and running this program prints the Euler approximation alongside the continuous-time analytical solution at the sampling instants. Plotting these values (e.g., using gnuplot) gives a visual comparison similar to the Python figure.

In Java, we may wrap the discrete-time recursion in a simple class. Libraries such as Apache Commons Math provide ODE solvers and linear algebra tools, but here we implement the recursion directly.

public class FirstOrderDiscreteSystem {
    private final double a;
    private final double b;
    private double x;

    public FirstOrderDiscreteSystem(double a, double b, double x0) {
        this.a = a;
        this.b = b;
        this.x = x0;
    }

    public double step(double u) {
        // x[k+1] = a x[k] + b u[k]
        x = a * x + b * u;
        return x;
    }

    public double getState() {
        return x;
    }

    public static void main(String[] args) {
        double tau = 1.0;
        double T = 0.1;
        // Euler coefficients
        double aEuler = 1.0 - T / tau;
        double bEuler = T / tau;

        FirstOrderDiscreteSystem sys = new FirstOrderDiscreteSystem(aEuler, bEuler, 0.0);

        int kMax = (int) (5.0 / T);
        double t = 0.0;
        System.out.println("k\t t\t x_euler");
        for (int k = 0; k <= kMax; ++k) {
            if (k > 0) {
                sys.step(1.0); // unit step input
            }
            System.out.println(k + "\t " + t + "\t " + sys.getState());
            t += T;
        }
    }
}
      

This object-oriented viewpoint will be helpful when building larger simulation frameworks or integrating discrete-time plant models with discrete-time controllers in later courses.

MATLAB provides built-in ODE solvers and direct support for discrete-time simulations. Simulink offers block diagrams for continuous and discrete integrators and zero-order hold components. The script below compares:

• Continuous-time solution using ode45.
• Discrete-time Euler approximation in a loop.

tau = 1.0;
T = 0.1;
t_final = 5.0;

% continuous-time dynamics for ode45
f = @(t, x) -(1 / tau) * x + (1 / tau) * 1.0;  % unit step input

[t_cont, x_cont] = ode45(f, [0 t_final], 0.0);

% discrete-time Euler simulation
k_max = floor(t_final / T);
t_k = (0:k_max)' * T;
x_euler = zeros(k_max + 1, 1);
a_euler = 1.0 - T / tau;
b_euler = T / tau;

for k = 1:k_max
    u = 1.0;
    x_euler(k + 1) = a_euler * x_euler(k) + b_euler * u;
end

% plot comparison
figure;
plot(t_cont, x_cont, 'LineWidth', 1.5); hold on;
stairs(t_k, x_euler, 'r--', 'LineWidth', 1.2);
xlabel('time t [s]');
ylabel('state x');
legend('continuous (ode45)', 'discrete Euler');
grid on;
title('Continuous-time vs. Euler discrete-time response');

% Simulink note:
% A Simulink model would contain:
% - a first-order Transfer Fcn or State-Space block (continuous),
% - a zero-order hold and discrete transfer function for the discrete model,
% both driven by the same step input for visual comparison.
      

Simulink's explicit continuous-time and discrete-time blocks make the conceptual distinction between the two models very clear and allow easy experimentation with different sampling periods and discretization methods.

Wolfram Mathematica has symbolic and numeric tools for both differential and difference equations. The following snippet compares the analytic solution with a discrete-time Euler approximation:

tau = 1.0;
T = 0.1;
tFinal = 5.0;

(* continuous-time solution of x'(t) = -(1/tau) x(t) + (1/tau), x(0) = 0 *)
sol = First@
   NDSolve[{x'[t] == -(1/tau) x[t] + (1/tau), x[0] == 0},
           x, {t, 0, tFinal}];

xCont[t_] := x[t] /. sol;

(* discrete-time Euler coefficients *)
aEuler = 1.0 - T/tau;
bEuler = T/tau;

kMax = Floor[tFinal/T];
tSamples = Range[0, kMax] T;

xEuler = NestList[aEuler #1 + bEuler &, 0.0, kMax];

ListLinePlot[
 {
  Table[{t, xCont[t]}, {t, 0, tFinal, T}],
  Transpose[{tSamples, xEuler}]
 },
 PlotLegends -> {"continuous", "Euler discrete"},
 AxesLabel -> {"t", "x"},
 PlotMarkers -> {None, Automatic},
 GridLines -> Automatic
]
      

Mathematica can also manipulate the exact exponential \( e^{-T/\tau} \) symbolically, derive discrete-time transfer functions (later chapters), and explore parameter dependencies analytically.

Problem 1 (Exact Discretization of a Stable First-Order System). Consider the continuous-time system \( \dot{x}(t) = -\alpha x(t) \) with \( \alpha > 0 \). Sample with period \( T > 0 \) under zero-order hold (no input). Derive the exact discrete-time model and compare its solution with the continuous-time solution at sampling instants.

Solution. The continuous-time solution is

\[ x(t) = e^{-\alpha t} x(0). \]

At sampling instants \( t_k = kT \),

\[ x(t_k) = e^{-\alpha kT} x(0). \]

Define \( x[k] := x(t_k) \). Then

\[ x[k+1] = e^{-\alpha (k+1)T} x(0) = e^{-\alpha T} e^{-\alpha kT} x(0) = e^{-\alpha T} x[k]. \]

Thus the exact discrete-time model is \( x[k+1] = a x[k] \) with \( a = e^{-\alpha T} \in (0,1) \). The discrete-time solution is \( x[k] = a^k x[0] = e^{-\alpha kT} x(0) \), which matches the continuous-time solution evaluated at sample times.

Problem 2 (Euler Approximation Error for a First-Order System). Using the same system as Problem 1, discretize it via forward Euler: \( x[k+1] = x[k] + T \dot{x}(t_k) \). Derive the Euler recursion and compute the one-step local truncation error at \( t_k \).

Solution. With \( \dot{x}(t) = -\alpha x(t) \), forward Euler gives

\[ x[k+1] = x[k] + T(-\alpha x[k]) = (1 - \alpha T)\,x[k]. \]

The exact solution satisfies

\[ x(t_{k+1}) = e^{-\alpha T} x(t_k). \]

Define the local truncation error (assuming exact \( x(t_k) \)) as

\[ \ell_k = x(t_{k+1}) - \bigl[x(t_k) + T \dot{x}(t_k)\bigr] = e^{-\alpha T} x(t_k) - \bigl[1 - \alpha T\bigr] x(t_k). \]

Expanding \( e^{-\alpha T} \) in a Taylor series,

\[ e^{-\alpha T} = 1 - \alpha T + \frac{\alpha^2 T^2}{2} + \mathcal{O}(T^3), \]

gives

\[ \ell_k = \Bigl(1 - \alpha T + \frac{\alpha^2 T^2}{2} + \mathcal{O}(T^3) - (1 - \alpha T)\Bigr) x(t_k) = \frac{\alpha^2 T^2}{2} x(t_k) + \mathcal{O}(T^3). \]

Hence the one-step local error is \( \mathcal{O}(T^2) \), confirming first-order accuracy.

Problem 3 (Step-Size Condition for Euler Stability). For the system in Problem 1, show that forward Euler is stable (i.e., \( x[k] \rightarrow 0 \) as \( k \rightarrow \infty \)) if and only if \( 0 < T < \frac{2}{\alpha} \).

Solution. From Problem 2, Euler gives

\[ x[k+1] = (1 - \alpha T)\,x[k]. \]

The discrete-time solution is \( x[k] = (1 - \alpha T)^k x[0] \). This converges to zero if and only if

\[ |1 - \alpha T| < 1. \]

Solving:

\[ -1 < 1 - \alpha T < 1 \quad \Longleftrightarrow \quad -2 < -\alpha T < 0 \quad \Longleftrightarrow \quad 0 < \alpha T < 2 \quad \Longleftrightarrow \quad 0 < T < \frac{2}{\alpha}. \]

Thus forward Euler imposes a maximum step size for stability, whereas the original continuous-time system is stable for all \( t \ge 0 \).

Problem 4 (Exact vs. Euler Discrete Parameters for a First-Order System). Consider the system \( \dot{x}(t) = -\frac{1}{\tau}\,x(t) + \frac{1}{\tau}\,u(t) \) with zero-order hold on \( u(t) \). For sampling period \( T \), derive:

1. The exact discrete-time coefficients \( a_{\text{exact}}, b_{\text{exact}} \) such that \( x[k+1] = a_{\text{exact}} x[k] + b_{\text{exact}} u[k] \).
2. The Euler coefficients \( a_{\text{Euler}}, b_{\text{Euler}} \) from forward Euler.
3. The leading-order difference \( a_{\text{exact}} - a_{\text{Euler}} \) for small \( T \).

Solution.

(1) The scalar system has \( a = -1/\tau \), \( b = 1/\tau \) in the notation of Section 4. Using the exact discretization formula with \( a \neq 0 \),

\[ a_{\text{exact}} = e^{aT} = e^{-T/\tau}, \qquad b_{\text{exact}} = \frac{e^{aT}-1}{a}\,b = \frac{e^{-T/\tau}-1}{-1/\tau} \frac{1}{\tau} = 1 - e^{-T/\tau}. \]

(2) Forward Euler yields

\[ x[k+1] = x[k] + T\Bigl(-\frac{1}{\tau} x[k] + \frac{1}{\tau} u[k]\Bigr) = \Bigl(1 - \frac{T}{\tau}\Bigr) x[k] + \frac{T}{\tau} u[k], \]

hence

\[ a_{\text{Euler}} = 1 - \frac{T}{\tau}, \qquad b_{\text{Euler}} = \frac{T}{\tau}. \]

(3) Expand \( a_{\text{exact}} = e^{-T/\tau} \) for small \( T \):

\[ a_{\text{exact}} = 1 - \frac{T}{\tau} + \frac{T^2}{2 \tau^2} + \mathcal{O}(T^3). \]

Therefore

\[ a_{\text{exact}} - a_{\text{Euler}} = \Bigl(1 - \frac{T}{\tau} + \frac{T^2}{2 \tau^2} + \mathcal{O}(T^3)\Bigr) - \Bigl(1 - \frac{T}{\tau}\Bigr) = \frac{T^2}{2 \tau^2} + \mathcal{O}(T^3), \]

confirming that the discrepancy between exact and Euler coefficients is of order \( T^2 \) for small \( T \).

In this lesson we distinguished continuous-time models (ODEs) from discrete-time models (difference equations) and showed how discrete-time descriptions arise naturally through sampling and zero-order hold. For linear time-invariant systems, exact discretization at sampling period \( T \) leads to matrices \( A_d = e^{AT} \) and \( B_d = \int_0^T e^{A\sigma} B\,d\sigma \), while numerical schemes such as forward Euler provide approximate recursions whose accuracy and stability depend on the step size. Through implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica, we illustrated how the same physical first-order system can be simulated both continuously and in discrete time, highlighting the interplay between analytical solutions and numerical approximations. These ideas are foundational for later chapters on transfer functions, block diagrams, state-space models, and sampled-data system dynamics.