diff --git a/interfaces/matlab/+ct/+impl/call.m b/interfaces/matlab/+ct/+impl/call.m
index 78cacab4ed..8f577a3c57 100644
--- a/interfaces/matlab/+ct/+impl/call.m
+++ b/interfaces/matlab/+ct/+impl/call.m
@@ -6,8 +6,6 @@
     % Convert output to double for better data compatibility
     output = double(output);
 
-    if ismember(output, ct.impl.errorCode)
-        error('Cantera:ctError', ct.impl.getError());
-    end
+    ct.impl.checkErrorCode(output);
 
 end
diff --git a/interfaces/matlab/+ct/+impl/checkErrorCode.m b/interfaces/matlab/+ct/+impl/checkErrorCode.m
new file mode 100644
index 0000000000..b3f9cf4f1e
--- /dev/null
+++ b/interfaces/matlab/+ct/+impl/checkErrorCode.m
@@ -0,0 +1,10 @@
+function checkErrorCode(code)
+    switch code
+        case {-1, -2}
+            error('Cantera:ctError', ct.impl.getError());
+        case -999.999
+            error('Cantera:ctError', ct.impl.getError());
+        case intmax('uint64')
+            error('Cantera:ctError', ct.impl.getError());
+    end
+end
diff --git a/interfaces/matlab/+ct/+impl/errorCode.m b/interfaces/matlab/+ct/+impl/errorCode.m
deleted file mode 100644
index 66854eaef9..0000000000
--- a/interfaces/matlab/+ct/+impl/errorCode.m
+++ /dev/null
@@ -1,3 +0,0 @@
-function err = errorCode()
-    err = [-1, -2, -999.999, double(intmax('uint64'))];
-end
diff --git a/interfaces/matlab/+ct/+impl/getArray.m b/interfaces/matlab/+ct/+impl/getArray.m
index f27065eea2..cdeb165bc1 100644
--- a/interfaces/matlab/+ct/+impl/getArray.m
+++ b/interfaces/matlab/+ct/+impl/getArray.m
@@ -12,14 +12,23 @@
         extraArgs = {};
     end
 
-    buf = clib.array.ctMatlab.Double(buflen);
-    iok = clib.ctMatlab.(funcName)(args{:}, buf, extraArgs{:});
+    persistent arrayCache
+    if isempty(arrayCache)
+        arrayCache = containers.Map('KeyType','double','ValueType','any');
+    end
 
-    iok = double(iok);
-    if ismember(iok, ct.impl.errorCode)
-        error('Cantera:ctError', ct.impl.getError());
+    % Get or create buffer for this size
+    if arrayCache.isKey(buflen)
+        buf = arrayCache(buflen);
+    else
+        buf = clib.array.ctMatlab.Double(buflen);
+        arrayCache(buflen) = buf;
     end
 
+    iok = clib.ctMatlab.(funcName)(args{:}, buf, extraArgs{:});
+
+    ct.impl.checkErrorCode(iok);
+
     % Convert output to double for better data compatibility
     output = buf.double;
 end
diff --git a/interfaces/matlab/+ct/+impl/getString.m b/interfaces/matlab/+ct/+impl/getString.m
index ed88f7fdb4..e0ce38c29d 100644
--- a/interfaces/matlab/+ct/+impl/getString.m
+++ b/interfaces/matlab/+ct/+impl/getString.m
@@ -2,29 +2,40 @@
     % Calls Cantera library functions with string outputs and returns
     % errors if necessary.
 
-    buf = clib.array.ctMatlab.Char(0);
-    buflen = clib.ctMatlab.(funcName)(varargin{:}, buf);
-
-    if buflen > 0
-        buf = clib.array.ctMatlab.Char(buflen);
-
-        nchar = sum(cellfun(@ischar, varargin));
-        if nchar == 0 || nchar == 1
-            iok = clib.ctMatlab.(funcName)(varargin{:}, buf);
-        else
-            error('not implemented - argument list contains more than one char array.')
-        end
-    else
-        error('Cantera:ctError', ct.impl.getError());
+    persistent emptyCache stringCache maxLen
+    if isempty(stringCache)
+        stringCache = clib.array.ctMatlab.Char(256);
+        maxLen = 256;
+    end
+    if isempty(emptyCache)
+        emptyCache = clib.array.ctMatlab.Char(0);
     end
 
