Передача переменных потока в функцию событий в Matlab ODE

Я разрабатываю динамическую модель, описывающую динамику растений и микробов. Рост растений может быть ограничен углеродом (светом), азотом или фосфором, и модель имеет различную динамику в зависимости от того, какой из этих трех элементов является лимитирующим. Растения могут получать азот либо из почвы, либо из атмосферы, и я хочу определить количество питательных веществ, получаемых каждым путем. Итак, на каждом временном шаге я хочу рассчитать, какое питательное вещество ограничивает, а затем рассчитать соответствующие потоки и изменения в пулах. Моя конечная цель — иметь таблицу с данными для всех пулов, потоков и ограничений.

До сих пор для решения этого вопроса я использовал ode45 в MATLAB. Моя функция имеет как события, так и ряд дополнительных значений, которые вычисляются в моем ODE (т. е. я хочу извлечь различные потоки с течением времени, а не только значения dY). Есть ли способ экспортировать все потоки/переменные и использовать функцию событий?

Код, который я использую для передачи дополнительных переменных (xvt, yvt, limvt), имеет форму ниже

(на основе этой публикации MATLAB Central):

function [dydx xvt yvt limvt] = myode(t,input,ps)
persistent xv yv limv

condition=input*ps;

if condition>1: lim=1; x=x*1; y=y/1;
elseif condition==1: lim=2; x=x*2; y=y/2;
elseif condition<1: lim=3; x=x*3; y=y/3;
end

limv = [limv; lim];
yv = [yv; y];
xv = [xv; x];

dydx = x + y;

if nargout>1
      xvt = xv; yvt = yv; limvt = limv;
end

Решение ОДУ можно решить с помощью команд:

[X Y] = ode45(@myode, [0 5], 1);

Затем xv, yv и limv можно получить с помощью:

[dY, xv, yv, limv]=myode([], input,ps);

Чтобы использовать функцию событий, где изменения в limv являются событиями, я установил параметры следующим образом (events_functions выполняет математические операции, чтобы события выдавали правильные значения в нужное время):

options = odeset('Events',@(t0,y,ps)events_function(t0,input,ps,limv,years));

Итак, я могу запустить функцию с событиями (в цикле while, чтобы запускать и останавливать всякий раз, когда происходят события):

[t,xx] = ode23(@myode,[t0:years],input,options,ps);

Однако приведенное выше позволяет мне выполнять либо события на основе limv, либо вычислять limv, но не то и другое одновременно. Я считаю, что этот плакат пытался сделать то же самое. Вы знаете, возможно ли это?