MATLAB: сохранение параметров внутри ode45 с использованием 'assignin'

Я запускаю набор ODE с ode45 в MATLAB, и мне нужно сохранить одну из переменных (это не производная) для дальнейшего использования. Я использую функцию assignin, чтобы назначить временную переменную в базовой рабочей области и обновлять ее на каждом этапе. Кажется, это работает, однако размер массива не соответствует размеру вектора решения, полученного из ode45. Например, у меня есть следующая вложенная функция:

function [Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0)

options = odeset('RelTol',1e-7,'AbsTol',1e-7);
[Z,Y] = ode45(@momentum,zspan,Y0,options);

function DY = momentum(z,y)

    DY = zeros(4,1);

    %Entrained Total Velocity
    Ve = sin(theta)*(y(4));

    %Total Relative Velocity
    Urs = sqrt((y(1) - y(4))^2 + (y(2) - Ve*cos(theta))^2 + (y(3))^2);

    %Coefficients
    PSI = K*Urs/y(1);
    PHI = P*Urs/y(1);

    %Liquid Axial Velocity
    DY(1) = PSI*sign(y(1) - y(4))*(1 + (1/6)*(abs(y(1) - y(4))*G)^(2/3));

    %Liquid Radial Velocity
    DY(2) = PSI*sign(y(2) - Ve*cos(theta))*(1 + (1/6)*(abs(y(2) - ...
        Ve*cos(theta))*G)^(2/3));

    %Liquid Tangential Velocity
    DY(3) = PSI*sign(y(3))*(1 + (1/6)*(abs(y(3))*G)^(2/3));

    %Gaseous Axial Velocity
    DY(4) = (1/z/y(4))*((PHI/z)*sign(y(1) - y(4))*(1 + ...
        (1/6)*(abs(y(1) - y(4))*G)^(2/3)) + Ve*Ve - y(4)*y(4));

    assignin('base','Ve_step',Ve);
    evalin('base','Ve_out(end+1) = Ve_step');
end

end

В приведенном выше коде тета (радианы), K (отрицательное значение), P и G являются константами и для данного примера могут быть приняты как любое значение. Zspan - это просто временной шаг интегрирования для решателя ODE, а Y0 - вектор начальных условий (4x1). Опять же, для этого примера они могут принимать любое разумное значение. Теперь в основном файле функция вызывается следующим образом:

Ve_out = 0;
[Z,Y] = droplet_momentum(theta,K,G,P,zspan,Y0);
Ve_out = Ve_out(2:end);

Этот метод работает без жалоб со стороны MATLAB, но проблема в том, что размер Ve_out не совпадает с размером Z или Y. Причина этого в том, что MATLAB вызывает функцию ODE несколько раз для своего алгоритма, поэтому решение будет немного меньше Ve_out. Как предложил am304, я мог бы просто вычислить DY, задав функции ode вектор Z и Y, такой как DY = momentum (Z, Y), однако мне нужно, чтобы это работало с помощью 'assignin' (или аналогичного метода), потому что другой версия этой проблемы имеет неявную зависимость между DY и Ve, и вычисление DY на каждой итерации было бы слишком затратным в вычислительном отношении (я буду запускать эту задачу для многих итераций).