function [t,u]=PC3N(fun,t0,y0,T,N)
h=T/N;
t=[t0:h:t0+T];
u(1)=y0;
for i=1:2
    f(i)=fun(t(i),u(i));
    K1=f(i);
    K2=fun(t(i)+h/2,u(i)+h/2*K1);
    K3=fun(t(i)+h/2,u(i)+h/2*K2);
    K4=fun(t(i)+h,u(i)+h*K3);
    u(i+1)=u(i)+h/6*(K1+2*K2+2*K3+K4);
end
for i=3:N
    f(i)=fun(t(i),u(i));
    u(i+1)=u(i)+h/2*(3*f(i)-f(i-1));
    for k=1:2
    f(i+1)=fun(t(i+1),u(i+1));
    u(i+1)=u(i)+h/24*(9*f(i+1)+19*f(i)-5*f(i-1)+f(i-2));
    end
end