Resolve o problema de valor inicial na forma autônoma para um sistema de EDOs usando o método previsor-corretor Adams-Bashforth

function Y=ab4_am4(F,x,h,t1)

    // Y=ab4_am4(F,x,h,t1)

    // Resolve o problema de valor inicial na forma autônoma para um

    // sistema de EDOs usando o método previsor-corretor Adams-Bashforth/

    // Adams-Moulton de 4a ordem.

    // Os valores iniciais encontram-se no vetor "x" passado como

    // argumento de entrada, e o primeiro elemento

    // de "x" contém o valor inicial de t, t_0.

    // A função F pode ser passada de duas formas:

    // a) como um vetor coluna contendo as funções expressas como "strings",

    //    com as variáveis de interesse sendo especificadas em termos dos

    //    elementos do vetor "x", como por exemplo:

    //    ['1' ; 'x(1)' ; 'x(1)*x(3)' ; 'x(1)*x(2)-x(3)'];

    // b) como uma função definida externamente no Scilab, utilizando-se

    //    um arquivo com extensão ".sci", a qual deve ter um vetor coluna

    //    como argumento de entrada, o qual contém os valores das variáveis

    //    de interesse, e devolver um outro vetor coluna, o qual conterá

    //    os resultados das funções avaliadas nos valores das variáveis de

    //    interesse. Ambos os vetores devem ter a mesma quantidade de elementos.

    //    As variáveis de interesse devem ser especificadas em termos dos

    //    elementos do vetor de entrada, como no exemplo abaixo (equivalente

    //    ao exemplo no item "a"):

    //       function Y=F(X)

    //           Y(1) = 1

    //           Y(2) = X(1)

    //           Y(3) = X(1)*X(3)

    //           Y(4) = X(1)*X(2)-X(3)

    //       endfunction

    // (C) Rudnei Dias da Cunha 2015

    // Definição da função F, se passada como vetor coluna de funções

    // expressas como strings

    n = max(size(x));

    if type(F)==10 then

        m = max(size(F));

        if n~=m then

            mprintf('Erro, dimensões incompatíveis entre F e x\n')

            Y = [];

            return

        end

        s = "v=[";

        for i=1:m-1

            s = s+F(i)+";";

        end

        s = s+F(m)+"]";

        deff('[v]=F(x)',s);

    end

    // Iteração

    F0 = h*F(x);

    x=x+F0;

    F1 = h*F(x);

    x=x+F1;

    F2 = h*F(x);

    x=x+F2;

    F3 = h*F(x);

    k = 1

    while x(1)<t1

        mprintf("k=%d ",k) ; disp(x')

        Y(k,:) = x'

        y = x+1.0/24.0*(55*F3-59*F2+37*F1-9*F0);

        F0 = F1;

        F1 = F2;

        F2 = F3;

        F3 = h*F(y);

        x = x+1.0/24.0*(9*F3+19*F2-5*F1+F0);

        k = k+1

    end

endfunction