function [y]=RK4(y0,x0,xf,f,Nmax)

// Esta função integra a EDO y'=f(x,y) com condição incial y(x0)=y0 
// através do método Runge-Kutta Clássico (4 etapas) no itervalo (x0,xf)
// em Nmax passos (pontos).
//
// Sistemas de EDO também podem ser integrados através desta função,
// basta definir a função f como uma função R^(n+1)->R^n
// A saída fica armazenada nas colunas da matriz y (a i-ésima coluna
// corresponde aos valores das componentes da aproximação em x0+i*h ).
// 
// Variáveis auxiliares
// h  : espaçamento entre os nós da integração.
// x_0: (i-1)-ésimo nó da integração.
// y_0: valor da aproximação em x_0.
// yi : valor da aproximação em no i-ésimo nó da integração.
// i  : contador.
// k1 : primeira etapa do método Runge-Kutta Clássico. 
// k2 : segunda etapa do método Runge-Kutta Clássico.
// k3 : terceira etapa do método Runge-Kutta Clássico.
// k4 : quarta etapa do método Runge-Kutta Clássico.

// incialização. 

h=(xf-x0)/Nmax

x_0=x0
y_0=y0
y=zeros(size(y0,1),Nmax) // inicialização da saída.

// laço principal

for i=1:Nmax

// 4 etapas
    k1=f(x_0,y_0)
    k2=f(x_0+0.5*h,y_0+0.5*h*k1)
    k3=f(x_0+0.5*h,y_0+0.5*h*k2)
    k4=f(x_0+h,y_0+h*k3)
    yi=y_0+h/6*(k1+2*(k2+k3)+k4)

    y(:,i)=yi
 
    x_0=x_0+h 
    y_0=yi
 
end

endfunction
