function [x,y] = rkpreg1_tema8(f,xspan,y0,h)
% rk en tema 8, test 2, 521230, S2-2011
xf = xspan(2);
x0 = xspan(1);

% numero de subintervalos
N = ceil((xf-x0)/h);
h = (xf-x0)/N;

% puntos donde se buscan aproximaciones a valores de solucion exacta del
% PVI con parte derecha f
x = x0:h:xf;
x = x';
% aproximaciones a solucion de PVI
% tantas columnas como componentes tenga y0
y = zeros(N+1,length(y0));
y0 = y0(:)';
y(1,:) = y0;

for i = 1: N
    k1 = h*feval(f,x(i),y(i,:));
    k2 = h*feval(f,x(i)+h/2,y(i,:)+k1/2);
    k3 = h*feval(f,x(i)+h/4,y(i,:)+(3*k1 + k2)/16);
    k4 = h*feval(f,x(i)+h/2,y(i,:)+ k3/2);
    k5 = h*feval(f,x(i)+3*h/4,y(i,:)+ (-3*k2 + 6*k3 + 9*k4)/16);
    k6 = h*feval(f,x(i)+h,y(i,:)+ (k1 + 4*k2 + 6*k3 - 12*k4 + 8*k5)/7);
    y(i+1,:) = y(i,:) + (7*k1 + 32*k3 + 12*k4 + 32*k5 + 7*k6)/90;
end