diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 009cfcf3e39..a2b14c60f3a 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -66,10 +66,25 @@ jobs:
           echo "LD_PRELOAD=$LIB_STDCXX:$LIB_OPENBLAS:$LIB_LAPACK" >> $GITHUB_ENV
       - name: Set up MATLAB
         uses: matlab-actions/setup-matlab@718d4320188c73c86eb94ce76b553cbf89778487 # v2.5.0
+      - name: Set MATLAB search paths for tests and build MLTBX toolbox # v2.2.1
+        uses: matlab-actions/run-command@67f012f1ee4bc1627a7a95a0ccdeec745c9f6f36
+        with:
+          command: |
+            toolboxPath = fullfile(getenv('CANTERA_ROOT'));
+            libPath = fullfile(toolboxPath, 'build', 'lib');
+            includePath = fullfile(toolboxPath, 'include');
+            addpath(fullfile(toolboxPath, 'interfaces', 'matlab_experimental', 'Utility'));
+            ctPaths(libPath, includePath, toolboxPath);
       - name: Run tests
         uses: matlab-actions/run-tests@f9fc3d8ca29fadef6227fa52884b144b9011fa2f # v2.1.1
         with:
           select-by-folder: /home/runner/work/cantera/cantera/test/matlab_experimental
+      - name: Run MATLAB test_examples
+        uses: matlab-actions/run-command@67f012f1ee4bc1627a7a95a0ccdeec745c9f6f36 # v2.2.1
+        with:
+          command: |
+            addpath('/home/runner/work/cantera/cantera/samples/matlab_experimental');
+            test_examples;
 
   ubuntu-multiple-pythons:
     name: ${{ matrix.os }} with Python ${{ matrix.python-version }}, Numpy ${{ matrix.numpy || 'latest' }}, Cython ${{ matrix.cython || 'latest' }}
diff --git a/interfaces/matlab_experimental/Base/ThermoPhase.m b/interfaces/matlab_experimental/Base/ThermoPhase.m
index 8f42f53db8d..78b7cfc1667 100644
--- a/interfaces/matlab_experimental/Base/ThermoPhase.m
+++ b/interfaces/matlab_experimental/Base/ThermoPhase.m
@@ -66,6 +66,8 @@
         %     Scalar double mean molecular weight. Units: kg/kmol.
         meanMolecularWeight
 
+        massDensity % Mass basis density. Units: kg/m^3.
+
         molarDensity % Molar basis density. Units: kmol/m^3.
 
         molecularWeights % Molecular weights of the species. Units: kg/kmol.
@@ -903,6 +905,10 @@ function display(tp)
             mmw = ctFunc('thermo_meanMolecularWeight', tp.tpID);
         end
 
+        function density = get.massDensity(tp)
+            density = ctFunc('thermo_density', tp.tpID);
+        end
+
         function density = get.molarDensity(tp)
             density = ctFunc('thermo_molarDensity', tp.tpID);
         end
diff --git a/interfaces/matlab_experimental/OneDim/Sim1D.m b/interfaces/matlab_experimental/OneDim/Sim1D.m
index b506f736b8b..83207ad0c27 100644
--- a/interfaces/matlab_experimental/OneDim/Sim1D.m
+++ b/interfaces/matlab_experimental/OneDim/Sim1D.m
@@ -54,14 +54,9 @@ function delete(s)
 
         %% Sim1D Utility Methods
 
-        function display(s, fname)
+        function display(s)
             % Show all domains.
-
-            if nargin == 1
-                fname = '-';
-            end
-
-            ctFunc('sim1D_show', s.stID, fname);
+            ctFunc('sim1D_show', s.stID);
         end
 
         function restore(s, fname, id)
diff --git a/interfaces/matlab_experimental/Utility/ctLoad.m b/interfaces/matlab_experimental/Utility/ctLoad.m
index cdedd769c0c..783910292b5 100644
--- a/interfaces/matlab_experimental/Utility/ctLoad.m
+++ b/interfaces/matlab_experimental/Utility/ctLoad.m
@@ -1,32 +1,39 @@
 function ctLoad()
+    % ctLoad
     % Load the Cantera C Library into Memory
 
+    paths = ctPaths();
+
+    if any(cellfun(@isempty, {paths.libPath, paths.includePath, paths.toolboxPath}))
+        error('ctLoad:MissingPath', ...
+              'Library, header, and toolbox paths must be configured with ctPaths(libPath,includePath, toolboxPath).');
+    end
+
     if ispc
-        ctName = '/bin/cantera_shared.dll';
+        libName = 'cantera_shared.dll';
     elseif ismac
-        ctName = '/Lib/libcantera_shared.dylib';
+        libName = 'libcantera_shared.dylib';
     elseif isunix
-        ctName = '/lib/libcantera_shared.so';
+        libName = 'libcantera_shared.so';
     else
-        error('Operating System Not Supported!');
-        return;
+        error('ctLoad:UnsupportedPlatform', 'Operating system not supported.');
     end
 
-    root = ctRoot;
-    if ispc
-        root = [ctRoot, '/Library'];
-    end
+    fullLibPath = fullfile(paths.libPath, libName);
+    fullHeaderPath = fullfile(paths.includePath, 'cantera', 'clib', 'ctmatlab.h');
 
     if ~libisloaded(ctLib)
-        [~, warnings] = loadlibrary([root, ctName], ...
-                                    [root, '/include/cantera/clib/ctmatlab.h'], ...
-                                    'includepath', [root, '/include'], ...
-                                    'addheader', 'ct', 'addheader', 'ctfunc', ...
-                                    'addheader', 'ctmultiphase', 'addheader', ...
-                                    'ctonedim', 'addheader', 'ctreactor', ...
-                                    'addheader', 'ctrpath', 'addheader', 'ctsurf');
+        [~, warnings] = loadlibrary(fullLibPath, fullHeaderPath, ...
+                                    'includepath', paths.includePath, ...
+                                    'addheader', 'ct', ...
+                                    'addheader', 'ctfunc', ...
+                                    'addheader', 'ctmultiphase', ...
+                                    'addheader', 'ctonedim', ...
+                                    'addheader', 'ctreactor', ...
+                                    'addheader', 'ctrpath', ...
+                                    'addheader', 'ctsurf');
     end
 
-    disp(sprintf('Cantera %s is ready for use.', ctVersion))
+    fprintf('Cantera %s is ready for use.\n', ctVersion);
 
 end
