A two-dimensional integral boundary layer method is developed to enable fast and economical computations of boundary layer flows. The ultimate goal is to provide some experience for the extension of this method in three dimensions. In this study, the unsteady momentum and kinetic energy integral equations are solved numerically, together with a set of closure relations based on assumed velocity profiles for laminar and turbulent flows. The robustness of the method is ensured by a Finite-Volume formulation based on an upwind scheme and a semi-implicit time discretization. The accuracy of the numerical method in the vicinity of the stagnation point is strongly improved by introducing a consistent corrective source term in the right-hand side of the equation system. The chosen closure relations are validated with test cases of self-similar flows. Numerical results are also compared with those of a full Prandtl equations code for NACA0012, GLC305 and MS317 airfoils test cases to demonstrate the capabilities of the method. Finally, preliminary results are shown proving the ability of the method to deal with iced airfoils even for complex glaze ice shapes.