Example 4.4-1
Boundary Layer "Theory" applied to flow near the Leading edge of a flat Plate.
>  restart;
>  econ:=D[1](vx)(x,y)+D[2](vy)(x,y)=0; Continuity Equation: 4.4-1 & 9
>  emotx:=(x,y)->vx(x,y)*D[1](vx)(x,y)+vy(x,y)*D[2](vx)(x,y)+D(P)(x)/rho-nu*D[2](D[2](vx))(x,y); x component of eq. of motion (4.4-2) neglecting the second derivative of vx wrt x. This agrees with eq. 4.4-10.
>  vy:=(x,y)->-int(D[1](vx)(x,y1),y1=0...y); This should be a solution to the continuity equation.

>  econ; It is!

>  BE:=P(x)+rho*ve(x)^2/2; Bernoulli Equation: 4.3-5 for potential flow and neglecting vy in the potential flow region.
>  P(x):=solve(BE,P(x));P:=unapply(P(x),x);
>  emotx(x,y); This gives 4.4-12.

>  vx:=(x,y)->vinf*phi(y/delta(x)); The assumed form of solution suggested in eq. 4.4-14

Here is eq. 14 after this substitution.
>  emotx(x,y);
>  eq19:=int(emotx(x,y),y=0...delta(x)); Now integrate wrt y from 0 to delta(x).
>  phi:=eta->(3/2)*eta-(1/2)*eta^3; One possible polynomial that vanishes at eta=0 and has a zero derivative at eta = 1. Lots of others could be used. Suggested in eq. 4.4-14
>  ve:=x->vinf;eq19; This is a DE that determines the variation of delta with x.
>  s:=dsolve({eq19=0,delta(0)=0},delta(x)); Solving for delta and using the BC that it should be zero at x=0.

>  assign(s[1]); This checks with eq. 4.4-25
>  evalf(delta(x)); So does this.