First-order System MATLAB scripts

This page allows you to download the MATLAB scripts for simulations in First-order Systems Lecture.

Just run the cell of the script you want or “Run all” to download everything and they will be saved under “Files”.

Back to Lecture.

1. Cruise Control

%%writefile cruise_control.m

% ======================================================
% Solve the first-order differential equation
%       e_dot = A * e + B * v_d(t)
% for a given desired speed profile v_d(t) and multiple initial conditions.
%
% User specifies:
%   - A, B: system parameters
%   - x0_vec: list of initial error values
% ======================================================

clear; close all; clc;

% ======================================================
% User-defined parameters
% ======================================================
A = -1.0;
B = 0.1;
x0_vec = [-1.0, -0.5, 0.0, 0.5, 1.0];   % initial error values [m/s]
t_span = [0 20];                         % time interval [s]
t_eval = linspace(t_span(1), t_span(2), 500);

% ======================================================
% Desired speed profile (as a function handle)
% ======================================================
v_d = @(t) (t < 5).*40.0 + (t >= 5).*50.0;

% ======================================================
% Differential equation (as a function handle)
% ======================================================
e_dot = @(t, e) A*e + B*v_d(t);

% ======================================================
% Solve for each initial condition
% ======================================================
solutions = cell(length(x0_vec),1);
for i = 1:length(x0_vec)
    x0 = x0_vec(i);
    % Solve using ode45 (similar to solve_ivp)
    [t_sol, e_sol] = ode45(e_dot, t_span, x0);
    % Interpolate to uniform t_eval points
    e_interp = interp1(t_sol, e_sol, t_eval, 'linear');
    solutions{i} = struct('t', t_eval, 'e', e_interp, 'x0', x0);
end

% ======================================================
% Plot 1: Desired speed profile
% ======================================================
figure('Name', 'Desired Speed Profile', 'Color', 'w');
plot(t_eval, v_d(t_eval), 'k', 'LineWidth', 2);
title('Desired Speed $v\_d(t)$', 'Interpreter', 'latex');
xlabel('Time [s]', 'Interpreter', 'latex');
ylabel('Desired Speed [m/s]', 'Interpreter', 'latex');
grid on; box on;
ylim([38 52]);

% ======================================================
% Plot 2: Error trajectories
% ======================================================
figure('Name', 'Error Trajectories', 'Color', 'w');
hold on;
for i = 1:length(solutions)
    plot(solutions{i}.t, solutions{i}.e, 'LineWidth', 1.5, ...
        'DisplayName', sprintf('$e(0) = %.1f$', solutions{i}.x0));
end
title('Evolution of Speed Error $e(t)$', 'Interpreter', 'latex');
xlabel('Time [s]', 'Interpreter', 'latex');
ylabel('Error $e(t)$ [m/s]', 'Interpreter', 'latex');
grid on; box on; legend('Interpreter', 'latex', 'Location', 'best');
xlim([-1,22])
ylim([-1.2,5.5])
hold off;

Overwriting cruise_control.m

2. Free Response

%%writefile free_response_first_order.m
% ======================================================
% Free Response of a First-Order System
%   x_dot = A * x
%   Solution: x(t) = e^(A t) * x0
%
% For multiple initial conditions x0 and different A values.
% ======================================================

clear; close all; clc;

% === Parameters ===
x0_values = [-2, -1, 0, 1, 2];
t = linspace(0, 10, 200);
A_values = [-1, 0, 1];

