Dynamic Simulink ODE Solver Generator

A MATLAB script dynamically generates Simulink models for nth-order ODEs with time-varying coefficients, automating model creation.

1. Architecture: The Model-Data Separation

This project provides a MATLAB script suite centered around a powerful Model-Data Separation architecture for solving any n-th order linear Ordinary Differential Equation (ODE).

This architecture consists of two main components:

1. generate_variable_coeff_model.m (The Model Generator - The Architect):

• This script acts as the "architect." Its sole job is to dynamically generate a structured Simulink model framework (.slx file).
• This generated framework is like a powerful, unconfigured calculator. It knows how to perform integration, multiplication, addition, and division, but it contains no specific numerical values or functions itself. It only defines the computational workflow.

1. test_ode_solver.m (The Simulation Script - The Operator):

• This script acts as the "operator." Its job is to define a specific mathematical problem.
• It defines all the equation's coefficients A_i(t) and the input term u(t) by creating a series of time-based vectors.
• It then "injects" this specific data into the generic Simulink model, drives the simulation, and retrieves the solution.

The advantage of this separation is immense: you can use the same, unchanged model generator to solve an infinite variety of different ODEs simply by defining new data in the simulation script.

2. Capabilities & Engineering Applications

The core script, generate_variable_coeff_model.m, is capable of uniformly handling all of the following types of n-th order linear ODEs:

An(t)·x⁽ⁿ⁾(t) + ... + A1(t)·x'(t) + A0(t)·x(t) = u(t)

• N-th Order: Solve second-order, third-order, or even higher-order equations.
• Variable Coefficients: Supports coefficients that are functions of time, such as sin(t), t^2+1, etc.
• Constant Coefficients: Handled as a special case of variable coefficients (a function that does not change with time).
• Non-homogeneous/Homogeneous: The forcing term u(t) can be any function or zero (for the homogeneous case).

This versatility allows it to solve a vast range of engineering problems, from basic RLC circuits and mechanical vibrations to complex simulations like aerospace vehicle dynamics (with varying mass and aerodynamic coefficients) and verification of advanced gain-scheduling control systems.

3. How to Use: Illustrated with Engineering Examples

Below are three examples demonstrating this process for problems in Biology, Physics, and Chemistry.

Example 1: Biology

Problem: Modeling drug concentration in the bloodstream. A drug is administered via an IV drip, and the body's elimination rate decreases over time.

Equation (1st Order): 1 * C'(t) + k(t) * C(t) = D(t)

• C(t): Drug concentration (mg/L).
• k(t) = 0.5e^{-0.1t} + 0.1: Time-varying elimination rate.
• D(t): Dosing rate (mg/L per hr).
• n=1, A1(t)=1, A0(t)=k(t), u(t)=D(t).
• Initial Condition: C(0) = 0.

run_biology_example.m

% Copyright (c) 2025, Jiarui Liang
% All rights reserved.
% Organization: [Macau University of Science and Technology]
% Email: [[e-mail]]
%
% THIS SCRIPT IS PART OF THE DYNAMIC SIMULINK ODE SOLVER PROJECT
% -----------------------------------------------------------------

%% --- Step 1: Generate the variable-coefficient model ---
n = 1;
initial_conditions = [0]; % C(0) = 0
model_name = 'Biology_DrugModel';
generate_variable_coeff_model(n, initial_conditions, model_name);

%% --- Step 2: Define ALL input signals for the simulation ---
% Define the time vector
t = (0:0.1:24)'; % Simulate for 24 hours

% Create the signal for u(t) (IV drip on for 8 hours)
u_signal = 5 * (t <= 8);

% Create signals for each coefficient A_i(t)
A1_signal = ones(size(t));
A0_signal = 0.5 * exp(-0.1 * t) + 0.1; % k(t)

%% --- Step 3: Create a SimulationInput object and run the simulation ---
% The order of columns must match the Inport block creation order:
% u_t, A1_t, A0_t
input_data = [t, u_signal, A1_signal, A0_signal];

simIn = Simulink.SimulationInput(model_name);
simIn = simIn.setExternalInput(input_data);
simIn = simIn.setModelParameter('StopTime', '24');

disp('Running Biology simulation...');
simOut = sim(simIn);
disp('Simulation finished.');

%% --- Step 4: Plot the results ---
figure;
plot(simOut.yout{1}.Values.Time, simOut.yout{1}.Values.Data, 'b-', 'LineWidth', 2);
title('Biology: Drug Concentration in Bloodstream');
xlabel('Time (hours)');
ylabel('Concentration C(t) (mg/L)');
grid on;
legend('Drug Concentration');

Example 2: Physics

Problem: Modeling a rocket launch with decreasing mass as fuel is burned.

Equation (2nd Order): m(t) * y''(t) + c * y'(t) + k * y(t) = F(t)

• y(t): Altitude (m).
• m(t) = 2000 - 50t: Time-varying mass (kg).
• c = 100, k = 0.
• F(t): Engine thrust (N).
• n=2, A2(t)=m(t), A1(t)=c, A0(t)=k, u(t)=F(t).
• Initial Conditions: y'(0) = 0, y(0) = 0.

run_physics_example.m

% Copyright (c) 2025, Jiarui Liang
% All rights reserved.
% Organization: [Macau University of Science and Technology]
% Email: [[e-mail]]
%
% THIS SCRIPT IS PART OF THE DYNAMIC SIMULINK ODE SOLVER PROJECT
% -----------------------------------------------------------------

