
function [time,w] = rk4(f,a,b,ic,N)



% a,b -- endpoints
% ic -- initial conditions
% f -- function we're integrating
% N -- partitions

w=ic;
h=(b-a)/N;
t=a;
time=t;




for i=2:N+1

   k1(i,:)=h*feval(f,t,w(i-1,:));
   k2(i,:)=h*feval(f,t+h/2,w(i-1,:)+k1(i,:)/2);
   k3(i,:)=h*feval(f,t+h/2,w(i-1,:)+k2(i,:)/2);
   k4(i,:)=h*feval(f,t+h,w(i-1,:)+k3(i,:));

   new=w(i-1,:) + (1/6)*(k1(i,:) +2*k2(i,:) + 2*k3(i,:) +k4(i,:));
   t=a+(i-1)*h;
   time=[time;t];
   w=[w;new];

end %for loop
