function [x,y] = pc_orden2(f,xspan,y0,yh,h)

% PC de orden 2, tema 5, test 2, 521230, S2-2011

% intervalo de integracion
xf = xspan(2);
x0 = xspan(1);

% numero de subintervalos en que dividir [t0,tf] segun valor de h
N = ceil((xf-x0)/h);
h = (xf-x0)/N;

% puntos donde se calcularan aproximaciones a los valores de la sol. exacta
% del problema
x = x0:h:xf;
x = x';

% convirtiendo valores iniciales en vectores fila
y0 = y0(:)';
yh = yh(:)';

% inicializando vectores donde se almacenaran las aproximaciones 
y = zeros(N+1,length(y0));

% valores iniciales 
y(1,:) = y0;
y(2,:) = yh;

for i = 3:N+1
    % evaluando en puntos anteriores
    fim1 = feval(f,x(i-1),y(i-1,:));
    fim1 = fim1';
    fim2 = feval(f,x(i-2),y(i-2,:));
    fim2 = fim2';
    % predecir
    y(i,:) = y(i-1,:) + h*((3/2)*fim1 - (1/2)*fim2);
    % evaluar
    fi = feval(f,x(i),y(i,:));
    fi = fi';
    % corregir
    y(i,:) = y(i-1,:) + h*((1/2)*fi + (1/2)*fim1);
end