% === Main loop over A values ===
for A = A_values
    figure('Color', 'w'); hold on;

    % Compute and plot responses for each initial condition
    for x0 = x0_values
        x = exp(A * t) * x0;
        plot(t, x, 'LineWidth', 1.5, 'DisplayName', sprintf('$x_0 = %.1f$', x0));
    end

    % Titles and labels
    title(sprintf('Free Response for $A = %.1f$: $x(t) = e^{A t} x_0$', A), ...
        'Interpreter', 'latex');
    xlabel('Time $t$', 'Interpreter', 'latex');
    ylabel('$x(t)$', 'Interpreter', 'latex');
    grid on; box on;

    % Legend placement depending on A
    if A > 0
        legend('Interpreter', 'latex', 'Location', 'northwest', ...
            'Box', 'on', 'Color', 'white', 'EdgeColor', [0.6 0.6 0.6]);
    else
        legend('Interpreter', 'latex', 'Location', 'northeast', ...
            'Box', 'on', 'Color', 'white', 'EdgeColor', [0.6 0.6 0.6]);
    end

    % === Add manual y-axis padding ===
    yl = ylim;
    dy = diff(yl);
    ylim([yl(1) - 0.05*dy, yl(2) + 0.05*dy]);

    hold off;
end

Overwriting free_response_first_order.m

3. Time constant response

%%writefile time_constant_response.m
% ======================================================
% Free Response and Time Constant Visualization
% ======================================================

clear; close all; clc;

x0 = 1.0;
t = linspace(0, 10, 400);
A_values = [-0.5, -1, -2];
colors = {'b', [0.8500 0.3250 0.0980], [0.4660 0.6740 0.1880]};

figure('Color', 'w'); hold on;

for i = 1:length(A_values)
    A = A_values(i);
    color = colors{i};
    T = -1 / A;
    x = exp(A * t) * x0;

    % --- main curve (legend entry) ---
    plot(t, x, 'Color', color, 'LineWidth', 1.5, ...
         'DisplayName', sprintf('$A = %.1f$, $T = %.2f$', A, T));

    % --- auxiliary markers (no legend entry) ---
    for k = [1, 5]
        tT = k * T;
        xT = exp(A * tT) * x0;
        plot([tT tT], ylim, '--', 'Color', color, 'LineWidth', 0.9, ...
             'HandleVisibility', 'off');
        plot(tT, xT, 'o', 'Color', color, 'MarkerFaceColor', color, ...
             'HandleVisibility', 'off');
        text(tT + 0.1, xT, sprintf('%dT', k), 'Color', color, ...
             'VerticalAlignment', 'middle', 'FontSize', 9);
    end
end

% --- horizontal reference lines (no legend entry) ---
for k = 1:3
    yref = exp(-k);
    yline(yref, ':', 'Color', 'r', 'LineWidth', 0.8, 'HandleVisibility', 'off');
    text(8.15, yref, sprintf('$e^{-%d}$', k), 'Color', 'r', ...
         'Interpreter', 'latex', 'VerticalAlignment', 'middle', 'FontSize', 9);
end

title('Free Response and Time Constant: $x(t) = e^{A t} x_0$', 'Interpreter', 'latex');
xlabel('Time $t$', 'Interpreter', 'latex');
ylabel('$x(t)$', 'Interpreter', 'latex');
grid on; box on;

legend('Interpreter', 'latex', 'Location', 'northeast'); % clean legend

% === Add manual y-axis padding ===
yl = ylim; dy = diff(yl);
ylim([yl(1) - 0.05*dy, yl(2) + 0.05*dy]);

hold off;

Overwriting time_constant_response.m

4. Step response with variable gain

%%writefile step_response_multi_B.m
% ======================================================
% Step Response of a First-Order System for Multiple B Values
%   System: x_dot = A*x + B*u(t)
%   Step input: u(t) = u_m for t >= 0
%
% Consistent notation:
%   x(t) = -(B/A) * u_m * (1 - e^{A t})
%
% Shows input step (t<0) and outputs for several B values.
% ======================================================

clear; close all; clc;

% === User parameters ===
A = -1.0;
B_values = [0.5, 1.0, 2.0];
u_m = 1.0;
t_pre = 2.0;
t_end = 6.0;
n_pts = 1200;

if A >= 0
    error('Use A < 0 for a stable first-order system.');
end
if isempty(B_values)
    error('Provide at least one value in B_values.');
