function [t,u]=RK4(f,tI,tF,yI,N)
h=(tF-tI)/N;
t=[tI:h:tF];
u(:,1)=yI;
for n=1:N
    K1=feval(f,t(n),u(:,n));
    K2=feval(f,t(n)+h/2,u(:,n)+h/2*K1);
    K3=feval(f,t(n)+h/2,u(:,n)+h/2*K2);
    K4=feval(f,t(n)+h,u(:,n)+h*K3); 
    u(:,n+1)=u(:,n)+h/6*(K1+2*K2+2*K3+K4);
end