function [t,x]=rk4(t0,tn,x0,f,h)

x(1) = x0;
t(1) = t0;
np = (tn-t0)/h;

for k=2:np+1
    t(k) = t(k-1) + h;
    k1 = f(t(k-1),x(k-1));
    k2 = f(t(k-1)+h/2, x(k-1) + (h/2)*k1);  
    k3 = f(t(k-1)+h/2, x(k-1) + (h/2)*k2);
    k4 = f(t(k-1)+h, x(k-1) + h*k3);
    x(k) = x(k-1) + (h/6)*(k1+2*k2+2*k3+k4);
end
    