end

% === Time setup ===
t_full = linspace(-t_pre, t_end, n_pts);
t = linspace(0, t_end, n_pts);
u = zeros(size(t_full));
u(t_full >= 0) = u_m;
T = 1 / abs(A); % time constant

% ======================================================
% Figure setup — wider layout
% ======================================================
figure('Color','w','Units','normalized','Position',[0.05 0.1 0.85 0.45]);
tiledlayout(1,2,'Padding','compact','TileSpacing','loose');

% ======================================================
% LEFT: Step input u(t)
% ======================================================
nexttile; hold on;
plot(t_full, u, 'LineWidth', 2, 'Color', [0 0.447 0.741], ...
    'DisplayName', sprintf('Step input $u_m=%.1f$', u_m));
yline(0,'k-','LineWidth',0.8,'HandleVisibility','off');
xline(0,'k-','LineWidth',0.8,'HandleVisibility','off');
xline(T,'--','LineWidth',1,'Color','k','HandleVisibility','off');
text(T, u_m*0.1, '$t=T$', 'Interpreter','latex', ...
    'HorizontalAlignment','center', 'VerticalAlignment','bottom');

xlim([-t_pre t_end]);
ylim([-0.2*u_m 1.2*u_m]);
title('Input step $u(t)$ (showing $t<0$ region)', 'Interpreter','latex');
xlabel('Time $t$', 'Interpreter','latex');
ylabel('$u(t)$', 'Interpreter','latex');
grid on; box on;
legend('Interpreter','latex','FontSize',9,'Location','best');

% --- Add y-axis padding ---
yl = ylim; dy = diff(yl);
ylim([yl(1) - 0.05*dy, yl(2) + 0.05*dy]);

% ======================================================
% RIGHT: Outputs x(t) for multiple B values
% ======================================================
nexttile; hold on;
cmap = lines(numel(B_values));

for i = 1:length(B_values)
    B = B_values(i);
    color = cmap(i,:);
    x_ss = -(B/A)*u_m;
    x = x_ss * (1 - exp(A*t));
    plot(t, x, 'Color', color, 'LineWidth', 2, ...
        'DisplayName', sprintf('$B=%.1f,\\ x_{ss}=%.1f$', B, x_ss));
    yline(x_ss, ':', 'LineWidth',1, 'Color', color, ...
        'HandleVisibility','off');
end

% 63.2% marker for first curve (no legend entry)
B0 = B_values(1);
xss0 = -(B0/A)*u_m;
plot(T, xss0*(1 - exp(-1)), 'o', 'MarkerSize',5, ...
    'MarkerFaceColor', cmap(1,:), 'Color', cmap(1,:), ...
    'HandleVisibility','off');
text(T*1.05, xss0*(1 - exp(-1)), '$63.2\%$ of $x_{ss}$ at $t=T$', ...
    'Interpreter','latex', 'HorizontalAlignment','left', ...
    'VerticalAlignment','middle', 'FontSize',9, 'Color', cmap(1,:));

yline(0,'k-','LineWidth',0.8,'HandleVisibility','off');
xline(T,'--','LineWidth',1,'Color','k','HandleVisibility','off');

title('Outputs $x(t)$ for different $B$ values', 'Interpreter','latex');
xlabel('Time $t$', 'Interpreter','latex');
ylabel('$x(t)$', 'Interpreter','latex');
grid on; box on;
legend('Interpreter','latex','FontSize',9,'Location','best');

% --- Add y-axis padding ---
yl = ylim; dy = diff(yl);
ylim([yl(1) - 0.05*dy, yl(2) + 0.05*dy]);

hold off;

Overwriting step_response_multi_B.m

5. Complete solution = Free Response + Forced Response

%%writefile decomposition_first_order.m
% ======================================================
% Decomposition of Total Response for a First-Order System
%
%   x_dot = A*x + B*u(t)
%   Total solution: x(t) = e^(A t)*x0 - (B/A)*u_m*(1 - e^(A t))
%
% Plots:
%   (1) Free and forced components
%   (2) Complete solution
% ======================================================

