input_case nil$

n:=4$

% Nonlinear water waves
% ---------------------
% perturbation theory up to n-th order in the amplitude a
% F - velocity potential, Z - water level
% boundary conditions:
% sub(z=Z, df(F,t) + df(F,x)^2/2 + df(F,z)^2/2 + z ) = 0
% df(Z,t) = sub(z=Z, df(F,z) - df(F,x)*df(Z,x) )
% wave vector k=1, xi = x-v*t

operator b; weight a=1$ wtlevel n+1$

% general form of the velocity potential
F:=a*cos(xi)*e^z+
    for i:=4:n+1 sum a^i*
        if evenp(i)
        then for j:=2 step 2 until i sum b(i,j)*sin(j*xi)*e^(j*z)
        else for j:=3 step 2 until i sum b(i,j)*cos(j*xi)*e^(j*z)$
fz:=df(F,z)$

weight z=1$ % Taylor expansion

for all x let e^x=
begin scalar s,u,i; s:=u:=i:=1;
    while (u:=u*x/i) neq 0 do << s:=s+u; i:=i+1 >>;
    return s
end;

f:=F$ fz:=fz$ for all x clear e^x;
depend xi,x,t; let df(xi,x)=1,df(xi,t)=-v; fx:=df(f,x)$

for all x let cos(x)^2=(1+cos(2*x))/2,
              sin(x)^2=(1-cos(2*x))/2;
for all x,y let cos(x)*cos(y)=(cos(x-y)+cos(x+y))/2,
                sin(x)*sin(y)=(cos(x-y)-cos(x+y))/2,
                sin(x)*cos(y)=(sin(x-y)+sin(x+y))/2;

% water level from the boundary condition
Z:=-df(f,t)-fz^2/2-fx^2/2$

repeat <<clear z;Z1:=sub(z=z1,Z);weight z=1;Z:=sub(z1=Z,Z1)>>
until freeof(Z,z);

clear z;

% this expression should vanish due to the boundary condition
q:=df(Z,t)-sub(z=Z,fz)+sub(z=Z,fx)*df(Z,x)$

clear f,fz,fx; clear df(xi,x),df(xi,t);
for all x clear cos(x)^2,sin(x)^2;
for all x,y clear cos(x)*cos(y),sin(x)*sin(y),sin(x)*cos(y);

weight u=2$
let v^2=1+u; q:=q$ Z:=Z$ clear v^2;

let v=
begin scalar a,b,j,i0; a:=1; b:=3/2; j:=0;
    i0:=if evenp(n) then n/2 else (n-1)/2;
    return 1+for i:=1:i0 sum <<b:=b-1;j:=j+1;a:=a*b/j*u>>
end;

q:=q$ Z:=Z$ clear v;

% resonant term
let z*cos(xi)=1; r:=q*z$ clear z*cos(xi);
r:=sub(z=0,r)$ q:=q-r*cos(xi)$ U:=u-r/a$

repeat <<clear u;U1:=sub(u=u1,U);weight u=2;U:=sub(u1=U,U1)>>
until freeof(U,u);

clear u; u:=U$ clear U; q:=q$ Z:=Z$
for all j let cos(j*xi)=z^j,sin(j*xi)=z^j; q:=q$
for all j clear cos(j*xi),sin(j*xi); clear a;
qa:=coeff(q,a)$ clear q; qa:=rest(rest(rest(rest(qa))))$

for i:=4:n+1 do
<<  qi:=first(qa); qa:=rest(qa);
    j0:=if evenp(i) then 2 else 3;
    qi:=sub(z=sqrt(z),qi/z^j0); qz:=coeff(qi,z); clear qi;
    for j:=j0 step 2 until i do
    <<  qj:=first(qz); qz:=rest(qz);
        b(i,j):=-sub(b(i,j)=0,qj)/df(qj,b(i,j));
        clear qj;
    >>
>>;

factor a; on revpri,rat; % printing results
1+u;                     % omega^2
F;                       % velocity potential
Z;                       % water level

showtime;

end;