function [t,x] = rk4m(f,tspan,x0,h),
a = tspan(1);
b = tspan(2);
n = (b-a)/h;
t = (a+h:h:b);
k1 = h*feval(f,a,x0);
k2 = h*feval(f,a+h/2,x0+k1/2);
k3 = h*feval(f,a+h/2,x0+k2/2);
k4 = h*feval(f,a+h,x0+k3);
x(1) = x0 + (k1)/6+ (k2)/3 + (k3)/3 + (k4)/6;
for i=1:n-1
k1 = h*feval(f,t(i),x(i));
k2 = h*feval(f,t(i)+h/2,x(i)+k1/2);
k3 = h*feval(f,t(i)+h/2,x(i)+k2/2);
k4 = h*feval(f,t(i)+h,x(i)+k3);
x(i+1) = x(i) + (k1)/6+ (k2)/3 + (k3)/3 + (k4)/6;
end
t = [a t];
x = [x0 x];