-    iok = double(iok);
-    if ismember(iok, ct.impl.errorCode)
+    buflen = clib.ctMatlab.(funcName)(varargin{:}, emptyCache);
+
+    if buflen <= 0
         error('Cantera:ctError', ct.impl.getError());
     end
 
+    % resize stringCache if needed
+    if buflen > maxLen
+        stringCache = clib.array.ctMatlab.Char(buflen);
+        maxLen = buflen;
+    end
+
+    nchar = sum(cellfun(@ischar, varargin));
+    if nchar == 0 || nchar == 1
+        buflen = clib.ctMatlab.(funcName)(varargin{:}, stringCache);
+    else
+        error('not implemented - argument list contains more than one char array.')
+    end
+
+    buflen = double(buflen);
+    ct.impl.checkErrorCode(buflen);
+
     % Discard the last character
-    s = buf.double;
+    s = stringCache.double;
+    s = s(1:buflen);
     if s(end) == 0
         s = s(1:end-1);
     end
diff --git a/interfaces/matlab/Base/+ct/Kinetics.m b/interfaces/matlab/Base/+ct/Kinetics.m
index 1536621ecc..6944d3203c 100644
--- a/interfaces/matlab/Base/+ct/Kinetics.m
+++ b/interfaces/matlab/Base/+ct/Kinetics.m
@@ -22,6 +22,13 @@
         kinID % ID of the Kinetics object.
     end
 
+    properties (Hidden)
+        kinCache struct = struct('nPhases', NaN, ...
+                                 'nTotalSpecies', NaN, ...
+                                 'nReactions', NaN)
+    end
+
+
     properties (SetAccess = protected)
 
         nPhases % Number of phases.
@@ -83,6 +90,31 @@
             obj.kinID = ct.impl.call('mSol_kinetics', id);
         end
 
+        %% Getter methods for cached properties.
+        % These properties rarely change, so we cache their values for efficiency.
+        % But make sure to update the cache if a certain method modifies them.
+
+        function n = get.nPhases(obj)
+            if isnan(obj.kinCache.nPhases)
+                obj.kinCache.nPhases = ct.impl.call('mKin_nPhases', obj.kinID);
+            end
+            n = obj.kinCache.nPhases;
+        end
+
+        function n = get.nTotalSpecies(obj)
+            if isnan(obj.kinCache.nTotalSpecies)
+                obj.kinCache.nTotalSpecies = ct.impl.call('mKin_nTotalSpecies', obj.kinID);
+            end
+            n = obj.kinCache.nTotalSpecies;
+        end
+
+        function n = get.nReactions(obj)
+            if isnan(obj.kinCache.nReactions)
+                obj.kinCache.nReactions = ct.impl.call('mKin_nReactions', obj.kinID);
+            end
+            n = obj.kinCache.nReactions;
+        end
+
         %% Get scalar attributes
 
         function n = kineticsSpeciesIndex(obj, name)
@@ -111,18 +143,6 @@
             n = ct.impl.call('mKin_multiplier', obj.kinID, irxn - 1);
         end
 
-        function n = get.nPhases(obj)
-            n = ct.impl.call('mKin_nPhases', obj.kinID);
-        end
-
-        function n = get.nReactions(obj)
-            n = ct.impl.call('mKin_nReactions', obj.kinID);
-        end
-
-        function n = get.nTotalSpecies(obj)
-            n = ct.impl.call('mKin_nTotalSpecies', obj.kinID);
-        end
-
         function n = phaseIndex(obj, phase)
             % The index of a specific phase. ::
             %
diff --git a/interfaces/matlab/Base/+ct/Mixture.m b/interfaces/matlab/Base/+ct/Mixture.m
index ea3960e727..1efe0399d6 100644
--- a/interfaces/matlab/Base/+ct/Mixture.m
+++ b/interfaces/matlab/Base/+ct/Mixture.m
@@ -40,6 +40,12 @@
         P % Pressure [Pa].
     end
 
