function [TOUT,YOUT]=rungekutta(ODEFUN,A,B,Y0,OPT)
    //adaptive Runge-Kutta (4th order) ODE solver with error step control
    // ODEFUN - ODE in form y'=ODEFUN(t,y)
    //          where t is independent variable
    //          ODEFUN must return column vector
    // A      - t0
    // B      - tFinal
    // Y0     - initial conditions for ODEFUN, column vector
    // OPT    - options for solver, vector
    //          OPT(1) - eps, error tolerance
    //          OPT(2) - hmax, max step allowed
    //          OPT(3) - h0, first step tried
    //          if not passed, default values = [1e-4 Inf 1e-6]

    if argn(2)<5 then
        OPT = [1e-4 %inf 1e-6];
        if argn(2)<4 then
            disp('rungekutta error: not enought input arguments')
            error(39)
        end
    end

    //test for valid ODEFUN
    [nr,nc]=size(ODEFUN(A,Y0));
    if nr>1 then
        disp('rungekutta error: ODEFUN must return column vector')
        error(15)
    end

    //test for valid Y0
    [nr,nc]=size(Y0);
    if nr>1 then
        disp('rungekutta error: Y0 must be column vector')
        error(15)
    end

    EPS  = OPT(1);
    HMAX = OPT(2);
    h    = OPT(3);

    y = Y0;
    x = A;

    Nsol = 1e4;

    TOUT = zeros(Nsol,1);
    YOUT = zeros(Nsol,nc);

    TOUT(1)   = A;
    YOUT(1,:) = Y0;

    isDone = %F
    count  = 1
    while ~isDone

        count = count+1

        k1=h*ODEFUN(x,y);
        k2=h*ODEFUN(x+h/2,y+k1/2);
        k3=h*ODEFUN(x+h/2,y+k2/2);
        k4=h*ODEFUN(x+h,y+k3);

        y1=y+(k1+2*k2+2*k3+k4)/6;

        isAcceptable = %F
        while ~isAcceptable

            k2=h*ODEFUN(x+h/4,y+k1/4)/2;
            k3=h*ODEFUN(x+h/4,y+k2/2)/2;
            k4=h*ODEFUN(x+h/2,y+k3)/2;

            y_=y+(k1/2+2*k2+2*k3+k4)/6;

            k1_=h*ODEFUN(x+h/2,y_)/2;
            k2 =h*ODEFUN(x+3*h/4,y_+k1_/2)/2;
            k3 =h*ODEFUN(x+3*h/4,y_+k2/2)/2;
            k4 =h*ODEFUN(x+h,y_+k3)/2;

            y2=y_+(k1_+2*k2+2*k3+k4)/6;

            y__=(16*y2-y1)/15;

            isAcceptable = and(abs(y2-y__)<EPS);
            if isAcceptable then
                isPrecise = and(abs(y2-y__)<(0.8*EPS/32));
                if isPrecise then
                    h_Next=2*h;
                else
                    h_Next=h;
                end
            else
                h=h/2;
                y1=y_;
                k1=k1/2;
            end
        end

        x=x+h;
        y=y__;

        if h_Next>HMAX then
            h_Next = HMAX;
        end
        h=h_Next;

        isDone = abs(B-x)<1e-6;
        if ~isDone then
            isNearFinal = abs(B-x)<h;
            if isNearFinal then
                h=B-x;
            end
        end

        TOUT(count)   = x;
        YOUT(count,:) = y;
    end

    TOUT = TOUT(1:count)
    YOUT = YOUT(1:count,:)

endfunction
