  In this article we develop a high order accurate method to solve the
incompressible boundary layer equations in a provably stable manner.~We first
derive continuous energy estimates,~and then proceed to the discrete
setting.~We formulate the discrete approximation using high-order finite
difference methods on summation-by-parts form and implement the boundary
conditions weakly using the simultaneous approximation term method.~By applying
the discrete energy method and imitating the continuous analysis,~the discrete
estimate that resembles the continuous counterpart is obtained proving
stability.~We also show that these newly derived boundary conditions removes
the singularities associated with the null-space of the nonlinear discrete
spatial operator.~Numerical experiments that verifies the high-order accuracy
of the scheme and coincides with the theoretical results are presented.~The
numerical results are compared with the well-known Blasius similarity solution
as well as that resulting from the solution of the incompressible Navier Stokes
equations.