clear; close all; clc;

% ---- Parameters ----
A = -1.0;     % stable system (A < 0)
B = 1.0;
x0 = 1.0;
u_m = 2.0;    % chosen so that -(B/A)*u_m = 2

% ---- Time and analytic responses ----
t = linspace(0, 5, 400);
x_free   = x0 * exp(A*t);
x_forced = -(B/A)*u_m*(1 - exp(A*t));
x_total  = x_free + x_forced;

% ---- Figure setup ----
figure('Color','w','Units','normalized','Position',[0.2 0.15 0.45 0.65]);

tiledlayout(2,1,'Padding','compact','TileSpacing','compact');

% ======================================================
% (1) Free + Forced components
% ======================================================
nexttile; hold on;
plot(t, x_free, '--', 'Color', [0 0.447 0.741], 'LineWidth', 1.5, ...
    'DisplayName', 'free response ($u\equiv0$)');
plot(t, x_forced, '--', 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 1.5, ...
    'DisplayName', 'forced response ($x_0=0$)');

yline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 0.8, 'HandleVisibility','off');
yline(2, 'Color', [0.5 0.5 0.5], 'LineStyle', ':', 'LineWidth', 0.8, 'HandleVisibility','off');

title('Decomposition of the total response', 'Interpreter', 'latex');
xlabel('');
ylabel('$x(t)$', 'Interpreter', 'latex');
grid on; box on;
legend('Interpreter','latex','FontSize',10,'Location','best','Box','off');

xlim([t(1) t(end)]);
yl = [min([x_free x_forced]) max([x_free x_forced])];
dy = diff(yl);
ylim([yl(1)-0.05*dy, yl(2)+0.05*dy]);  % padding

% ======================================================
% (2) Complete solution
% ======================================================
nexttile; hold on;
plot(t, x_total, 'Color', [0.8500 0.3250 0.0980], 'LineWidth', 1.8, ...
    'DisplayName', 'complete solution');

yline(2, 'Color', [0.5 0.5 0.5], 'LineStyle', ':', 'LineWidth', 0.8, 'HandleVisibility','off');

title('Complete Solution: $x(t) = 2 - e^{A t}$', 'Interpreter', 'latex');
xlabel('Time $t$', 'Interpreter', 'latex');
ylabel('$x(t)$', 'Interpreter', 'latex');
legend('Interpreter','latex','FontSize',10,'Location','best','Box','off');
grid on; box on;

xlim([t(1) t(end)]);
yl = [min(x_total) max(x_total)];
dy = diff(yl);
ylim([yl(1)-0.05*dy, yl(2)+0.05*dy]);  % padding

hold off;

Overwriting decomposition_first_order.m

6. Effect of the variables time constants and stead-state gains on cruise control example

%%writefile first_order_error_model.m
% ======================================================
% Demonstrate the impact of time constant (T = 1/(-A))
% and steady-state gain (-B/A) on the first-order error model:
%
%     e_dot = A * e + B * v_d(t)
%
% Desired speed:
%     v_d(t) = 40 for t < 5
%            = 50 for t >= 5
% ======================================================

clear; close all; clc;

% ======================================================
% Simulation parameters
% ======================================================
t_span = [0, 40];
t_eval = linspace(t_span(1), t_span(2), 1200);
x0_vec = 0.0;           % single IC to isolate parameter effects

A_list = [-0.25, -0.5, -1.0, -2.0];   % affects time constant T = 1/(-A)
B_list = [0.05, 0.1, 0.2, 0.4];       % affects steady-state gain -B/A

A_fixed = -1.0;
B_fixed = 0.1;

% ======================================================
% Desired speed profile
% ======================================================
v_d = @(t) (t < 5)*40 + (t >= 5)*50;

