function [t,x]=rk23(f,x0,t0,tf,reltol)
    abstol=1e-9;
    hmin=1e-9*(tf-t0);
    hmax=(tf-t0)/50;   
    h=hmin;
    N=2;
    t=zeros(1,100000);
    x=zeros(length(x0),length(t));
    t(1)=t0;
    x(:,1)=x0;
    k=1;
    while t(k)<tf
        [x1,x2]=rk23_step(f,x(:,k),t(k),h);
        e=norm(x1-x2)+1e-12;
        tol=reltol*norm(x1)+abstol;
        if e<tol then
            x(:,k+1)=x2;
            t(k+1)=t(k)+h;
            k=k+1;
        end
        h=0.8*h*(tol/e)^(1/(N+1));
        if h<hmin then 
            h=hmin;
        end
        if t(k)+h>tf then 
            h=tf-t(k);
        end
        if h>hmax then 
            h=hmax;
        end
        
    end
    t=t(1:k);
    x=x(:,1:k);
endfunction

function [x1,x2]=rk23_step(f,xk,tk,h)
    k1=f(xk,tk);
    k2=f(xk+h*k1,tk+h);
    k3=f(xk+h*k1/4+h*k2/4,tk+h/2);
    x1=xk+h*(k1/2+k2/2);
    x2=xk+h*(k1/6+k2/6+k3*2/3);
endfunction

    
            
    