diff --git a/interfaces/matlab_experimental/Utility/ctPaths.m b/interfaces/matlab_experimental/Utility/ctPaths.m
new file mode 100644
index 00000000000..d5d95c795fc
--- /dev/null
+++ b/interfaces/matlab_experimental/Utility/ctPaths.m
@@ -0,0 +1,73 @@
+function paths = ctPaths(varargin)
+    % ctPaths ::
+    % Configure or retrieve the library/header/toolbox paths for Cantera.
+    % The paths are stored as MATLAB preferences.
+    %
+    %   >> ctPaths()                       % Get current config as struct
+    %   >> ctPaths(libPath, includePath)   % Set library/header/toolbox paths
+    %   >> ctPaths('clear')                % Clear saved paths
+    %
+    % :return:
+    %   paths: struct with fields 'libPath', 'includePath', and
+    %   'toolboxPath'
+
+    if nargin == 1 && strcmp(varargin{1}, 'clear')
+        if ispref('Cantera', 'Paths')
+            paths = getpref('Cantera', 'Paths');
+            subDirs = strsplit(genpath(paths.toolboxPath), pathsep);
+            currentPaths = strsplit(path, pathsep);
+            subdirsToRemove = intersect(subDirs, currentPaths);
+            if ~isempty(subdirsToRemove)
+                rmpath(subdirsToRemove{:});
+            end
+            rmpref('Cantera', 'Paths');
+        end
+        paths = struct('libPath', '', 'includePath', '', 'toolboxPath', '');
+        return
+    elseif nargin == 3 && all(cellfun(@ischar, varargin))
+        paths = struct('libPath', varargin{1}, ...
+                       'includePath', varargin{2}, ...
+                       'toolboxPath', varargin{3});
+        setpref('Cantera', 'Paths', paths);
+    else
+        % Load from saved preferences if available
+        if ispref('Cantera', 'Paths')
+            paths = getpref('Cantera', 'Paths');
+            return
+        else
+            paths = struct('libPath', '', 'includePath', '', 'toolboxPath', '');
+        end
+    end
+
+    if any(cellfun(@isempty, {paths.libPath, paths.includePath, paths.toolboxPath}))
+        error('ctPaths:MissingPath', ...
+              'Library, header, and toolbox paths must be configured with ctPaths(libPath,includePath, toolboxPath).');
+    end
+
+    mapping = {
+        'interfaces/matlab_experimental', 'toolbox';
+        'samples/matlab_experimental', 'samples';
+        'test/matlab_experimental', 'test/matlab_toolbox';
+        'test/data', 'test/data';
+        'data', 'data'
+        };
+
+    % Check whether user is using Cantera source code or MLTBX based on folder structure
+    if isfolder(fullfile(paths.toolboxPath, 'interfaces'))
+        col = 1;
+    else
+        col = 2;
+    end
+
+    for i = 1:size(mapping, 1)
+        subdir = fullfile(paths.toolboxPath, mapping{i, col});
+        if isfolder(subdir)
+            addpath(genpath(subdir));
+        else
+            warning('ctPaths:MissingDirectory', ...
+                    'Directory not found: %s', subdir);
+        end
+    end
+
+    savepath();
+end
diff --git a/samples/matlab_experimental/conhp.m b/samples/matlab_experimental/conhp.m
index 8483a40d0a0..d4776825063 100644
--- a/samples/matlab_experimental/conhp.m
+++ b/samples/matlab_experimental/conhp.m
@@ -6,6 +6,7 @@
     % It assumes that the ``gas`` object represents a reacting ideal gas mixture.
 
     % Set the state of the gas, based on the current solution vector.
+    gas.basis = 'mass';
     gas.Y = y(2:end);
     gas.TP = {y(1), gas.P};
     nsp = gas.nSpecies;
@@ -13,7 +14,6 @@
     % energy equation
     wdot = gas.netProdRates;
     H = gas.partialMolarEnthalpies';
-    gas.basis = 'mass';
     tdot =- 1 / (gas.D * gas.cp) .* wdot * H;
 
     % set up column vector for dydt
diff --git a/samples/matlab_experimental/conuv.m b/samples/matlab_experimental/conuv.m
index 84f47687777..b183f25e5bf 100644
--- a/samples/matlab_experimental/conuv.m
+++ b/samples/matlab_experimental/conuv.m
@@ -6,6 +6,7 @@
     % It assumes that the ``gas`` object represents a reacting ideal gas mixture.
 
     % Set the state of the gas, based on the current solution vector.
+    gas.basis = 'mass';
     gas.Y = y(2:end);
     gas.TD = {y(1), gas.D};
     nsp = gas.nSpecies;
@@ -13,7 +14,6 @@
     % energy equation
     wdot = gas.netProdRates;
     U = gas.partialMolarIntEnergies';
-    gas.basis = 'mass';
     tdot =- 1 / (gas.D * gas.cv) .* wdot * U;
 
     % set up column vector for dydt
diff --git a/samples/matlab_experimental/crit_properites.m b/samples/matlab_experimental/crit_properites.m
new file mode 100644
index 00000000000..0b8e3d20044
--- /dev/null
+++ b/samples/matlab_experimental/crit_properites.m
@@ -0,0 +1,61 @@
+%% Critical state properties
+% Print the critical state properties for the fluids for which Cantera has
+% built-in liquid/vapor equations of state.
+%
+% .. tags:: Matlab, thermodynamics, multiphase, non-ideal fluid
+
+clear all;
+close all;
+
+tic
+help crit_properites
+
+%% Create a pure fluid object
+fluids = struct('water', Water(), ...
+    'nitrogen', Nitrogen(), ...
+    'methane', Methane(), ...
+    'hydrogen', Hydrogen(), ...
+    'oxygen', Oxygen(), ...
+    'carbon_dioxide', CarbonDioxide(), ...
+    'heptane', Heptane(), ...
+    'HFC_134a', HFC134a());
+
+names = fieldnames(fluids);
+
+%% Plot critical properties and print tabulated values
+
+figure;
+hold on;
+title('Critical Properties of Pure Fluids');
+xlabel('Critical Temperature [K]');
+ylabel('Critical Pressure [Pa]');
+xlim([0 750]);
+grid on;
+
+% Print header
+fprintf('Critical State Properties\n');
+fprintf('%-16s   %-7s   %-10s   %-7s\n', 'Fluid', 'Tc [K]', 'Pc [Pa]', 'Zc');
+fprintf('%s   %s   %s   %s\n', repmat('-',1,16), repmat('-',1,7), repmat('-',1,10), repmat('-',1,7));
+
+%% Loop through fluids
+for i = 1:length(names)
+    name = names{i};
+    f = fluids.(name);
+
+    Tc = f.critTemperature;
+    Pc = f.critPressure;
+    rhoc = f.critDensity;
+    mw = f.meanMolecularWeight;
+    R = GasConstant;
+
+    Zc = Pc * mw / (rhoc * R * Tc);
+
+    % Plot
+    plot(Tc, Pc, 'o');
+    text(Tc + 4, Pc + 2e5, strrep(name, '_', ' '), 'FontSize', 9);
+
+    % Print table row
+    fprintf('%-16s   %7.2f   %10.4g   %7.4f\n', strrep(name, '_', ' '), Tc, Pc, Zc);
+end
+
+toc
diff --git a/samples/matlab_experimental/diamond_cvd.m b/samples/matlab_experimental/diamond_cvd.m
index 651d4ec6f62..dd9499e5af2 100644
--- a/samples/matlab_experimental/diamond_cvd.m
+++ b/samples/matlab_experimental/diamond_cvd.m
@@ -69,7 +69,7 @@
     r = surf_phase.netProdRates;
     carbon_dot = r(iC);
     mdot = mw * carbon_dot;