% ======================================================
% Helper: ODE for given A, B
% ======================================================
make_e_dot = @(A, B) @(t, e) A*e + B*v_d(t);

% ======================================================
% (1) Plot v_d(t)
% ======================================================
figure('Color','w','Units','normalized','Position',[0.25 0.4 0.45 0.35]);
plot(t_eval, v_d(t_eval), 'k', 'LineWidth', 2);
xline(5.0, '--', 'Color', 'k', 'LineWidth', 1);
text(5.0, 41, 'step @ 5 s', 'HorizontalAlignment','left', ...
    'VerticalAlignment','bottom', 'FontSize', 10);
title('Desired Speed $v_d(t)$', 'Interpreter','latex');
xlabel('Time [s]', 'Interpreter','latex');
ylabel('Desired Speed [m/s]', 'Interpreter','latex');
grid on; box on;
ax = gca; ax.GridAlpha = 0.35;
yl = ylim; dy = diff(yl);
ylim([yl(1)-0.05*dy, yl(2)+0.05*dy]);  % padding

% ======================================================
% (2) Vary A (time constant effect), B fixed
% ======================================================
figure('Color','w','Units','normalized','Position',[0.15 0.3 0.65 0.45]);
cmap = lines(numel(A_list));

for i = 1:length(A_list)
    A = A_list(i);
    if A >= 0
        error('Use A < 0 for stability.');
    end
    T = 1 / (-A);
    color = cmap(i,:);
    e_dot = make_e_dot(A, B_fixed);
    [t_sol, e_sol] = ode45(e_dot, t_span, x0_vec);
    e_interp = interp1(t_sol, e_sol, t_eval, 'linear');
    label = sprintf('$A=%.2f\\ (T=1/(-A)=%.2f\\,\\mathrm{s})$', A, T);
    plot(t_eval, e_interp, 'Color', color, 'LineWidth', 1.8, ...
        'DisplayName', label); hold on;
end

title('Error $e(t)$ for different $A$ (time constant effect), $B$ fixed', ...
    'Interpreter','latex');
xlabel('Time [s]', 'Interpreter','latex');
ylabel('$e(t)$ [m/s]', 'Interpreter','latex');
grid on; box on;
ax = gca; ax.GridAlpha = 0.35;
legend('Interpreter','latex','FontSize',9,'Location','eastoutside');
yl = ylim; dy = diff(yl);
ylim([yl(1)-0.05*dy, yl(2)+0.05*dy]);  % padding

% ======================================================
% (3) Vary B (steady-state gain effect), A fixed
% ======================================================
figure('Color','w','Units','normalized','Position',[0.15 0.3 0.65 0.45]);
cmap = lines(numel(B_list));

for i = 1:length(B_list)
    B = B_list(i);
    color = cmap(i,:);
    e_dot = make_e_dot(A_fixed, B);
    [t_sol, e_sol] = ode45(e_dot, t_span, x0_vec);
    e_interp = interp1(t_sol, e_sol, t_eval, 'linear');
    label = sprintf('$B=%.2f,\\ -B/A=%.2f$', B, -B/A_fixed);
    plot(t_eval, e_interp, 'Color', color, 'LineWidth', 1.8, ...
        'DisplayName', label); hold on;

    % steady-state level (v_d → 50)
    e_ss = -(B / A_fixed) * 50.0;
    yline(e_ss, ':', 'Color', color, 'LineWidth', 1, ...
        'HandleVisibility','off');
end

title('Error $e(t)$ for different $B$ (steady-state gain effect), $A$ fixed', ...
    'Interpreter','latex');
xlabel('Time [s]', 'Interpreter','latex');
ylabel('$e(t)$ [m/s]', 'Interpreter','latex');
grid on; box on;
ax = gca; ax.GridAlpha = 0.35;
legend('Interpreter','latex','FontSize',9,'Location','eastoutside');
yl = ylim; dy = diff(yl);
ylim([yl(1)-0.05*dy, yl(2)+0.05*dy]);  % padding

