Comment: Check the slowly-varying normal form of thin film fluids.  
Make a near identity coordinate change: p=pp+..., v=vv+..., u=uu+...,
eta=hh+...  where pp=vv=0 and uu_t=guu and hh_t=ghh
Have put Weber and Reynolds numbers in explicitly.

Use a scaling where uu=O(d^2) and d/dx=O(d).

This is in Reduce, see http://www.reduce-algebra.com/

August 2004, Tony Roberts; 
% improves printing
on div; off allfac; on revpri;
factor d,re,we; 

% Use d to count the number of derivatives of x,
% and the number of uu's depending upon relative scaling
depend xs,x;
let df(~a,x) => d*df(a,xs);

% ii denotes \int_0^y ... dy
operator ii; linear ii;
let { ii(~a,y,1) => ii(a,y)
    , ii(~a,y,~n) => ii(ii(a,y),y,n-1) when n>1
    , ii(y*~a,y) => y*ii(a,y)-ii(ii(a,y),y)
    , ii(~a*y^~p,y) => y^p*ii(a,y)-ii(p*y^(p-1)*ii(a,y),y)
    , ii(1,y)=>y, ii(y,y)=>y^2/2, ii(y^~p,y)=>y^(p+1)/(p+1)
    };

% use operator hh(m) to denote df(hh,x,m)
operator hh;
depend hh,xs,t;
hs:=hh(0)$ hx:=hh(1)*d$ hxx:=hh(2)*d^2$ hxxx:=d^3*hh(3)$
let { df(hh(~m),xs) => hh(m+1)
    , df(hh(~m),t) => df(ghh,xs,m) };

% use operator uu(m,n) to denote df(uu,x,m,y,n)
operator uu;
depend uu,xs,y,t;
let { df(uu(~m,~n),xs) => uu(m+1,n)
    , df(uu(~m,~n),y) => uu(m,n+1)
    , df(uu(~m,~n),t) => df(guu,xs,m,y,n) when n>=0 
    , df(uu(~m,~n),t) => ii(df(guu,xs,m),y,-n) when n<=-1 
    , ii(uu(~m,~n),y) => uu(m,n-1)-bed(uu(m,n-1),y)
    };

% fs denotes evaluation on the psuedo-free surface y=hh
operator fs; linear fs;
let { fs(fs(~a,y),y)=>fs(a,y)
    , df(fs(~a,y),y)=>0
    , df(fs(~a,y),xs)=>fs(df(a,xs),y)+fs(df(a,y),y)*hh(1)
    , df(fs(~a,y),t)=>fs(df(a,t),y)+fs(df(a,y),y)*ghh
    , fs(1,y)=>1, fs(y,y)=>hs, fs(y^~p,y)=>hs^p 
    , fs(~a*~b,y)=>fs(a,y)*fs(b,y) 
    , ii(fs(~a,y),y)=>fs(a,y)*ii(1,y)
    , ii(~b*fs(~a,y),y)=>fs(a,y)*ii(b,y)
    };

% putative terms to help check orders of errors OK
depend etad,xs,t;
depend ghd,xs,t;
operator error;
depend error,xs,y,t;

% the PDE and BC on quasi-velocity uu give these FS relations
let { fs(uu(0,1),y)=>d^2*fs(uu(2,-1)+4*uu(1,0)*hh(1),y)
        +d^3*error(f)  
    , fs(uu(1,1),y)=>-fs(uu(0,2)*hh(1),y)
        +d^2*fs(uu(3,-1)+5*uu(2,0)*hh(1)+4*(-uu(0,2)*hh(1))*hh(1)^2
        +4*uu(1,0)*hh(2),y) +d^3*error(a) % from d/dx 
    , fs(uu(2,1),y)=>-fs(2*uu(1,2)*hh(1)+uu(0,3)*hh(1)^2
        +uu(0,2)*hh(2),y)+d^2*error(b) % from d^/dx^2 
    , fs(uu(3,1),y)=>-fs(3*uu(2,2)*hh(1)+3*uu(1,3)*hh(1)^2
        +uu(0,4)*hh(1)^3+2*uu(1,2)*hh(2)+2*uu(0,3)*hh(1)*hh(2)
        +uu(0,2)*hh(3),y) +d^2*error(c) % from d^3/dx^3 
    , fs(uu(0,3),y)=>d^2*error(d) % from d/dt 
    , fs(uu(1,3),y)=>-fs(uu(0,4)*hh(1),y)
        +d^2*error(e) % from d^2/dx/dt 
    };

