MATLAB

MATLAB is the chemical engineer's numerical workbench. Learn to write mole and energy balances as ODEs, solve them with ode45/ode15s (stiff systems), fit kinetic parameters to data, solve nonlinear equations with fzero/fsolve, and simulate control loops with Simulink — with complete scripts.

Most chemical engineering models are differential equations: concentrations changing with time in a batch reactor, with length in a PFR, temperatures in a heat exchanger. MATLAB solves them in a few lines and plots the results instantly.

Python (SciPy) offers the same tools and is increasingly popular; the ideas transfer directly.

An analogy

MATLAB is a super-calculator that understands equations of change: tell it how things change, and it tells you what happens.

Scripts, functions and plots

Vectors and matrices are native; element-wise operators (.*, ./, .^) apply formulas to whole arrays. Anonymous functions (@(x) …) define models in one line; plot, xlabel and legend make publication-ready figures.

ODE solvers

ode45 (explicit Runge–Kutta) for non-stiff problems; ode15s/ode23s for stiff problems where time scales differ greatly. Events stop integration at a condition (e.g. 95 % conversion). PFRs use volume instead of time as the independent variable.

Algebraic equations and optimisation

fzero (one equation — e.g. CSTR steady states, bubble points), fsolve (systems — e.g. multicomponent flash), fminsearch/fmincon (optimisation — e.g. optimal reactor temperature), lsqcurvefit (parameter fitting).

Simulink and toolboxes

Simulink models dynamic systems as block diagrams — ideal for control loops (PID blocks, transfer functions, disturbances). Control System Toolbox (step, bode, pidtune), Optimization and Global Optimization toolboxes, and Model Predictive Control Toolbox support advanced process control.

Formulas

Worked example: Optimum batch time for intermediate

Given: A → B → C, k₁ = 0.5, k₂ = 0.2 min⁻¹, A₀ = 1 mol/L

Find: t_max and B_max.

Answer: 3.05 min, 0.543 mol/L (ode45 agrees)

Worked example: Time for 90 % conversion

Given: Fitted k = 0.100 min⁻¹

Find: t for X = 0.9.

Answer: 23 min

Series reactions in a batch reactor (ode45)

% A -> B -> C, first-order steps in an isothermal batch reactor
k1 = 0.5;  k2 = 0.2;            % 1/min
C0 = [1; 0; 0];                 % initial A, B, C (mol/L)

rhs = @(t, C) [ -k1*C(1);
                 k1*C(1) - k2*C(2);
                 k2*C(2) ];

[t, C] = ode45(rhs, [0 20], C0);

[Bmax, i] = max(C(:,2));
fprintf('B_max = %.3f mol/L at t = %.2f min\n', Bmax, t(i));
% -> B_max = 0.543 mol/L at t = 3.05 min

plot(t, C, 'LineWidth', 1.5); grid on
legend('A', 'B', 'C'); xlabel('t (min)'); ylabel('C (mol/L)')

Write the mole balances as a function returning dC/dt, pick a solver and time span, and MATLAB integrates. For stiff systems (fast and slow reactions, or reactors with energy balances), switch to ode15s.

Fit a rate constant to data (lsqcurvefit)

% Batch data: first-order decay C = C0*exp(-k t)
t    = [0 5 10 15 20 30];              % min
Cexp = [2.00 1.21 0.74 0.45 0.27 0.10]; % mol/L

model = @(p, t) p(1) * exp(-p(2) * t);  % p = [C0, k]
p0    = [2, 0.05];
p     = lsqcurvefit(model, p0, t, Cexp);

fprintf('C0 = %.2f mol/L, k = %.3f 1/min\n', p(1), p(2));
% -> k ≈ 0.100 1/min, so 90 % conversion needs ln(10)/k = 23 min

Nonlinear regression finds kinetic parameters directly from concentration–time data, avoiding the distortion of linearised plots. Report confidence intervals (nlparci) too.

Common mistakes

Applications

Key takeaways

Going further