%% --- Step 1: Generate the variable-coefficient model ---
n = 2;
initial_conditions = [0, 0]; % [y'(0), y(0)]
model_name = 'Physics_RocketModel';
generate_variable_coeff_model(n, initial_conditions, model_name);

%% --- Step 2: Define ALL input signals for the simulation ---
t = (0:0.05:30)'; % Simulate for 30 seconds

% u(t) is thrust F(t), on for 20s
burn_duration = 20;
u_signal = 35000 * (t <= burn_duration);

% A2(t) is mass m(t)
m0 = 2000;
burn_rate = 50;
mass_after_burn = m0 - burn_rate * burn_duration;
A2_signal = (m0 - burn_rate * t) .* (t <= burn_duration) + mass_after_burn * (t > burn_duration);

% A1(t) is damping c, A0(t) is spring k
A1_signal = 100 * ones(size(t));
A0_signal = zeros(size(t));

%% --- Step 3: Create a SimulationInput object and run the simulation ---
% Order: u_t, A2_t, A1_t, A0_t
input_data = [t, u_signal, A2_signal, A1_signal, A0_signal];

simIn = Simulink.SimulationInput(model_name);
simIn = simIn.setExternalInput(input_data);
simIn = simIn.setModelParameter('StopTime', '30');

disp('Running Physics simulation...');
simOut = sim(simIn);
disp('Simulation finished.');

%% --- Step 4: Plot the results ---
figure;
plot(simOut.yout{1}.Values.Time, simOut.yout{1}.Values.Data, 'r-', 'LineWidth', 2);
title('Physics: Rocket Altitude with Variable Mass');
xlabel('Time (seconds)');
ylabel('Altitude y(t) (meters)');
grid on;
legend('Rocket Altitude');

Example 3: Chemistry

Problem: Modeling a contaminant in a reactor where processes follow a daily cycle.

Equation (3rd Order): 1*C''' + A2(t)*C'' + A1(t)*C' + A0(t)*C = S(t)

• C(t): Contaminant concentration.
• A2(t), A1(t): Cyclical diffusion/dissipation coefficients.
• A0(t): Constant reaction rate.
• S(t): Source term from a spill.
• n=3, A3=1, A2(t), A1(t), A0(t).
• Initial Conditions: C''(0)=0, C'(0)=0, C(0)=0.

run_chemistry_example.m

% Copyright (c) 2025, Jiarui Liang
% All rights reserved.
% Organization: [Macau University of Science and Technology]
% Email: [[e-mail]]
%
% THIS SCRIPT IS PART OF THE DYNAMIC SIMULINK ODE SOLVER PROJECT
% -----------------------------------------------------------------

%% --- Step 1: Generate the variable-coefficient model ---
n = 3;
initial_conditions = [0, 0, 0]; % [C''(0), C'(0), C(0)]
model_name = 'Chemistry_ReactorModel';
generate_variable_coeff_model(n, initial_conditions, model_name);

%% --- Step 2: Define ALL input signals for the simulation ---
t = (0:0.1:72)'; % Simulate for 72 hours

% u(t) is the source term S(t) from a spill
spill_time = 5;
u_signal = (t >= spill_time) .* 50 .* exp(-0.5 * (t - spill_time));

% Coefficients for A3*C''' + A2*C'' + A1*C' + A0*C
A3_signal = ones(size(t));
A2_signal = 2 + cos(2 * pi * t / 24);
A1_signal = 5 * (1 + 0.5 * sin(2*pi*t/24));
A0_signal = 10 * ones(size(t));

%% --- Step 3: Create a SimulationInput object and run the simulation ---
% Order: u_t, A3_t, A2_t, A1_t, A0_t
input_data = [t, u_signal, A3_signal, A2_signal, A1_signal, A0_signal];

simIn = Simulink.SimulationInput(model_name);
simIn = simIn.setExternalInput(input_data);
simIn = simIn.setModelParameter('StopTime', '72');

disp('Running Chemistry simulation...');
simOut = sim(simIn);
disp('Simulation finished.');

%% --- Step 4: Plot the results ---
figure;
plot(simOut.yout{1}.Values.Time, simOut.yout{1}.Values.Data, 'g-', 'LineWidth', 2);
title('Chemistry: Contaminant Concentration in a Cyclical Reactor');
xlabel('Time (hours)');
ylabel('Concentration C(t)');
grid on;
legend('Contaminant Concentration');

4. Important Notes and Best Practices

• Singularity Warning: The highest-order coefficient An(t) must never be zero within the simulation time domain. This would cause a "division by zero" error, which reflects a mathematical singularity in the equation itself, not a bug in the code.
• Dimensional Consistency: Ensure all signal vectors (u_signal, A_i_signal) have the exact same number of rows as the time vector t. For constants, use C * ones(size(t)) or zeros(size(t)).
• Simulation Parameters: The recommended way to set parameters like stop time or solver is programmatically via simIn = simIn.setModelParameter('ParameterName', 'value').
• Data Extraction: To access simulation results from the output object, use the format simOut.yout{1}.Values, which contains the .Time and .Data properties.

5. Acknowledgements

This project builds upon the foundational work of Nobel Mondal's "Automated Simulink Model Creator from Ordinary Differential Equation", a popular entry on the MATLAB File Exchange.

While Mondal's tool established the powerful concept of programmatic model generation for constant-coefficient systems, this project extends that concept to solve ODEs with time-varying coefficients. Our implementation is therefore a complementary tool designed to address a different and more complex class of dynamic systems.

The original File Exchange submission can be found here: https://www.mathworks.com/matlabcentral/fileexchange/53160-automated-simulink-model-creator-from-ordinary-differential-equation

License

This source code is licensed under the MIT License. See the header of the generate_variable_coeff_model.m file for full details.

Copyright (c) 2025, Jiarui Liang

All rights reserved.

Macau University of Science and Technology

Email: [e-mail]