如何从ODE45接收更多输出

Question

我有N = 2个耦合的非线性动力学系统，其耦合由一个2 x 2矩阵W给出。它们中的每一个都由n = 8个一阶ode描述。下面的代码针对参数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的所有值，这应该是尺寸为3维的矩阵

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的输出

Answer 1

这应该工作

for i=1:length(p)
    [t,y(:,:,i)] = ode45(@(t,y) ode(t,y,n,N,W,p(i,:)), t, y0);
end
yout = y(:,3:8:end,:) - y(:,4:8:end,:);
WW = repmat(W(:)',[size(y,1) 1 size(y,3)]);
F_global = WW .* yout;