hold off;

Overwriting first_order_error_model.m

7. Complex number plot

%%writefile plot_H.m
% ======================================================
% Plot a complex number H on the complex plane (Python-style)
% ======================================================
function plot_H(H, show_degrees, save_path)
    if nargin < 2, show_degrees = false; end
    if nargin < 3, save_path = ''; end

    % Extract components
    a = real(H);
    b = imag(H);
    r = abs(H);
    theta = angle(H);

    % --- Figure setup ---
    figure('Position', [100, 100, 600, 600], 'Color', 'w');
    hold on; grid on;
    ax = gca;
    ax.Box = 'on';
    ax.XColor = 'k';
    ax.YColor = 'k';
    ax.LineWidth = 1.2;
    axis equal;

    % --- Plot main vector ---
    plot([0, a], [0, b], 'Color', [0, 0, 0.5], 'LineWidth', 2.5);
    scatter(a, b, 50, 'filled', 'MarkerFaceColor', [0, 0, 0.5]);

    % --- Projections ---
    plot([a, a], [0, b], 'k--', 'LineWidth', 1);
    plot([0, a], [b, b], 'k--', 'LineWidth', 1);

    % --- Labels ---
    text(a/2, -0.3, '$a$', 'Interpreter', 'latex', 'FontSize', 12, ...
         'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
    text(0.15, b/2, '$b$', 'Interpreter', 'latex', 'FontSize', 12, ...
         'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle');
    text(a*0.6, b*0.55, '$r$', 'Interpreter', 'latex', 'FontSize', 12, ...
         'Color', [0, 0, 0.5]);
    text(a*1.05, b*1.05, '$H$', 'Interpreter', 'latex', 'FontSize', 14, ...
         'Color', [0, 0, 0.5]);

    % --- Angle arc ---
    arc_radius = 0.3 * max([1, abs(a), abs(b)]);
    theta_deg = rad2deg(theta);
    theta_arc = linspace(0, theta, 100);
    plot(arc_radius*cos(theta_arc), arc_radius*sin(theta_arc), ...
         'Color', [0.5 0.5 0.5], 'LineWidth', 1.5);

    label_angle = theta / 2; if label_angle == 0, label_angle = 1e-4; end
    if show_degrees
        angle_label = sprintf('$\\theta = %.1f^\\circ$', theta_deg);
    else
        angle_label = sprintf('$\\theta = %.2f\\,\\text{rad}$', theta);
    end
    text(arc_radius*1.2*cos(label_angle), arc_radius*1.2*sin(label_angle), ...
         angle_label, 'Interpreter', 'latex', 'FontSize', 12);

    % --- Axes limits ---
    lim = max([abs(a), abs(b), r]) * 1.3;
    xlim([-lim, lim]);
    ylim([-lim, lim]);

    % --- Full-length axes lines (draw *after* limits are set) ---
    plot([-lim lim], [0 0], 'k', 'LineWidth', 1.5); % x-axis
    plot([0 0], [-lim lim], 'k', 'LineWidth', 1.5); % y-axis

    % --- Axis labels (fixed placement like Python) ---
    text(0, -lim*1.1, 'Real axis', 'FontSize', 12, ...
         'HorizontalAlignment', 'center', 'VerticalAlignment', 'top');
    text(-lim*1.1, 0, 'Imag axis', 'FontSize', 12, ...
         'Rotation', 90, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');

    % --- Title ---
    if show_degrees
        title(sprintf('$H = a + bj,\\; |H| = %.2f,\\; \\theta = %.2f^\\circ$', ...
              r, theta_deg), 'Interpreter', 'latex');
    else
        title(sprintf('$H = a + bj,\\; |H| = %.2f,\\; \\theta = %.2f\\,\\text{rad}$', ...
              r, theta), 'Interpreter', 'latex');
    end

    % --- Save or show ---
    if ~isempty(save_path)
        saveas(gcf, save_path);
        fprintf('Saved plot as %s\n', save_path);
    end
end

% Example usage:
plot_H(3 + 4j, true);

Overwriting plot_H.m

8. Sine Response

%%writefile sine_response.m
% ======================================================
% Time-domain response to sinusoidal input (matches Python output)
% ======================================================

% Define state-space system
A = -1;
B = 1;
C = 2;
D = 0;
sys_ss = ss(A, B, C, D);

% Simulation parameters
omega = 3.0;
t = linspace(0.0, 20.0, 1000);
u = sin(omega * t);
x0 = 3.0;  % initial condition on the state

% Simulate
[y, t_out, x] = lsim(sys_ss, u, t, x0);

% --- Plot ---
figure;
plot(t_out, u, 'b', 'LineWidth', 1.2); hold on;
plot(t_out, y, 'r', 'LineWidth', 1.5);
xlabel('Time [s]');
ylabel('Amplitude');
title('Input and Output');
legend('u(t) = sin(3t)', 'y(t)', 'Location', 'best');
grid on;

% Add space around limits
xlim([min(t)*0.95, max(t)*1.05]);
ylim padded;

Overwriting sine_response.m

9. Bode Plot

%%writefile bode_comparison.m
% ======================================================
% Frequency response comparison: numeric vs analytic
% ======================================================

% Define transfer function
s = tf('s');
sys_tf = 2 / (s + 1);

% Frequency grid
w = logspace(-2, 2, 800);   % rad/s
[H, wout] = freqresp(sys_tf, w);   % complex frequency response

% Extract magnitude and phase
mag = abs(squeeze(H));
phase = angle(squeeze(H));

% Target frequency
omega = 3;

% Analytical reference values
mag_ref = 2.0 / sqrt(1.0 + omega^2);
phase_ref = -atan(omega);

% Numeric interpolation (for plotting)
mag_at_3 = interp1(w, mag, omega);
phase_at_3 = interp1(w, unwrap(phase), omega);

fprintf('|G(j3)|  numeric = %.6f,  analytic = %.6f\n', mag_at_3, mag_ref);
fprintf('∠G(j3)  numeric = %.4f°,  analytic = %.4f°\n', rad2deg(phase_at_3), rad2deg(phase_ref));

% --- Plot magnitude ---
figure;
loglog(w, mag, 'b', 'LineWidth', 1.5); hold on;
scatter(omega, mag_ref, 50, 'r', 'filled');
xline(omega, '--r', 'LineWidth', 1.0);
yline(mag_ref, '--r', 'LineWidth', 1.0);
text(omega*1.05, mag_ref*1.1, sprintf('|G(j3)| ≈ %.3f', mag_ref), 'Color', 'r');

xlabel('Frequency \omega [rad/s]');
ylabel('|G(j\omega)|');
title('Bode Magnitude');
legend('|G(j\omega)|', 'Location', 'best');
grid on; grid minor;

% Add space around limits
xlim([min(w)*0.8, max(w)*1.2]);
ylim([min(mag)*0.8, max(mag)*1.2]);

% --- Plot phase ---
figure;
semilogx(w, rad2deg(phase), 'b', 'LineWidth', 1.5); hold on;
scatter(omega, rad2deg(phase_ref), 50, 'r', 'filled');
xline(omega, '--r', 'LineWidth', 1.0);
yline(rad2deg(phase_ref), '--r', 'LineWidth', 1.0);
text(omega*1.05, rad2deg(phase_ref)+5, sprintf('%.1f°', rad2deg(phase_ref)), 'Color', 'r');

xlabel('Frequency \omega [rad/s]');
ylabel('Phase [deg]');
title('Bode Phase');
legend('\angle G(j\omega)', 'Location', 'best');
grid on; grid minor;

% Add space around limits
xlim([min(w)*0.8, max(w)*1.2]);
ylim padded;

Overwriting bode_comparison.m