-    rate = mdot / dbulk.D;
+    rate = mdot / dbulk.massDensity;
     xx = [xx; x(ih)];
     rr = [rr; rate * 1.0e6 * 3600.0];
     cov = [cov; surf_phase.coverages];
diff --git a/samples/matlab_experimental/diff_flame.m b/samples/matlab_experimental/diff_flame.m
index 021e54f0ef9..769fa42b3c0 100644
--- a/samples/matlab_experimental/diff_flame.m
+++ b/samples/matlab_experimental/diff_flame.m
@@ -104,7 +104,7 @@
 % ``help setRefineCriteria``.
 
 f.energyEnabled = true;
-fl.setRefineCriteria(2, 200.0, 0.1, 0.2);
+fl.setRefineCriteria(2, 4, 0.2, 0.3, 0.04);
 fl.solve(loglevel, refine_grid);
 
 %% Show statistics of solution and elapsed time
diff --git a/samples/matlab_experimental/equations_of_state.m b/samples/matlab_experimental/equations_of_state.m
new file mode 100644
index 00000000000..c29d1670ac2
--- /dev/null
+++ b/samples/matlab_experimental/equations_of_state.m
@@ -0,0 +1,227 @@
+%% Non-ideal equations of state
+%
+% This example demonstrates a comparison between ideal and non-ideal equations of
+% state (EoS) using Cantera.
+%
+% The following equations of state are used to evaluate thermodynamic properties
+% in this example:
+%
+% 1. Ideal-gas EoS from Cantera
+% 2. Non-ideal Redlich-Kwong EoS (R-K EoS) from Cantera
+%
+% .. tags:: Matlab, thermodynamics, non-ideal fluid
+
+clear all
+close all
+
+tic
+help equations_of_state
+
+%% Helper Functions
+% This examples uses \mathrm{ CO_2 } as the only species. The function
+% "get_thermo_Cantera" calculates thermodynamic properties based on the
+% thermodynamic state (T, p) of the species using Cantera. Applicable phases
+% are "Ideal-gas" and "Redlich-Kwong". The ideal-gas equation can be stated
+% as
+%
+% .. math:: 
+%     \mathrm{ pv = RT },
+%
+% where p, v, and T represent thermodynamic pressure, molar volume, and the
+% temperature of the gas-phase. R is the universal gas constant. The
+% Redlich-Kwong equation is a cubic, non-ideal equation of state,
+% represented as
+%
+% .. math::
+%     \mathrm{ p = \frac{RT}{v-b^\ast}-\frac{a^\ast}{v\sqrt{T}(v+b^\ast)} }.
+%
+% In this expression, R is the universal gas constant and v is the molar
+% volume. The temperature-dependent van der Waals attraction parameter
+% \mathrm{ a^\ast } and volume correction parameter (repulsive parameter)
+% \mathrm{ b^\ast } represent molecular interactions.
+%
+% To plot the comparision of thermodynamic properties among the tree EoS,
+% the "plotEoS" function is used.
+
+function [h, u, s, cp, cv] = get_thermo_Cantera(phase, T, p)
+    phase.basis = {'mass'};
+    n = length(p);
+    u = zeros(1, n);
+    h = zeros(1, n);
+    s = zeros(1, n);
+    cp = zeros(1, n);
+    cv = zeros(1, n);
+
+    for i = 1:n
+        phase.TPX = {T, p(i), 'CO2:1.0'};
+        u(i) = phase.U / 1000;
+        h(i) = phase.H / 1000;
+        s(i) = phase.S / 1000;
+        cp(i) = phase.cp / 1000;
+        cv(i) = phase.cv / 1000;
+    end
+
+    % Reference to first point
+    u = u - u(1);
+    h = h - h(1);
+    s = s - s(1);
+end
+
+function plotEoS(p, ideal, rk, y_label)
+    figure;
+    plot(p / 1e5, ideal, 'b-', 'LineWidth', 2); 
+    hold on;
+    plot(p / 1e5, rk, 'r-', 'LineWidth', 2);
+    xlabel('Pressure [bar]');
+    ylabel(y_label);
+    legend('Ideal EoS', 'R-K EoS');
+    grid on;
+end
+
+%% EoS Comparison based on thermodynamic properties
+%
+% This is th emain subroutine that compares the plots and thermodynamic
+% values obtained using three equations of state.
+
+% Input parameters
+T = 300; % K
+p = 1e5 * linspace(1, 100, 1000); % Pa
+
+% Read the ideal-gas phase
+idealGasPhase = Solution('example_data/co2-thermo.yaml', 'CO2-Ideal');
+[h_ideal, u_ideal, s_ideal, cp_ideal, cv_ideal] = get_thermo_Cantera(idealGasPhase, T, p);
+
+% Read the Redlich-Kwong phase
+redlichKwongPhase = Solution('example_data/co2-thermo.yaml', 'CO2-RK');
+[h_RK, u_RK, s_RK, cp_RK, cv_RK] = get_thermo_Cantera(redlichKwongPhase, T, p);
+
+% Plot the results
+plotEoS(p, u_ideal, u_RK, "Relative Internal Energy [kJ/kg]");
+plotEoS(p, h_ideal, h_RK, "Relative Enthalpy [kJ/kg]"); 
+plotEoS(p, s_ideal, s_RK, "Relative Entropy [kJ/kg-K]");
+
+%%
+% The thermodynamic properties such as internal energy, enthalpy, and 
+% entropy are plotted against the operating pressure at a constant 
+% temperature T = 300 K. The three equations follow each other closely at low 
+% pressures (P < 10bar). However, the ideal gas EoS departs significantly 
+% from the observed behavior of gases near the critical regime 
+% (Pcrit = 73.77 bar).
+%
+% The ideal gas EoS does not consider inter-molecular interactions and the 
+% volume occupied by individual gas particles. At low temperatures and high
+% pressures, inter-molecular forces become particularly significant due to 
+% a reduction in inter-molecular distances. Additionally, at high density, 
+% the volume of individual molecules becomes significant. Both of these 
+% factors contribute to the deviation from ideal behavior at high pressures. 
+% The cubic Redlich-Kwong EoS, on the other hand, predicts thermodynamic 
+% properties accurately near the critical regime.
+
+% Specific heat at constant pressure 
+plotEoS(p, cp_ideal, cp_RK, "C_p [kJ/kg-K]");
+% Specific heat at constant volume
+plotEoS(p, cv_ideal, cv_RK, "C_v [kJ/kg-K]");
+%%
+% In the case of Ideal gas EoS, the specific heats at constant pressure
+% (Cp) and constant volume (Cv) are independent of the pressure. Hence, Cp
+% and Cv for ideal EoS do not change as the pressure is varied from 1 bar
+% to 100 bar in this study.
+%
+% Cp for the R-K EoS follows the trend closely with the Helmholtz EoS from
+% CoolProp up to the critical regime. Alhtough Cp shows reasonable
+% agreement with the Helmholtz EoS in sub-critical and supercritical
+% reginmes, it inaccurately predicts a very high value near the critical
+% point. However, Cp at the critical point is finite for the real fluid.
+% The sudden rise in Cp in the case of the R-K EoS is just a numerical
+% artifact, due to the EoS yielding infinite values in the limiting case,
+% and not a real singularity.
+%
+% Cv, on the other hand, predicts smaller values in the subcritical and
+% critical regime. However, it shows completely incorrect values in the
+% super-critical region, making it invalid at very high pressures. It is
+% well known that the cubic equations typically fail to predics accurate
+% constant-volume heat capacity in the transcritical region [2]_. Certain
+% cubic EoS models have been extended to resolve the discrepancy using
+% crossover models. For further information sese the work of Span [2]_ and
+% Saeed et al. [3]_.
+
+%% Temperature-Density plots
+%
+% The following function plots the T-'rho' diagram over a wide pressure and
+% temperature range. The temperature is varied from 250 K to 400 K. The
+% pressure is changed from 1 bar to 600 bar.
+
+% Input parameters
+% Set up arrays for pressure and temperature
+p_arr = logspace(1, log10(600), 10);
+T_arr = linspace(250, 401, 20);
+p_arr = 1e5 * p_arr(:);
+
+% Initialize matrices to hold densities
+density_ideal = zeros(length(p_arr), length(T_arr));
+density_RK = zeros(length(p_arr), length(T_arr));
+
+% Loop over pressure and temperature to compute densities
+for i = 1:length(p_arr)
+    for j = 1:length(T_arr)
+        T = T_arr(j);
+        P = p_arr(i);
+
+        % Ideal gas
+        idealGasPhase.TP = {T, P};
+        density_ideal(i, j) = idealGasPhase.D;
+
+        % Redlich-Kwong gas
+        redlichKwongPhase.TP = {T, P};
+        density_RK(i, j) = redlichKwongPhase.D;
+
+    end
+end
+
+% Plotting
+figure;
+hold on;
+colors = parula(length(p_arr));
+
+for i = 1:length(p_arr)
+    % Ideal gas lines (dashed)
+    ideal_line = plot(density_ideal(i,:), T_arr, '--', 'Color', colors(i,:), 'HandleVisibility', 'on');
+    % RK EoS lines (circles)
+    RK_line = plot(density_RK(i,:), T_arr, 'o', 'Color', colors(i,:), 'HandleVisibility', 'on');
+end
+xlabel('Density [kg/m^3]');
+
+ylabel('Temperature [K]');
+legend([ideal_line(1), RK_line(1)], {"Ideal EoS", "R-K EoS"}, 'Location', 'northeast'); 
+title('T vs Density for CO2 at Various Pressures');
+grid on;
+
+% Text annotations
+text(30, 320, 'p = 1 bar', 'Color', colors(1,:), 'Rotation', 90);
+text(430, 318, 'p = 97 bar', 'Color', colors(6,:), 'Rotation', -12);
+text(960, 320, 'p = 600 bar', 'Color', colors(10,:), 'Rotation', -68);
+
+toc
+%%
+% The figure compares T-\rho plots for ideal, R-K, and Helmholtz EoS at
+% different operating pressures. All three EoS yield the same plots at low pressures (0
+% bar and 10 bar). However, the Ideal gas EoS departs significantly at high pressures
+% (P > 10 bar), where non-ideal effects are prominent. The R-K EoS closely
+% matches the Helmholtz EoS at supercritical pressures (P > 70 bar). However,
+% it does depart in the liquid-vapor region that exists at P < P_crit
+% and low temperatures (T_crit).
+%
+% .. [1] I.H. Bell, J. Wronski, S. Quoilin, V. Lemort, "Pure and Pseudo-pure Fluid
+%    Thermophysical Property Evaluation and the Open-Source Thermophysical Property
+%    Library CoolProp," Industrial & Engineering Chemistry Research 53 (2014),
+%    https://pubs.acs.org/doi/10.1021/ie4033999
+%
+% .. [2] R. Span, "Multiparameter Equations of State - An Accurate Source of
+%    Thermodynamic Property Data," Springer Berlin Heidelberg (2000),
+%    https://dx.doi.org/10.1007/978-3-662-04092-8
+%
+% .. [3] A. Saeed, S. Ghader, "Calculation of density, vapor pressure and heat capacity
+%    near the critical point by incorporating cubic SRK EoS and crossover translation,"
+%    Fluid Phase Equilibria (2019) 493, https://doi.org/10.1016/j.fluid.2019.03.027
+%
+% sphinx_gallery_thumbnail_number = -1
diff --git a/samples/matlab_experimental/equil.m b/samples/matlab_experimental/equil.m
index e06aaae1dd5..21e87baaa1b 100644
--- a/samples/matlab_experimental/equil.m
+++ b/samples/matlab_experimental/equil.m
@@ -20,7 +20,7 @@ function equil(g)
 
     nsp = gas.nSpecies;
 
