Using a value calculated by 'OutputFcn' as a parameter in the set-up of the system of ODEs solver

Question

Is it possible, without global variables, to have the function that describes the system of ODEs that uses an output calclated by from an OutputFcn in order to decide how the ODE is structured?

Problem

There is a system of ODEs that describes the heat transfer through a wall. At a given node in the wall, there is a heat loss that is dependent on the node's temperature. I would like to inregrate this total heat loss. Once the heat loss reaches a certain threshold this heat loss should stop and the heat loss term disappears from the local ODE. I have a working solution, but I want to be able to improve it so that I can run many cases using parfor loops. At the moment, this isn't possible because I have used global variables to track the heat flow that I am integrating.

Current Solution

I have done this by declaring a series of global variables. I use an OutputFcn to calculate the heat flow at each time step. In the function that describes the ODEs, there is an IF statement that checks whether or not the current time is greater than the "timeEnd" value. If it is greater, then the ODE does not include the heat loss term.

function [status] = myOutputFcn(t,y,flag)
global t_i T_1 T_2 Q_i threshold timeEnd index_1 index_2

switch flag
    case 'init'
    case ''
        %%% Save the temperature and time matricies
        T_1 = [T_1;y(index_1)'];
        T_2 = [T_2;y(index_2)'];
        t_i = [t_i;t];
        %%% Integrate the difference between T_2 and T_1
        Q = trapz(t,T_2-T_1);
        %%% check if the integrated Q is above the threshold for the first time
        if (max(Q)>threshold)&&(timeEnd==0)
               timeEnd=t;
        end
    case 'done'
end
status = 0;
end

As you can see, the size of the time and temperature matricies are changing with each time step. This becomes quite cumbersome as I have upwards of 40k timesteps and hundreds of equations that are solved simultaneously.