Matlab - Método Adams-Bashforth-Moulton para sistemas EDO

Buenas,

Tengo el código de Adams-Bashforth-Moulton (Método multipaso) para solución de UNA ecuación diferencial. Quiero transformarlo para usarlo en un SISTEMA de ecuaciones diferenciales. El problema es ejecutar el código llamando a un archivo F.m donde tengo el sistema.

Archivo F.m;
function Z=F(t,Z)
x=Z(1);
y=Z(2);
Z=[y,-2*y-2*x];
end
(Cada una de las ecuaciones del sistema [ y , -2y-2x ] debe evaluarse en el programa Adams-Bashforth-Moulton, y almacenar el resultado en una matriz donde las filas corresponden al numero de pasos que ejecuta el programa y las columnas al número de ecuaciones del sistema [2 en este caso])

Adams-Bashforth-Moulton:
function A=abm(f,T,Y)

%Input - f is the function entered as a string 'f'
% - T is the vector of abscissas
% - Y is the vector of ordinates
%Remark. The first four coordinates of T and Y must
% have starting values obtained with RK4
%Output - A=[T' Y'] where T is the vector of abscissas and
% Y is the vector of ordinates

n=length(T);

F=zeros(1,4);
F=feval(f,T(1:4),Y(1:4));
h=T(2)-T(1);

for k=4:n-1
if n<5, break, end;
%Predictor
p=Y(k)+(h/24)*(F*[-9 37 -59 55]');
T(k+1)=T(1)+h*k;
F=[F(2) F(3) F(4) feval(f,T(k+1),p)];
%Corrector
Y(k+1)=Y(k)+(h/24)*(F*[1 -5 19 9]');
F(4)=feval(f,T(k+1),Y(k+1));
end

A=[T' Y'];

(En este programa la f se sustituiría por F (sistema) para obtener cada función del sistema evaluada y multiplicar esas evaluaciones por los coeficientes del método (-9 37 -59 55))

¿Alguien tiene alguna idea?
Muchas gracias por la ayuda

Saludos!!