Unsteady boundary layer equations are solved by a noniterative implicit numerical scheme which is second-order accurate both in time and space. The present numerical method provides initial data within the method itself by solving a reduced initial boundary value problem applicable at or very near the origin of a streamwise coordinate. The method has been applied to unsteady laminar boundary layers with large disturbance parameters such as transition to Falkner-Skan flow, oscillating Blasius flow, oscillating Stagnation-Point flow, constant accelerated flow and harmonically fluctuating boundary layer past a circular cylinder. The thermal effects are also considered in the boundary layers of the last application. Comparison with the existing data has demonstrated excellency of the present method both in accuracy and computer-time economy.