+    properties (Hidden)
+        mixCache struct = struct('nElements', NaN, ...
+                                 'nSpecies', NaN, ...
+                                 'nPhases', NaN)
+    end
+
     properties (SetAccess = protected)
 
         phases % Phases in the mixture
@@ -150,28 +156,46 @@ function addPhase(obj, phase, moles)
 
             ct.impl.call('mMix_addPhase', obj.mixID, phase.tpID, moles);
 
-        end
-
-        %% Mixture Get methods
+            % Reset cached properties
+            obj.mixCache.nElements = NaN;
+            obj.mixCache.nPhases = NaN;
+            obj.mixCache.nSpecies = NaN;
 
-        function temperature = get.T(obj)
-            temperature = ct.impl.call('mMix_temperature', obj.mixID);
         end
 
-        function pressure = get.P(obj)
-            pressure = ct.impl.call('mMix_pressure', obj.mixID);
-        end
+        %% Getter methods for cached properties.
+        % These properties rarely change, so we cache their values for efficiency.
+        % But make sure to update the cache if a certain method modifies them.
 
         function n = get.nElements(obj)
-            n = ct.impl.call('mMix_nElements', obj.mixID);
+            if isnan(obj.mixCache.nElements)
+                obj.mixCache.nElements = ct.impl.call('mMix_nElements', obj.mixID);
+            end
+            n = obj.mixCache.nElements;
         end
 
         function n = get.nPhases(obj)
-            n = ct.impl.call('mMix_nPhases', obj.mixID);
+            if isnan(obj.mixCache.nPhases)
+                obj.mixCache.nPhases = ct.impl.call('mMix_nPhases', obj.mixID);
+            end
+            n = obj.mixCache.nPhases;
         end
 
         function n = get.nSpecies(obj)
-            n = ct.impl.call('mMix_nSpecies', obj.mixID);
+            if isnan(obj.mixCache.nSpecies)
+                obj.mixCache.nSpecies = ct.impl.call('mMix_nSpecies', obj.mixID);
+            end
+            n = obj.mixCache.nSpecies;
+        end
+
+        %% Mixture Get methods
+
+        function temperature = get.T(obj)
+            temperature = ct.impl.call('mMix_temperature', obj.mixID);
+        end
+
+        function pressure = get.P(obj)
+            pressure = ct.impl.call('mMix_pressure', obj.mixID);
         end
 
         function mu = get.chemPotentials(obj)
diff --git a/interfaces/matlab/Base/+ct/ThermoPhase.m b/interfaces/matlab/Base/+ct/ThermoPhase.m
index 9710b27305..3eddea9d22 100644
--- a/interfaces/matlab/Base/+ct/ThermoPhase.m
+++ b/interfaces/matlab/Base/+ct/ThermoPhase.m
@@ -10,6 +10,17 @@
     % :param id:
     %     Integer ID of the solution holding the :mat:class:`ct.ThermoPhase` object.
 
+    properties (Hidden)
+        thermoCache struct = struct('nElements', NaN, ...
+                                    'nSpecies', NaN, ...
+                                    'minTemp', NaN, ...
+                                    'maxTemp', NaN, ...
+                                    'refPressure', NaN, ...
+                                    'molecularWeights', [], ...
+                                    'atomicWeights', [], ...
+                                    'charges', [])
+    end
+
     properties (SetAccess = public)
 
         X % Mole fractions.
@@ -745,18 +756,74 @@ function display(obj)
             str = ct.impl.getString('mThermo_report', obj.tpID, 1, threshold);
         end
 
