I am implementing Runge-Kutta method with adaptive step in matlab. I get different results as compared to matlab's own ode45 and my own implementation of Runge-Kutta method with fixed step. What am I doing wrong in my code? Is it possible?

function [ result ] = rk4_modh( f, int, init, h, h_min )
%
%   f - function handle
%   int - interval - pair (x_min, x_max)
%   init - initial conditions - pair (y1(0),y2(0))
%   h_min - lower limit for h (step length)
%   h - initial step length
%   x - independent variable ( for example time )
%   y - dependent variable - vertical vector - in our case ( y1, y2 )
    function [ k1, k2, k3, k4, ka, y ] = iteration( f, h, x, y )
    %   core functionality performed within loop
        k1 = h * f(x,y);
        k2 = h * f(x+h/2, y+k1/2);
        k3 = h * f(x+h/2, y+k2/2);
        k4 = h * f(x+h, y+k3);
        ka = (k1 + 2*k2 + 2*k3 + k4)/6;
        y = y + ka;
    end

    % constants
    % relative error
    eW = 1e-10;
    % absolute error
    eB = 1e-10;
    s = 0.9;
    b = 5;

    % initialization
    i = 1;
    x = int(1);
    y = init;
    while true
        hy = y; hx = x;

        %algorithm
        [ k1, k2, k3, k4, ka, y ] = iteration( f, h, x, y );
        % error estimation
        for j=1:2
            [ hk1, hk2, hk3, hk4, hka, hy ] = iteration( f, h/2, hx, hy );
            hx = hx + h/2;
        end
        err(:,i) = abs(hy - y);
        % step adjustment
        e = abs( hy ) * eW + eB;
        a = min( e ./ err(:,i) )^(0.2);
        mul = a * s;
        if mul >= 1
        % step length admitted
            keepH(i) = h;
            k(:,:,i) = [ k1, k2, k3, k4, ka ];
            previous(i,:) = [ x+h, y' ]; %'
            i = i + 1;
            if floor( x + h + eB ) == int(2)
                break;
            else
                h = min( [mul*h, b*h, int(2)-x] );
                x = x + keepH(i-1);
            end

        else
        % step length requires further adjustments
            h = mul * h;
            if ( h < h_min )
                error('Computation with given precision impossible');
            end
        end
    end
    result = struct( 'val', previous, 'k', k, 'err', err, 'h', keepH );
end

The function in question is:

function [ res ] = fun( x, y )
%   
    res(1) = y(2) + y(1) * ( 0.9 - y(1)^2 - y(2)^2 );
    res(2) = -y(1) + y(2) * ( 0.9 - y(1)^2 - y(2)^2 );
    res = res'; %'
end

The call is:

res = rk4( @fun, [0,20], [0.001; 0.001], 0.008 );

The resulting plot for x1 :

The result of ode45( @fun, [0, 20], [0.001, 0.001] ) is: