Função para Cálculo pelo Método Dernum para Gretl

function d=dernum(Func,u,k)

    // d=dernum(F,x,k)

    // Calcula uma aproximação numérica para a derivada da função F em relação à k-ésima variável, no

    // ponto x, retornando-a em d.

    // F deve ser transmitida como um vetor de "strings", cada uma contendo a função escrita de

    // acordo com a sintaxe do Scilab. P.ex., para o sistema de equações não-lineares

    // {x^2-y^2 = 1,y = e^x}, escreve-se: ["x(1)^2-x(2)^2-1";"exp(x(1))-x(2)"].

    // Baseada em C.T.Kelley, "Iterative Methods for Linear and Nonlinear Equations", Def. 5.4.1, pp.80

    // (C) Rudnei Dias da Cunha 2011

    n = max(size(u))

    s = "v=["

    for i=1:n-1

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

    end

    s = s+Func(n)+"]"

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

    h=sqrt(number_properties('eps'))

    normax=norm(x,'fro')

    if (abs(normax)<sqrt(number_properties("tiny")))

        e=zeros(n,1); e(k)=1

        F1=F(h*e)

        F2=F(x)

        d=(1.0/h)*(F1-F2)

    else

        e=zeros(n,1); e(k)=1

        F1=F(x+h*normax*e)

        F2=F(x)

        d=(1.0/(h*normax))*(F1-F2)

    end

endfunction