function [x,t]=rk4(x,t,h,fun)
  w1=feval(fun,x,t);
  w2=x+(h/2)*w1;
  w3=x+(h/6)*w1;
  w1=feval(fun,w2,t+h/2);
  w2=x+(h/2)*w1;
  w3=w3+(h/3)*w1;
  w1=feval(fun,w2,t+h/2);
  w2=x+h*w1;
  w3=w3+(h/3)*w1;
  w1=feval(fun,w2,t+h);  
  t=t+h;
  x=w3+(h/6)*w1;