function [x,y]=rk4(dfun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1)=y0;
for k=1:length(x)-1
k1=feval(dfun,x(k),y(k));
k2=feval(dfun,x(k)+h/2,y(k)+h/2*k1);
k3=feval(dfun,x(k)+h/2,y(k)+h/2*k2);
k4=feval(dfun,x(k+1),y(k)+h*k3);
y(k+1)=y(k)+h*(k1+2*k2+2*k3+k4)/6;
end
x=x';y=y';