function varargout = stepest(varargin)
    [lhs,rhs] = argn(0)
    //checking the number of inputs
    if rhs > 2  then
        error(msprintf(gettext("%s: Unexpected number of input arguments "),"stepest"))
    end
    modelData = varargin(1)
    if typeof(modelData) <> "idpoly" then
        error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"stepest"))
    end
    //adding noise
    if rhs == 2 then
        noiseFlag = varargin(2)
        if typeof(noiseFlag) <> 'boolean' then
            error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"stepest"))
        end
    else
        noiseFlag = %F
    end
    z = poly(0,'z')
    aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff')
    bPoly = poly(modelData.b,'z','coeff')
    fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff')
    afPoly = aPoly*fPoly
    
    bCoeff = modelData.b
    extra = 1
    if ~bCoeff(1,1) then
        afCoeff = coeff(afPoly)
        bLength = length(bCoeff);afLength = length(afCoeff)
        if bLength == afLength then
            extra = 1
        else
            extra = z^-(bLength-afLength)
        end
    end
    
    bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff')
    sys = syslin('d',bPoly,afPoly)*extra
    if size(find(gsort(abs(roots(sys.den)))>1),'*') then
        tempRoot = clean(roots(sys.den))
        [sorted index] = gsort(abs(real(tempRoot)))
        tempRoot = tempRoot(index)
        n = 26/log10(abs(real(tempRoot(1))))
         
    else
        finalValue = horner(sys,1)*1
        n = 10
        tempValue = sum(ldiv(sys.num,sys.den,n))
        if finalValue >= 0 then
            while finalValue*0.99 >= tempValue || n >=1000
                n = n+10
                tempValue = sum(ldiv(sys.num,sys.den,n))
            end
        elseif finalValue < 0 then
            while finalValue*0.99 <tempValue || n >=1000
                n = n+10
                tempValue = sum(ldiv(sys.num,sys.den,n))
            end  
        end
    end
    uData = ones(1,n)
    yData = flts(uData,sys)
    timeData = (0:(n-1))*modelData.Ts
    if noiseFlag then
        varargout(1) = yData'
    else
        plot2d2(timeData,yData,2)
        varargout(1) = []
    end
endfunction