function [t,x]=rk308h(fun,t0,x0,N,h)
f=fcnchk(fun,'t','x');
t=t0+(0:N)*h;
x=[x0,zeros(1,N)];
for nn=1:N,
    t_old=t(nn); x_old=x(nn);
    L1=feval(f,t_old,x_old);
    L2=feval(f,t_old+h/2,x_old +(h/2)*L1);
    L3=feval(f,t_old+h/2,x_old +(h/2)*L2);
    L4=feval(f,t_old+h,x_old + h*L3);
     x(nn+1)=x_old+(h/6)*(L1+2*L2+2*L3+L4);
end