function [t,u]=PCAM2(fun,tI,tF,yI,N)
h=(tF-tI)/N;
t=[tI:h:tF];
u(:,1)=yI;
f(:,1)=feval(fun,t(1),u(:,1));
K1=f(:,1);
K2=feval(fun,t(1)+h,u(:,1)+h*K1);
u(:,2)=u(:,1)+h/2*(K1+K2);
for n=2:N
    f(:,n)=feval(fun,t(n),u(:,n));
    aux=u(:,n)+h/2*(3*f(:,n)-f(:,n-1));
    faux=feval(fun,t(n+1),aux);
    u(:,n+1)=u(:,n)+h/12*(5*faux+8*f(:,n)-f(:,n-1));
end