function [y,x,hc,rp]=dopri5(y0,x0,xf,%fun,rtol,atol)
    // A função dopri5 implementa o método Dormand-Prince5(4)
    // com passo variável e adptado à tolerância do erro estimado
    // em norma sup a partir das aproximações com erro global 
    // O(h^5) e O(h^4) para o P.V.I.:
    //
    // y'=f(x,y), x>x0,
    // y(x0)=y0.
    //
    // O passo inicial é calculado a partir do algoritmo proposto
    // em Hairer,E;Nørsett,S.P.;Wanner,G. "Solving Ordinary
    // Differential Equations I. Nonstiff Problems" 2nd rev.ed.(1993)
    // pg. 169.
    // As estimativas para o cálculo dos passos seguem as ideias 
    // expostas nessa mesma referência às páginas 167 e 168.
    
    // Variáveis de entrada
    // y0  : condição inicial.
    // x0  : valor inicial da variável independente.
    // xf  : valor final da variável independente.
    // f   : lado direito da EDO.
    // rtol: tolerância na estimativa do erro relativo.
    // atol: tolerância na estimativa do erro absoluto.
    // As os parâmetros rtol e atol podem ser matrizes linha.
    // Neste caso, as tolerâncias aplicam-se a cada componente.
    
    // Variáveis de saída
    // y : aproximações para a solução em x.
    // x : valor da variável independente em cada nó.
    // hc: valor do espaçamento utilizado.
    // rp: iteradas nas quais houve rejeição de passo.
    
    // Variáveis auxiliares
    // a        : matriz 7x7 com os parâmetros do método. 
    // b        : matriz 1x7 com os parâmetros do método (b'). 
    // c        : matriz 7x1 com os parâmetros do método. 
    // fac      : fator de "segurança" no cálculo dos passos.
    // facmax   : fator máximo no cálculo dos passos.
    // facmin   : fator mínimo no cálculo dos passos.
    // dim      : dimensão da EDO.
    // xi       : varável indepente.
    // yi       : aproximação para a solução em xi.
    // prossegue: varável booleana para o controle do laço.
    // i        : contador do número de aproximações aceitas.
    // r        : contador do número de passos rejeitados.
    // d0       : valor da norma (sup) de y0.
    // d1       : valor da norma (sup) de y'(x0).
    // d2       : valor da norma (sup) de uma aprox. para y''(x0).
    // h0       : estimativa para o passo incial.
    // h1       : estimativa para o passo incial.
    // h        : valor do passo.
    // hnew     : estimativa para o passo.
    // k1 a k6  : estágios do método dopri5(4).
    // k7       : corresponde à matriz b' da estimativa O(h^4)
    // y4       : estimativa O(h^4).
    // y5       : estimativa O(h^5).


    // Inicialização
    xi=x0        
    yi=y0        
    dim=size(y0,1)
    fac=0.8    
    facmax=5 
    facmin=0.2
    
    a=zeros(7,7)
    a(2,1)=1/5
    a(3,1:2)=[3/40,9/40]
    a(4,1:3)=[44/45,-56/15,32/9]
    a(5,1:4)=[19372/6561,-25360/2187,64448/6561,-212/729]
    a(6,1:5)=[9017/3168,-355/33,46732/5247,49/176,-5103/18656]
    a(7,1:6)=[35/384,0,500/1113,125/192,-2187/6784,11/84]
    b=[5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40]
    c=[0 1/5 3/10 4/5 8/9 1 1]'

    i=0
    r=1
    rp=[]
    prossegue=%T 

    // estimativa para o passo inicial.
    d0=norm(y0,'inf')
    d1=norm(%fun(x0,y0),'inf')
    if min(d0,d1)<1d-5 then
        h0=1d-6
    else
        h0=0.01*d0/d1
    end
    d2=norm(%fun(x0+h0,y0+h0*%fun(x0,y0))-%fun(x0,y0),'inf')/h0
    if max(d1,d2)<=1d-15 then
        h1=max(1d-6,h0/1000)
    else
        h1=(0.01/max(d1,d2))^(1/5)
    end
    h=min(100*h0,h1)
    
    
    // laço principal
    while prossegue

        k1=%fun(xi,yi)
        k2=%fun(xi+h*c(2),yi+h*a(2,1)*k1)
        k3=%fun(xi+h*c(3),yi+h*(a(3,1)*k1+a(3,2)*k2))
        k4=%fun(xi+h*c(4),yi+h*(a(4,1)*k1+a(4,2)*k2+a(4,3)*k3))
        k5=%fun(xi+h*c(5),yi+h*(a(5,1)*k1+a(5,2)*k2+a(5,3)*k3+a(5,4)*k4))
        k6=%fun(xi+h*c(6),yi+h*(a(6,1)*k1+a(6,2)*k2+a(6,3)*k3+a(6,4)*k4+a(6,5)*k5))
        
        y5=yi + h*(a(7,1)*k1+a(7,3)*k3+a(7,4)*k4+a(7,5)*k5+a(7,6)*k6)
        
        k7=%fun(xi+h*c(7),y5)
        
        y4=yi + h*(b(1)*k1+b(3)*k3+b(4)*k4+b(5)*k5+b(6)*k6+b(7)*k7)
        
        // Estimativa para o novo passo, passo "ótimo": h*(1/err)^(1/5) 
        sc=atol + max(norm(y5,'inf'),norm(yi,'inf')).*rtol
        err=norm((y5-y4)./sc,'inf')
        hnew=h*min(facmax,max(facmin,fac*(1/err)^(1/5))) // 
        
        // O passo é aceito se err<=1.
        if err<=1 then
            // O passo é aceito se não ultrapassar
            // o final da integração.
            if (xi+h)<=xf then
                i=i+1 
                xi=xi+h
                yi=y5
                x(i)=xi
                y(:,i)=yi
                hc(i)=h
                h=hnew
            else
            // Se o passo for maior do que o necesssário para
            // encerra ra integração, a aproximação é  recalculada 
            // com o passo correto.
                h=xf-xi
            end             
        else
            // Se err>1 o passo é rejeitado e nova aproximação 
            // é calculada com h=hnew.
            h=hnew
            r=r+1
            rp(r)=i
        end
        prossegue = xi<xf 
    end
endfunction