-        %% Single-property getter methods
+        %% Getter methods for cached properties.
+        % These properties rarely change, so we cache their values for efficiency.
+        % But make sure to update the cache if a certain method modifies them.
+
+        function n = get.nSpecies(obj)
+            if isnan(obj.thermoCache.nSpecies)
+                obj.thermoCache.nSpecies = ct.impl.call('mThermo_nSpecies', obj.tpID);
+            end
+            n = obj.thermoCache.nSpecies;
+        end
+
+        function n = get.nElements(obj)
+            if isnan(obj.thermoCache.nElements)
+                obj.thermoCache.nElements = ct.impl.call('mThermo_nElements', obj.tpID);
+            end
+            n = obj.thermoCache.nElements;
+        end
+
+        function t = get.maxTemp(obj)
+            if isnan(obj.thermoCache.maxTemp)
+                obj.thermoCache.maxTemp = ct.impl.call('mThermo_maxTemp', obj.tpID, -1);
+            end
+            t = obj.thermoCache.maxTemp;
+        end
+
+        function t = get.minTemp(obj)
+            if isnan(obj.thermoCache.minTemp)
+                obj.thermoCache.minTemp = ct.impl.call('mThermo_minTemp', obj.tpID, -1);
+            end
+            t = obj.thermoCache.minTemp;
+        end
+
+        function p = get.refPressure(obj)
+            if isnan(obj.thermoCache.refPressure)
+                obj.thermoCache.refPressure = ct.impl.call('mThermo_refPressure', obj.tpID);
+            end
+            p = obj.thermoCache.refPressure;
+        end
+
+        function mw = get.molecularWeights(obj)
+            if isempty(obj.thermoCache.molecularWeights)
+                nsp = obj.nSpecies;
+                obj.thermoCache.molecularWeights = ...
+                    ct.impl.getArray('mThermo_getMolecularWeights', nsp, obj.tpID);
+            end
+            mw = obj.thermoCache.molecularWeights;
+        end
 
         function amu = get.atomicWeights(obj)
-            nel = obj.nElements;
-            amu = ct.impl.getArray('mThermo_atomicWeights', nel, obj.tpID);
+            if isempty(obj.thermoCache.atomicWeights)
+                nel = obj.nElements;
+                obj.thermoCache.atomicWeights = ...
+                    ct.impl.getArray('mThermo_atomicWeights', nel, obj.tpID);
+            end
+            amu = obj.thermoCache.atomicWeights;
         end
 
         function e = get.charges(obj)
-            nsp = obj.nSpecies;
-            e = ct.impl.getArray('mThermo_getCharges', nsp, obj.tpID);
+            if isempty(obj.thermoCache.charges)
+                nsp = obj.nSpecies;
+                obj.thermoCache.charges = ...
+                    ct.impl.getArray('mThermo_getCharges', nsp, obj.tpID);
+            end
+            e = obj.thermoCache.charges;
         end
 
+        %% Single-property getter methods
+
         function c = get.cv(obj)
             if strcmp(obj.basis, 'molar')
                 c = ct.impl.call('mThermo_cv_mole', obj.tpID);
@@ -801,14 +868,6 @@ function display(obj)
             b = ct.impl.call('mThermo_isothermalCompressibility', obj.tpID);
         end
 
-        function t = get.maxTemp(obj)
-            t = ct.impl.call('mThermo_maxTemp', obj.tpID, -1);
-        end
-
-        function t = get.minTemp(obj)
-            t = ct.impl.call('mThermo_minTemp', obj.tpID, -1);
-        end
-
         function mmw = get.meanMolecularWeight(obj)
             mmw = ct.impl.call('mThermo_meanMolecularWeight', obj.tpID);
         end
@@ -826,23 +885,6 @@ function display(obj)
             c = ct.impl.getArray('mThermo_getConcentrations', nsp, obj.tpID);
         end
 
-        function mw = get.molecularWeights(obj)
-            nsp = obj.nSpecies;
-            mw = ct.impl.getArray('mThermo_getMolecularWeights', nsp, obj.tpID);
-        end
-
-        function nel = get.nElements(obj)
-            nel = ct.impl.call('mThermo_nElements', obj.tpID);
-        end
-
-        function nsp = get.nSpecies(obj)
-            nsp = ct.impl.call('mThermo_nSpecies', obj.tpID);
-        end
-
-        function p = get.refPressure(obj)
-            p = ct.impl.call('mThermo_refPressure', obj.tpID);
-        end
-
         function p = get.satPressure(obj)
             p = ct.impl.call('mThermo_satPressure', obj.tpID, obj.T);
         end
diff --git a/interfaces/matlab/Utility/+ct/cleanUp.m b/interfaces/matlab/Utility/+ct/cleanUp.m
index 917b96e6f5..1d9b0d9fdc 100644
--- a/interfaces/matlab/Utility/+ct/cleanUp.m
+++ b/interfaces/matlab/Utility/+ct/cleanUp.m
@@ -3,10 +3,13 @@ function cleanUp()
 
     ct.isLoaded(true);
 