-    % find methane, nitrogen, and oxygen indices
+    %% find methane, nitrogen, and oxygen indices
     ich4 = gas.speciesIndex('CH4');
     io2 = gas.speciesIndex('O2');
     in2 = gas.speciesIndex('N2');
@@ -41,7 +41,7 @@ function equil(g)
         xeq(:, i) = gas.X;
     end
 
-    % make plots
+    %% make plots
     clf;
     subplot(1, 2, 1);
     plot(phi, tad);
diff --git a/samples/matlab_experimental/flame.m b/samples/matlab_experimental/flame.m
index 12ba9beeb56..ded43474c3b 100644
--- a/samples/matlab_experimental/flame.m
+++ b/samples/matlab_experimental/flame.m
@@ -32,7 +32,7 @@
     f = Sim1D({left flow right});
 
     % set default initial profiles.
-    rho0 = gas.D;
+    rho0 = gas.massDensity;
 
     % find the adiabatic flame temperature and corresponding
     % equilibrium composition
diff --git a/samples/matlab_experimental/ignite.m b/samples/matlab_experimental/ignite.m
index 239b6c0aa3d..13f7e00625a 100644
--- a/samples/matlab_experimental/ignite.m
+++ b/samples/matlab_experimental/ignite.m
@@ -25,7 +25,7 @@
     gas.TPX = {1001.0, OneAtm, 'H2:2,O2:1,N2:4'};
     gas.basis = 'mass';
     y0 = [gas.U
-          1.0 / gas.D
+          1.0 / gas.massDensity
           gas.Y'];
 
     time_interval = [0 0.001];
diff --git a/samples/matlab_experimental/isentropic.m b/samples/matlab_experimental/isentropic.m
index d0ac15ddf13..fee820e3acb 100644
--- a/samples/matlab_experimental/isentropic.m
+++ b/samples/matlab_experimental/isentropic.m
@@ -1,5 +1,5 @@
 function isentropic(g)
-    %% Isentropic, adiabatic flow
+    %% Adiabatic, isentropic flow
     %
     % In this example, the area ratio vs. Mach number curve is computed for a
     % hydrogen/nitrogen gas mixture.
diff --git a/samples/matlab_experimental/periodic_cstr.m b/samples/matlab_experimental/periodic_cstr.m
index 380bc513006..6d5ddd34cfb 100644
--- a/samples/matlab_experimental/periodic_cstr.m
+++ b/samples/matlab_experimental/periodic_cstr.m
@@ -30,7 +30,7 @@
     tic
     help periodic_cstr
 
-    % create the gas mixture
+    %% create the gas mixture
     gas = Solution('h2o2.yaml', 'ohmech');
 
     % pressure = 60 Torr, T = 770 K
@@ -39,26 +39,27 @@
 
     gas.TPX = {300, p, 'H2:2, O2:1'};
 
-    % create an upstream reservoir that will supply the reactor.  The
-    % temperature, pressure, and composition of the upstream reservoir are
+    %%
+    % Create an upstream reservoir that will supply the reactor.
+    % The temperature, pressure, and composition of the upstream reservoir are
     % set to those of the 'gas' object at the time the reservoir is
     % created.
     upstream = Reservoir(gas);
-
+    %%
     % Now set the gas to the initial temperature of the reactor, and create
     % the reactor object.
     gas.TP = {t, p};
     cstr = IdealGasReactor(gas);
-
+    %%
     % Set its volume to 10 cm^3. In this problem, the reactor volume is
     % fixed, so the initial volume is the volume at all later times.
     cstr.V = 10.0 * 1.0e-6;
-
+    %%
     % We need to have heat loss to see the oscillations. Create a
     % reservoir to represent the environment, and initialize its
     % temperature to the reactor temperature.
     env = Reservoir(gas);
-
+    %%
     % Create a heat-conducting wall between the reactor and the
     % environment. Set its area, and its overall heat transfer
     % coefficient. Larger U causes the reactor to be closer to isothermal.
@@ -67,27 +68,27 @@
     w = Wall(cstr, env);
     w.area = 1.0;
     w.heatTransferCoeff = 0.02;
-
+    %%
     % Connect the upstream reservoir to the reactor with a mass flow
     % controller (constant mdot). Set the mass flow rate to 1.25 sccm.
     sccm = 1.25;
     vdot = sccm * 1.0e-6/60.0 * ((OneAtm / gas.P) * (gas.T / 273.15)); % m^3/s
-    mdot = gas.D * vdot; % kg/s
+    mdot = gas.massDensity * vdot; % kg/s
     mfc = MassFlowController(upstream, cstr);
     mfc.massFlowRate = mdot;
-
+    %%
     % now create a downstream reservoir to exhaust into.
     downstream = Reservoir(gas);
-
+    %%
     % connect the reactor to the downstream reservoir with a valve, and
     % set the coefficient sufficiently large to keep the reactor pressure
     % close to the downstream pressure of 60 Torr.
     v = Valve(cstr, downstream);
     v.valveCoeff = 1.0e-9;
-
+    %%
     % create the network
     network = ReactorNet({cstr});
-
+    %%
     % now integrate in time
     tme = 0.0;
     dt = 0.1;
@@ -104,11 +105,19 @@
         y(3, n) = cstr.massFraction('H2O');
     end
 
+    for i = 2:length(tm)
+        if abs(y(3, i-1) - y(3, i)) > (y(3, i-1) * 5)
+            disp(tm(i))
+        end
+    end
+    %% plot
     clf
     figure(1)
     plot(tm, y)
     legend('H2', 'O2', 'H2O')
     title('Mass Fractions')
+    ylabel('Mass Fractions')
+    xlabel('Time (s)')
 
     toc
 end
diff --git a/samples/matlab_experimental/plug_flow_reactor.m b/samples/matlab_experimental/plug_flow_reactor.m
index d026fcf0315..94e16b03519 100644
--- a/samples/matlab_experimental/plug_flow_reactor.m
+++ b/samples/matlab_experimental/plug_flow_reactor.m
@@ -99,7 +99,7 @@
 
 T_calc(1) = gas_calc.T;
 Y_calc(1, :) = gas_calc.Y;
-rho_calc(1) = gas_calc.D;
+rho_calc(1) = gas_calc.massDensity;
 
 for i = 2:length(x_calc)
 
diff --git a/samples/matlab_experimental/prandtl1.m b/samples/matlab_experimental/prandtl1.m
index a6c038cdc59..2b9b031f374 100644
--- a/samples/matlab_experimental/prandtl1.m
+++ b/samples/matlab_experimental/prandtl1.m
@@ -56,7 +56,7 @@ function prandtl1(g)
 
     disp(['CPU time = ' num2str(cputime - t0)]);
 
-    % plot results
+    %% plot results
 
     clf;
     %figure(1);
diff --git a/samples/matlab_experimental/prandtl2.m b/samples/matlab_experimental/prandtl2.m
index bc93a6992bc..1b56b63d987 100644
--- a/samples/matlab_experimental/prandtl2.m
+++ b/samples/matlab_experimental/prandtl2.m
@@ -55,7 +55,7 @@ function prandtl2(g)
 
     disp(['CPU time = ' num2str(cputime - t0)]);
 
-    % plot results
+    %% plot results
 
     clf;
     subplot(2, 2, 1);
diff --git a/samples/matlab_experimental/rankine.m b/samples/matlab_experimental/rankine.m
index 44de1aedd45..311c597e952 100644
--- a/samples/matlab_experimental/rankine.m
+++ b/samples/matlab_experimental/rankine.m
@@ -10,35 +10,39 @@
 
 help rankine
 
-% Initialize parameters
+%% Initialize parameters
 eta_pump = 0.6;
 eta_turbine = 0.8;
 p_max = 8.0 * OneAtm;
 t1 = 300.0;
 
-% create an object representing water
+%% create an object representing water
 w = Water;
 
+%%
 % start with saturated liquid water at t1
 basis = 'mass';
 w.TQ = {t1, 0};
 h1 = w.H;
 p1 = w.P;
 
+%%
 % pump it to p2
 pump_work = pump(w, p_max, eta_pump);
 h2 = w.H;
 
+%%
 % heat to saturated vapor
 w.PQ = {p_max, 1.0};
 h3 = w.H;
 
 heat_added = h3 - h2;
 
+%%
 % expand adiabatically back to the initial pressure
 turbine_work = expand(w, p1, eta_turbine);
 
-% compute the efficiency
+%% compute the efficiency
 efficiency = (turbine_work - pump_work) / heat_added;
 disp(sprintf('efficiency = %.2f%%', efficiency*100));
 
diff --git a/samples/matlab_experimental/reactor1.m b/samples/matlab_experimental/reactor1.m
index f6d2232a6dc..cd52e6ed232 100644
--- a/samples/matlab_experimental/reactor1.m
+++ b/samples/matlab_experimental/reactor1.m
@@ -20,7 +20,7 @@ function reactor1(g)
     end
 
     P = 101325.0;
-    % set the initial conditions
+    %% set the initial conditions
     gas.TP = {1001.0, P};
     nsp = gas.nSpecies;
     xx = zeros(nsp, 1);
@@ -29,26 +29,29 @@ function reactor1(g)
     xx(48) = 0.573;
     gas.X = xx;
 
-    % create a reactor, and insert the gas
+    %% create a reactor, and insert the gas
     r = IdealGasReactor(gas);
 
-    % create a reservoir to represent the environment
+    %% create a reservoir to represent the environment
     a = Solution('air.yaml', 'air', 'none');
     a.TP = {a.T, P};
     env = Reservoir(a);
 
+    %%
     % Define a wall between the reactor and the environment and
     % make it flexible, so that the pressure in the reactor is held
     % at the environment pressure.
     w = Wall(r, env);
 
+    %%
     % set expansion parameter. dV/dt = KA(P_1 - P_2)
     w.expansionRateCoeff = 1.0e6;
 
+    %%
     % set wall area
     w.area = 1.0;
 
-    % create a reactor network and insert the reactor:
+    %% create a reactor network and insert the reactor:
     network = ReactorNet({r});
 
     nSteps = 100;
diff --git a/samples/matlab_experimental/reactor2.m b/samples/matlab_experimental/reactor2.m
index 06cc112b26c..bc26abd2391 100644
--- a/samples/matlab_experimental/reactor2.m
+++ b/samples/matlab_experimental/reactor2.m
@@ -19,13 +19,13 @@ function reactor2(g)
         gas = Solution('gri30.yaml', 'gri30', 'none');
     end
 
-    % set the initial conditions
+    %% set the initial conditions
     gas.TPX = {1001.0, OneAtm, 'H2:2,O2:1,N2:4'};
 
-    % create a reactor, and insert the gas
+    %% create a reactor, and insert the gas
     r = IdealGasReactor(gas);
 
-    % create a reactor network and insert the reactor
+    %% create a reactor network and insert the reactor
     network = ReactorNet({r});
 
     nSteps = 100;
diff --git a/samples/matlab_experimental/shock_tube.m b/samples/matlab_experimental/shock_tube.m
new file mode 100644
index 00000000000..aa546cbd4fa
--- /dev/null
+++ b/samples/matlab_experimental/shock_tube.m
@@ -0,0 +1,116 @@
+%% Shock-tube species profiles as a function of time
+% Simulate species profiles for a shock tube as a function of time, and
+% observe the impact of incorporating the reduced-pressure linear mixture
+% rule (LMR-R) in such calculations.
+%
+% Here we predict the H2O mole fraction time profiles for a mixture of 1163
+% ppm H2O2/1330 ppm H2O/665 ppm O2/20% CO2/Ar following reflected shock
+% waves (1196 K, 2.127 atm) and compare results against the experimental
+% measurements of Shao et al. [1] Two models are compared in this example:
+%
+% 1. A 2023 model of H2 and NH3 chemistry published by Alzueta et al. [2]
+% 2. An adapted version of this model that has applied the reduced-pressure
+% linear mixture rule (LMR-R) and ab initio third-body efficiencies. [3]
+% References:
+%
+% [1] J. Shao, R. Choudhary, D. F. Davidson, R. K. Hanson, Shock tube/laser
+% absorption measurement of the rate constant of the reaction: H2O2+CO2 =
+% 2OH+CO2, Proc. Combust. Inst. 39 (2023) 735 – 743.
+%
+% [2] M. U. Alzueta, I. Salas, H. Hashemi, P. Glarborg, CO-assisted NH3
+% oxidation, Combust. Flame 257 (2023) 112438.
+%
+% [3] P. J. Singal, J. Lee, L. Lei, R. L. Speth, M. P. Burke,
+% Implementation of New Mixture Rules Has a Substantial Impact on
+% Combustion Predictions for H2 and NH3, Proc. Combust. Inst. 40 (2024)
+% 105779.
+%
+% .. tags:: Matlab, shock tube, kinetics, combustion
+
+clear all;
+close all;
+
+tic
+help shock_tube
+
+file = 'example_data/ammonia-CO-H2-Alzueta-2023.yaml';
+%% Models and colors
+models = struct('Original', 'baseline', 'LMR_R', 'linear-Burke');
+colors = struct('Original', [0.6, 0.6, 0.6], 'LMR_R', [0.5, 0, 0.5]);
+
+results = struct();
+
+%% Experimental data from Shao et al.
+expData.t = [12.3, 20.3, 26.4, 39.6, 58.5, 79.2, 96.1, 113.8, 131.6, ...
+             145.7, 161.2, 181.6, 195.3, 219.9, 237.2, 248.6, 262.4, ...
+             272.2, 280.9];  % microseconds
+expData.X_H2O = [1.47E-03, 1.59E-03, 1.66E-03, 1.78E-03, 1.98E-03, ...
+                 2.06E-03, 2.15E-03, 2.22E-03, 2.26E-03, 2.30E-03, ...
+                 2.39E-03, 2.38E-03, 2.40E-03, 2.42E-03, 2.47E-03, ...
+                 2.53E-03, 2.51E-03, 2.50E-03, 2.47E-03];
+
+%% Loop over mechanisms
+model_names = fieldnames(models);
+for k = 1:length(model_names)
+    model_key = model_names{k};
+    mech_name = models.(model_key);
+
+    % Define mixture composition
+    X_H2O2 = 1163e-6;
+    X_H2O = 1330e-6;
+    X_O2   = 665e-6;
+    X_CO2  = 0.2 * (1 - X_H2O2 - X_H2O - X_O2);
+    X_Ar   = 1 - X_CO2;
+
+    gasComp = 'H2O2: 1163e-6, H2O: 1330e-6, O2: 665e-6, CO2: 0.1994, AR: 0.8006';
+
+    gas = Solution(file, mech_name);
+    gas.TPX = {1196, 2.127 * OneAtm, gasComp};
+
+    r = Reactor(gas,'Reactor');
+    r.energy = 'on';
+    sim = ReactorNet({r});
+
+    time = 0.0;
+    estIgnDelay = 1.0; % seconds
+    counter = 0;
+
+    time_data = [];
+    H2O_X_data = [];
+
+    while time < estIgnDelay
+        time = sim.step();
+        if mod(counter, 10) == 0
+            time_data(end+1) = time * 1e6;  % convert to µs
+            H2O_X_data(end+1) = gas.moleFraction({'H2O'});
+        end
+        counter = counter + 1;
+    end
+
+    % Store results
+    results.(model_key).t = time_data;
+    results.(model_key).X_H2O = H2O_X_data;
+end
+
+%% Plotting
+figure;
+hold on;
+
+for k = 1:length(model_names)
+    model_key = model_names{k};
+    plot(results.(model_key).t, 100 * results.(model_key).X_H2O, ...
+         'LineWidth', 2, 'Color', colors.(model_key), 'DisplayName', strrep(model_key, '_', '-'));
+end
+
+plot(expData.t, 100 * expData.X_H2O, 'ko', 'MarkerSize', 6, ...
+     'MarkerFaceColor', 'white', 'DisplayName', 'Shao et al.');
+
+xlabel('Time [\mus]');
+ylabel('H_2O Mole Fraction [%]');
+legend('Location', 'southwest');
+xlim([0 300]);
+title('Ignition Simulation vs. Experiment');
+grid on;
+
+
+toc
diff --git a/samples/matlab_experimental/test_examples.m b/samples/matlab_experimental/test_examples.m
index 00f389af144..1bcc1b10948 100644
--- a/samples/matlab_experimental/test_examples.m
+++ b/samples/matlab_experimental/test_examples.m
@@ -1,31 +1,64 @@
 % runs selected examples without pausing
-
+clear all
+close all
 ctLoad
+if run_test_examples()
+    disp('test_examples run successfully');
+end
+ctUnload
 clear all
 close all
 
-equil();
-isentropic();
-reactor1();
-reactor2();
-surf_reactor;
-periodic_cstr;
-% plug_flow_reactor;
-lithium_ion_battery
-rankine;
-prandtl1();
-prandtl2();
-% flame1;
-flame2;
-catcomb;
-diff_flame;
-ignite;
-ignite_hp;
-ignite_uv;
-diamond_cvd;
+function success = run_test_examples()
 
-clear all
-close all
-ctUnload
+    success = false;
+
+    examples = {
+        'equil', 'isentropic', 'reactor1', 'reactor2', 'surf_reactor', ...
+        'periodic_cstr', 'plug_flow_reactor', 'lithium_ion_battery', ...
+        'rankine', 'prandtl1', 'prandtl2', 'flame1', 'flame2', ...
+        'catcomb', 'diff_flame', 'ignite', 'ignite_hp', 'ignite_uv', 'diamond_cvd'
+    };
+
+    passed = {};
+    failed = {};
+
+    for idx = 1:length(examples)
+        scriptName = examples{idx};
+        try
+            % Run script in base workshop to protect local variables like
+            % examples
+            evalin('base', sprintf('run(''%s.m'')', scriptName));
+            passed{end+1} = scriptName;
+        catch ME
+            fprintf('An error occurred while running %s: %s\n', scriptName, ME.message);
+            fprintf('Identifier: %s\n', ME.identifier);
+            if strcmp(ME.identifier, 'Cantera:ctError')
+                disp('Caught a CanteraError. Continuing execution...\n');
+            end
+            failed{end+1} = scriptName;
+        end
+    end
+
+        % Summary report
+    disp(' ');
+    disp('============================');
+    disp('Summary of example runs');
+    disp('============================');
+
+    if ~isempty(passed)
+        fprintf('✅ Passed: %s\n', strjoin(passed, ', '));
+    else
+        disp('✅ Passed: (none)');
+    end
+
+    if ~isempty(failed)
+        fprintf('❌ Failed: %s\n', strjoin(failed, ', '));
+    else
+        disp('❌ Failed: (none)');
+        success = true;
+    end
+
+    disp(' ');
 
-disp('Test example run successful.')
+end
diff --git a/samples/matlab_experimental/vapor_dome.m b/samples/matlab_experimental/vapor_dome.m
new file mode 100644
index 00000000000..f8eb21d0f05
--- /dev/null
+++ b/samples/matlab_experimental/vapor_dome.m
@@ -0,0 +1,103 @@
+%% Vapor Dome
+% This example generates a saturated steam table and plots the vapor dome.
+% The steam table corresponds to data typically found in thermodynamic text
+% books and uses the same customary units
+%
+% .. tags:: Matlab, thermodynamics, non-ideal fluid, plotting
+
+clear all
+close all
+
+tic
+help vapor_dome
+
+
+w = Water;
+w.basis = {'mass'};
+
+%% Temperature vector in degC
+degC = [w.minTemp - 273.15, 4, 5, 6, 8, ...
+        10:36, 38, ...
+        40:5:95, 100:10:290, ...
+        300:20:360, 370:1:373, ...
+        w.critTemperature - 273.15];
+
+T_vec = degC + 273.15;
+n = length(T_vec);
+
+%% Preallocate arrays
+T = zeros(n,1); P = zeros(n,1);
+vf = zeros(n,1); vg = zeros(n,1); vfg = zeros(n,1);
+uf = zeros(n,1); ug = zeros(n,1); ufg = zeros(n,1);
+hf = zeros(n,1); hg = zeros(n,1); hfg = zeros(n,1);
+sf = zeros(n,1); sg = zeros(n,1); sfg = zeros(n,1);
+
+%% Loop through temperatures and get saturated liquid/vapor properties
+for i = 1:n
+    T_i = T_vec(i);
+    T(i) = T_i - 273.15;
+
+    % Saturated vapor
+    w.TQ = {T_i, 1};
+    vg(i) = w.V;
+    ug(i) = w.U/1e3;
+    hg(i) = w.H/1e3;
+    sg(i) = w.S/1e3;
+
+    % Saturated liquid
+    w.TQ = {T_i, 0};
+    vf(i) = w.V;
+    uf(i) = w.U/1e3;
+    hf(i) = w.H/1e3;
+    sf(i) = w.S/1e3;
+
+    % Pressure (same for both states)
+    P(i) = w.P/1e5;
+end
+
+%% Delta values
+vfg = vg - vf;
+ufg = ug - uf;
+hfg = hg - hf;
+sfg = sg - sf;
+
+%% Reference state: triple point liquid
+w.TQ = {w.minTemp, 0};
+uf0 = w.U/1e3;
+hf0 = w.H/1e3;
+sf0 = w.S/1e3;
+pv0 = w.P * w.V/1e3;
+
+%% Adjust reference state
+uf = uf - uf0;
+ug = ug - uf0;
+hf = hf - hf0 + pv0;
+hg = hg - hf0 + pv0;
+sf = sf - sf0;
+sg = sg - sf0;
+
+%% Print and write saturated steam table to csv file
+data = table(T, P, vf, vfg, vg, uf, ufg, ug, hf, hfg, hg, sf, sfg, sg);
+disp(data);
+writetable(data, 'saturated_steam_T.csv');
+%% Plot P-v Diagram
+figure;
+semilogx(vf, P, 'b', 'LineWidth', 1.5); hold on;
+semilogx(vg, P, 'r', 'LineWidth', 1.5);
+plot(vg(end), P(end), 'ko', 'MarkerFaceColor', 'w');
+xlabel('Specific Volume v [m^3/kg]');
+ylabel('Pressure P [bar]');
+legend('Saturated Liquid', 'Saturated Vapor', 'Critical Point');
+title('Vapor Dome - P-v Diagram');
+grid on;
+
+%% Plot T-s Diagram
+figure;
+plot(sf, T, 'b', 'LineWidth', 1.5); hold on;
+plot(sg, T, 'r', 'LineWidth', 1.5);
+plot(sg(end), T(end), 'ko', 'MarkerFaceColor', 'w');
+xlabel('Specific Entropy s [kJ/kg-K]');
+ylabel('Temperature T [°C]');
+legend('Saturated Liquid', 'Saturated Vapor', 'Critical Point');
+title('Vapor Dome - T-s Diagram');
+grid on;
diff --git a/test/matlab_experimental/ctTestPath.m b/test/matlab_experimental/ctTestPath.m
deleted file mode 100644
index f7fa212efc0..00000000000
--- a/test/matlab_experimental/ctTestPath.m
+++ /dev/null
@@ -1,23 +0,0 @@
-% ctTestPath.m
-% Set up environment for testing the Cantera Matlab interface
-% from within the Cantera source tree. Run this file from the
-% root of the Cantera source tree, for example:
-%
-%    cd ~/src/cantera
-%    run interfaces/matlab_experimental/testpath.m
-
-% get the list of directories on the Matlab path
-dirs = split(path, pathsep);
-
-% if 'cantera' is already in the path, remove it
-for i = 1:length(dirs)
-    if contains(dirs{i}, 'CANTERA', 'IgnoreCase', true)
-        rmpath(dirs{i});
-    end
-end
-
-cantera_root = getenv('CANTERA_ROOT');
-
-% Add the Cantera toolbox to the Matlab path
-addpath(genpath([cantera_root, '/interfaces/matlab_experimental']));
-addpath(genpath([cantera_root, '/test/matlab_experimental']));
diff --git a/test/matlab_experimental/ctTestSetUp.m b/test/matlab_experimental/ctTestSetUp.m
index 9590f149da7..c58f82d71e1 100644
--- a/test/matlab_experimental/ctTestSetUp.m
+++ b/test/matlab_experimental/ctTestSetUp.m
@@ -1,24 +1,3 @@
-clear all
-
-% Copy library to test folder
-ctTestPath;
-
-if ispc
-    ctName = '/build/lib/cantera_shared.dll';
-elseif ismac
-    ctName = '/build/lib/libcantera_shared.dylib';
-elseif isunix
-    ctName = '/build/lib/libcantera_shared.so';
+function ctTestSetUp()
+    ctLoad();
 end
-% Load Cantera
-if ~libisloaded('libcantera_shared')
-    [~, warnings] = loadlibrary([cantera_root, ctName], ...
-                                [cantera_root, '/include/cantera/clib/ctmatlab.h'], ...
-                                'includepath', [cantera_root, '/include'], ...
-                                'addheader', 'ct', 'addheader', 'ctfunc', ...
-                                'addheader', 'ctmultiphase', 'addheader', ...
-                                'ctonedim', 'addheader', 'ctreactor', ...
-                                'addheader', 'ctrpath', 'addheader', 'ctsurf');
-end
-
-disp('Cantera is loaded for test');
diff --git a/test/matlab_experimental/ctTestTearDown.m b/test/matlab_experimental/ctTestTearDown.m
index bab612e0866..2174435ec89 100644
--- a/test/matlab_experimental/ctTestTearDown.m
+++ b/test/matlab_experimental/ctTestTearDown.m
@@ -1,10 +1,3 @@
-clear all
-% Unload Cantera
-if ispc
-    lib = 'cantera_shared';
-else
-    lib = 'libcantera_shared';
+function ctTestTearDown()
+    ctUnload();
 end
-
-unloadlibrary(lib);
-disp('Cantera has been unloaded');