A stable numerical method for the solution of the adsorption integral equation (SAIEUS) is proposed. The method combines regularization principle with the B-spline representation of the distribution function. Such a representation is convenient and in the case of smooth functions requires fewer variables as compared to the discreet representation. The stability of the solution is imposed by the regularization method. The inherent problem of the method, which is the choice of the optimal degree of smoothing, is solved using a comprehensive analysis of the variance of the solution and the effective bias due to the smoothing effect. This approach ensures that the chosen solution contains all the information which can be extracted from the data while the artifacts are excluded. The problem of the resolution of the method and the accuracy of the solution in context of the error in the data and the complexity of the true adsorption energy distribution are discussed using several simulated examples of adsorption isotherms.