-classList = {'ct.Interface', 'ct.Kinetics', 'ct.Mixture', 'ct.ThermoPhase', ...
-             'ct.Transport', 'ct.Solution', 'ct.Func1', 'ct.oneD.Domain', ...
-             'ct.oneD.Sim1D', 'ct.zeroD.Connector', 'ct.zeroD.ReactorBase', ...
-             'ct.zeroD.ReactorNet'};
+    % Clear the persistent cache
+    clear getArray getString
+
+    classList = {'ct.Interface', 'ct.Kinetics', 'ct.Mixture', 'ct.ThermoPhase', ...
+                 'ct.Transport', 'ct.Solution', 'ct.Func1', 'ct.oneD.Domain', ...
+                 'ct.oneD.Sim1D', 'ct.zeroD.Connector', 'ct.zeroD.ReactorBase', ...
+                 'ct.zeroD.ReactorNet'};
 
     varList = evalin('base', 'whos');
 
diff --git a/interfaces/matlab/Utility/+ct/gitCommit.m b/interfaces/matlab/Utility/+ct/gitCommit.m
index e97676f893..a0fff7dc4a 100644
--- a/interfaces/matlab/Utility/+ct/gitCommit.m
+++ b/interfaces/matlab/Utility/+ct/gitCommit.m
@@ -7,5 +7,5 @@
     %     A string containing the Git commit hash for the current version of Cantera.
 
     ct.isLoaded(true);
-    v = ct.impl.getString('mCt_getGitCommit');
+    v = ct.impl.getString('mCt_gitCommit');
 end
diff --git a/interfaces/matlab/Utility/+ct/unload.m b/interfaces/matlab/Utility/+ct/unload.m
index ea12fa199f..39e2e281d3 100644
--- a/interfaces/matlab/Utility/+ct/unload.m
+++ b/interfaces/matlab/Utility/+ct/unload.m
@@ -11,6 +11,9 @@ function unload()
                 "cleanUp failed (%s).", ME.message);
     end
 
+    % Clear the persistent cache even if cleanUp failed
+    clear getString getArray
+
     if ct.executionMode() == "inprocess"
         warning("ct.unload:UnloadFailed", ...
                 ("Unloading of `ctMatlab` library is not supported for " + ...
diff --git a/samples/matlab/plug_flow_reactor.m b/samples/matlab/plug_flow_reactor.m
index c463ab37fb..6533a0b0bf 100644
--- a/samples/matlab/plug_flow_reactor.m
+++ b/samples/matlab/plug_flow_reactor.m
@@ -94,6 +94,7 @@
 
 x_calc = 0:dx:L;
 nsp = gas_calc.nSpecies;
+Ru = ct.GasConstant;
 
 % Initialize arrays for T, Y, and rho at each location:
 T_calc = zeros(length(x_calc), 1);
@@ -125,7 +126,7 @@
     %--------------------------------------------------------------------------
     % These values are passed onto the ode15s solver
     [~, y] = ode15s(@PFR_solver, limits, inlet_soln, options, ...
-                    gas_calc, mdot_calc, A_in, dAdx, k, nsp);
+                    gas_calc, mdot_calc, A_in, dAdx, k, nsp, Ru);
 
     T_calc(i) = y(end, 2);
     rho_calc(i) = y(end, 1);
@@ -191,7 +192,7 @@
 % The integrator integrates the derivatives spatially, to solve the density,
 % temperature, and species mass fraction profiles as a function of distance x.
 
-function F = PFR_solver(x, soln_vector, gas, mdot, A_in, dAdx, k, nsp)
+function F = PFR_solver(x, soln_vector, gas, mdot, A_in, dAdx, k, nsp, Ru)
     % Plug flow reactor governing equations
 
     rho = soln_vector(1);
@@ -211,7 +212,6 @@
     gas.TDY = {T, rho, Y};
 
     MW_mix = gas.meanMolecularWeight;
-    Ru = ct.GasConstant;
     R = Ru / MW_mix;
     vx = mdot / (rho * A);
     P = rho * R * T;