function [x,t]=rk4(fun,x0,t0,tk,n)
%[x,t]=rk4(fun,x0,t0,tk,n)
%finds an approximate solution the differential equation
%x'(t)=fun(t,x(t))
%x(t0)=x0
%on the interval [t0,tk] using the Runge-Kutta method of the 4. order
%n is the number of steps
d=length(x0);%the dimension of the system (number of equations)
t=linspace(t0,tk,n);%independent variable
h=t(2)-t(1);%step size
x=zeros(d,n);%space for the solution
x(:,1)=x0;%initial value 
for i=1:n-1
	k1=h*feval(fun,t(i),x(:,i));
	k2=h*feval(fun,t(i)+h/2,x(:,i)+k1/2);
	k3=h*feval(fun,t(i)+h/2,x(:,i)+k2/2);
	k4=h*feval(fun,t(i)+h,x(:,i)+k3);
	x(:,i+1)=x(:,i)+(k1+2*k2+2*k3+k4)/6;
end