Как получить больше результатов от ODE45

У меня есть N = 2 связанные нелинейные динамические системы, связь которых задается матрицей 2 на 2 W. Каждая из них описывается n = 8 одой первого порядка. Приведенный ниже код решает эту связанную систему для многих значений параметра p:

for i=1:length(p)
    [t,y(:,:,i)] = ode45(@(t,y) ode(t,y,n,N,W,p(i,:)), t, y0);
end

function [dydt] = ode(t,y,n,N,W,p)
    dydt = NaN(n, N);
    y = reshape(y,[n, N]);    
    y_out = zeros(N,1);
    F_Global = zeros(N,1);
    for i = 1:N
        y_out(i) = y(3,i)-y(4,i);
    end
    for i = 1:N
        F_Global(i) = W(i,:)*sigm(y_out);
    end
    for i = 1:N
        dydt(1,i) = y(5,i);
        dydt(2,i) = y(6,i);
        dydt(3,i) = y(7,i);
        dydt(4,i) = y(8,i);
        dydt(5,i) = sigm(y(3,i) - y(4,i)) - y(5,i) - y(1,i) + F_Global(i);
        dydt(6,i) = sigm(y(3,i) - y(4,i)) - y(6,i) - y(2,i);
        dydt(7,i) = C2*sigm(y(1,i)) + p(i) - y(7,i) - y(3,i);
        dydt(8,i) = sigm(y(2,i)) - y(8,i) - y(4,i);
    end
    dydt = reshape(dydt,n*N, 1);
end

function X = sigm(u)
    ...
end

Внутри функции уже вычислена разница:

y_out(i) = y(3,i)-y(4,i);

Для i = 1,...,N и для всех времен и всех значений p предполагается, что это трехмерная матрица измерений

y_out = size(length(time), length(p), N);

Также внутри функции вычисляются:

F_Global(i) = W(i,:)*sigm(y_out);

который для i = 1,...,N и для всех значений p, но после усреднения по времени, предполагается, что это двумерная матрица измерений

F_Global = size(length(p),N);

Мне нужна помощь, чтобы извлечь y_out и F_Global как выходные данные ode45