% bed denotes evaluation on the substrate y=0
operator bed; linear bed;
let { bed(bed(~a,y),y)=>bed(a,y)
    , bed(fs(~a,y),y)=>fs(a,y)
    , fs(bed(~a,y),y)=>bed(a,y)
    , df(bed(~a,y),y)=>0
    , df(bed(~a,y),xs)=>bed(df(a,xs),y)
    , df(bed(~a,y),t)=>bed(df(a,t),y)
    , bed(1,y)=>1, bed(y,y)=>0, bed(y^~p,y)=>0 
    , bed(~a*~b,y)=>bed(a,y)*bed(b,y) 
    , bed(uu(~m,~n),y)=>0 when n<=0
    , ii(bed(~a,y),y)=>bed(a,y)*ii(1,y)
    , ii(~b*bed(~a,y),y)=>bed(a,y)*ii(b,y)
    };
 
% ifs=\int_0^hh ... dy
procedure ifs(a); fs(ii(a,y),y);

% scale the quasi-velocity with (d/dx)^2
us:=d^2*uu(0,0);

% specify the normal form as found to date
u:=us+we*(hs*y-y^2/2)*hxxx +d^5*error(uu);
v:=-ii(df(us,x),y)-we*df((hs*y^2/2-y^3/6)*hxxx,x) +d^6*error(vv);
p:=we*(-hxx+(3/2)*hx^2*hxx-(1/2)*df((hs^2+2*hs*y-y^2)*hxxx,x))
   -2*df(fs(us,y),x) +bed(df(df(us,y),x),y)*(hs-y) +d^5*error(pp);
eta:=hs+re*( df(ifs((hs*y-y^2/2)*us),x) 
+ifs((1/24*y^4-1/6*hs*y^3-1/4*hs^2*y^2+3/2*hs^3*y-1/8*hs^4)*df(us,x,3))
+df(hs,x)*ifs((-1/6*y^3-1/2*hs*y^2+13/2*hs^2*y-5/6*hs^3)*df(us,x,x))
+hs*df(hs,x)^2*ifs((4*y-hs)*df(us,x))
+hs^2*df(hs,x,x)*ifs((2*y-1/3*hs)*df(us,x))
+12*hs^3*df(hs,x)^2*fs(df(us,x),y)
+26/3*hs^3*df(hs,x)*df(hs,x,x)*fs(us,y)
+7*hs^2*df(hs,x)^3*fs(us,y)
+3*hs^4*df(hs,x)*fs(df(us,x,x),y)
+2*hs^4*df(hs,x,x)*fs(df(us,x),y)
+hs^4*df(hs,x,x,x)*fs(us,y)
-2*hs^4*df(hs,x)^3*fs(df(us,y,y),y)
)+d^6*etad;

% quasi-evolution, scaling is relative to uu
guu:=(1/d^2)*((df(us,y,2)+df(us,x,2))/re
    -((hs-y)*bed(df(us,x,x,y),y)+bed(df(us,x,y),y)*hx)/re
    -(-2*fs(df(us,x,x),y) +2*fs(df(us,y,y),y)*hx^2)/re
    +d^5*error(gu));
ghh:=-(we/3)*df(hs^3*hxxx,x) +d^6*ghd;

% continuity equation
write ceq:=(-df(u,x)-df(v,y) where d^6=>0);
write kin:=(fs(df(eta,t)-v+u*df(eta,x),y) where d^6=>0);
write con:=(-df(eta,t)-df(ifs(u),x) where d^6=>0);
% normal momentum & normal stress
write veq:=(re*( df(v,t)+u*df(v,x)+v*df(v,y) )
    +df(p,y) -df(v,x,2)-df(v,y,2) where d^5=>0);
curv:=hxx*(1-hx^2/2+3*hx^4/8-5*hx^6/16)$
write tn:= (fs(-p*(1+hx^2) +2*(df(v,y)
    +hx^2*df(u,x)-hx*(df(u,y)+df(v,x)))
    -we*curv ,y) where d^5=>0);
% lateral momentum & tangential stress
write ueq:=(re*( df(u,t)+u*df(u,x)+v*df(u,y) )
    +df(p,x) -df(u,x,2)-df(u,y,2) where d^5=>0);
write tt:=(fs( (1-hx^2)*(df(u,y)+df(v,x))
    +2*hx*(df(v,y)-df(u,x)) ,y) where d^5=>0);
showtime;
end;