P(z) := (1 + z^(-1))^2*(1 + z)^2*(r1*z + r0 + r1*z^(-1));
series(P(z), z);

Should be  arctan  instead of  arctg  .

f:=(x,y)->piecewise(x=0 and y=0,139 ,x=1 and y=0,144,x=2 and y=0,149,x=3 and y=0,155,x=4 and y=0,155,x=5 and y=0,155,x=6 and y=0,155,x=7 and y=0,155,x=0 and y=1,144,x=1 and y=1,151,x=2 and y=1,153,x=3 and y=1,156,x=4 and y=1,159,x=5 and y=1,159,x=6 and y=1,159,x=7 and y=1,159,x=0 and y=2,150,x=1 and y=2,155,x=2 and y=2,160,x=3 and y=2,163,x=4 and y=2,158,x=5 and y=2,156,x=6 and y=2,156,x=7 and y=2,156,x=0 and y=3,159,x=1 and y=3,161,x=2 and y=3,162,x=3 and y=3,160,x=4 and y=3,160,x=5 and y=3,159,x=6 and y=3,159,x=7 and y=3,159,x=0 and y=4,159,x=1 and y=4,160,x=2 and y=4,161,x=3 and y=4,162,x=4 and y=4,162,x=5 and y=4,155,x=6 and y=4,155,x=7 and y=4,155,x=0 and y=5,161,x=1 and y=5,161,x=2 and y=5,161,x=3 and y=5,161,x=4 and y=5,160,x=5 and y=5,157,x=6 and y=5,157,x=7 and y=5,157,x=0 and y=6,162,x=1 and y=6,162,x=2 and y=6,161,x=3 and y=6,163,x=4 and y=6,162,x=5 and y=6,157,x=6 and y=6,157,x=7 and y=6,157,x=0 and y=7,162,x=1 and y=7,162,x=2 and y=7,161,x=3 and y=7,161,x=4 and y=7,163,x=5 and y=7,158,x=6 and y=7,158,x=7 and y=7,158):

C:=t->piecewise(t=0, 1/sqrt(2), t>0 or t<0, 1):

F:=(u,v)->1/4*C(u)*C(v)*add(add(f(x,y)*cos((2*x+1)*u*Pi/16)*cos((2*y+1)*v*Pi/16), y=0..7), x=0..7):

F(0,0);

                                                                             1261

Your equation has too many parameters. If you set their values​​, the equation can be solved exactly sometimes, but more often numerically as follows:

restart;

p:=a*sin(omega[P]*t): q:=b*sin(omega[Q]*t+phi[Q]): r:=c*sin(omega[R]*t+phi[R]): s:=d*sin(omega[S]*t+phi[S]):

z:=p+q+r+s:

a,b,c,d,omega[P],omega[Q],omega[R],omega[S],phi[Q],phi[R],phi[S]:=1,2,3,4, Pi/6,Pi/3,Pi/4,2*Pi/3, Pi/2,Pi/4,Pi/4:  # values of parameters

fsolve(z, t); 

               -0.6371104099

f:=a->fsolve({0.5704+0.5211*exp(-0.6569*x)-0.03849*y=a, 0.561+0.5284*exp(-0.6479*x)-0.03929*y=a}):

f(1);

A:=plot(a->rhs(f(a)[1]),2..10, labels=[a, x(a)]):

B:=plot(a->rhs(f(a)[2]),2..10, labels=[a, y(a)]):

plots[display](Array([A, B]));

f:=(x,y)->piecewise(x=0 and y=0,139 ,x=1 and y=0,144,x=2 and y=0,149,x=3 and y=0,155,x=4 and y=0,155,x=5 and y=0,155,x=6 and y=0,155,x=7 and y=0,155,x=0 and y=1,144,x=1 and y=1,151,x=2 and y=1,153,x=3 and y=1,156,x=4 and y=1,159,x=5 and y=1,159,x=6 and y=1,159,x=7 and y=1,159,x=0 and y=2,150,x=1 and y=2,155,x=2 and y=2,160,x=3 and y=2,163,x=4 and y=2,158,x=5 and y=2,156,x=6 and y=2,156,x=7 and y=2,156,x=0 and y=3,159,x=1 and y=3,161,x=2 and y=3,162,x=3 and y=3,160,x=4 and y=3,160,x=5 and y=3,159,x=6 and y=3,159,x=7 and y=3,159,x=0 and y=4,159,x=1 and y=4,160,x=2 and y=4,161,x=3 and y=4,162,x=4 and y=4,162,x=5 and y=4,155,x=6 and y=4,155,x=7 and y=4,155,x=0 and y=5,161,x=1 and y=5,161,x=2 and y=5,161,x=3 and y=5,161,x=4 and y=5,160,x=5 and y=5,157,x=6 and y=5,157,x=7 and y=5,157,x=0 and y=6,162,x=1 and y=6,162,x=2 and y=6,161,x=3 and y=6,163,x=4 and y=6,162,x=5 and y=6,157,x=6 and y=6,157,x=7 and y=6,157,x=0 and y=7,162,x=1 and y=7,162,x=2 and y=7,161,x=3 and y=7,161,x=4 and y=7,163,x=5 and y=7,158,x=6 and y=7,158,x=7 and y=7,158);

add(f(0, y), y = 0 .. 7); 

Suggest significantly faster way to solve the original problem. Code is designed as a procedure. Formal argument - the search range. Solutions based on the idea of ​​a two-stage approach.

1) First, in the specified range, we find vectors of integer length.

2) Next, we will find already the desired triangles of these vectors.

Because we use the sets, then repeats excluded. Output contains coordinates of the vertices of the triangles and also the lengths of the all sides.

IntegerTriangle:=proc(R::range)

local L, x, y, z, d, P, Q, T, a, b, c, p;

L:={}:

for x from lhs(R) to rhs(R) do

for y from lhs(R) to rhs(R) do

for z from lhs(R) to rhs(R) do

d:=sqrt(x^2+y^2+z^2):

if type(d,posint) then L:={op(L), [x,y,z,d]} fi:

od: od: od:

T:={}:

for P in L do

if P[1]<=P[2] and P[2]<=P[3] then

for Q in (L minus {P}) do

a:=P[4]: b:=Q[4]: c:=sqrt((P[1]-Q[1])^2+(P[2]-Q[2])^2+(P[3]-Q[3])^2):

p:=(a+b+c)/2:

if type(2*p, integer) and type(sqrt(p*(p-a)*(p-b)*(p-c)), posint) then

T:={op(T), {sort([a,b,c]), {[0,0,0], P[1..3], Q[1..3]}}} fi:

od: fi; od:

T;

end proc:

Examples.

Triangles in the original range:

t:=time():

L:=IntegerTriangle(-4..4):

nops(L);

L[i] $ i=31..40;

time()-t;

Integer triangles in the 1st octant:

t:=time():

L:=IntegerTriangle(0..30):

nops(L);

time()-t;

                       364

                     1.750

restart;

C:=proc(n, k)

option remember;

if n::even then C(n,0):=(-1)^(n/2)*n!/(n/2)!; C(n,1):=0 else

C(n,0):=0; C(n,1):=2*(-1)^((n-1)/2)*n!/((n-1)/2)! fi;

2*(k-2-n)/(k*(k-1))*C(n,k-2);

end proc:

Examples:

C(10,5);  C(11,7);

                  0

             506880

I think that with your statement of the problem, in principle, there is no single solution. The fact the condition  D(f)(infinity) =0  will be satisfied if  limit(f(eta), eta=infinity)=Const . However, Const  may be different, ie solutions will differ. In reality you are interested in solutions on a finite interval  eta=0..T . Assume that for  eta>T  the function does not change substantially. Therefore, you solve the equation for this interval, and for the values   eta>T  you consider  f(eta)=a  where  a=f(T) . Let  T=20 .

restart;

eq1 := diff(f(eta), eta, eta, eta) = (diff(f(eta), eta))^2-f(eta)*(diff(f(eta), eta, eta));

bcs1 := (D(f))(0) = 1, f(0) = 0, (D(f))(20) = 0:

sol1 := dsolve({bcs1, eq1}, numeric):

g := proc(eta)

rhs(sol1(eta)[2]) end proc:

a := g(20);

f := eta->piecewise(0 <= eta and eta <= 20, g(eta), 20 < eta, a);

plot(f, 0 .. 40, labels=[eta, 'f(eta)']);

The procedure  MinAreaRect  returns a list that contains the smallest rectangle and the output of  DirectSearch[Optima]  command (the area of this rectangle, etc.).

MinAreaRect:=proc(P)

local f, Left, Right, Bottom, Top, Rot, Area, phi0, n, k, X, Y, PosXmin, PosXmax, PosYmin, PosYmax, P1, P2, P3, P4;

uses ListTools;

f:=(x,y)->[x*cos(phi)-y*sin(phi), x*sin(phi)+y*cos(phi)];

Left:=L->min(seq(L[i,1], i=1..nops(L)));

Right:=L->max(seq(L[i,1], i=1..nops(L)));

Bottom:=L->min(seq(L[i,2], i=1..nops(L)));

Top:=L->max(seq(L[i,2], i=1..nops(L)));

Rot:=L->map(t->f(op(t)), L);

Area:=DirectSearch[GlobalOptima]((Right(Rot(P))-Left(Rot(P)))*(Top(Rot(P))-Bottom(Rot(P))), [phi=-Pi/2..Pi/2]);

phi0:=rhs(op(Area[2]));

n:=nops(P); k:=-tan(phi0);

X:=[seq(solve(-P[i,2]=k*(x-P[i,1])), i=1..n)];

PosXmin:=Search(min(X), X);

PosXmax:=Search(max(X), X);

Y:=[seq(P[i,2]-1/k*(-P[i,1]), i=1..n)];

PosYmin:=Search(min(Y), Y);

PosYmax:=Search(max(Y), Y);

P1:=eval([x,y], solve({y-P[PosXmin,2]=k*(x-P[PosXmin,1]), y-P[PosYmin,2]=-1/k*(x-P[PosYmin,1])}));

P2:=eval([x,y], solve({y-P[PosYmin,2]=-1/k*(x-P[PosYmin,1]), y-P[PosXmax,2]=k*(x-P[PosXmax,1])}));

P3:=eval([x,y], solve({y-P[PosXmax,2]=k*(x-P[PosXmax,1]), y-P[PosYmax,2]=-1/k*(x-P[PosYmax,1])}));

P4:=eval([x,y], solve({y-P[PosYmax,2]=-1/k*(x-P[PosYmax,1]), y-P[PosXmin,2]=k*(x-P[PosXmin,1])}));

[[P1,P2,P3,P4], Area];

end proc: 

Examples:

Model example. It shows that the statement "The rectangle has a side parallel to the longest diagonal of P "  is false.

P:=[[0,1],[1,0],[0,-1],[-1,0]]:
MinAreaRect(P);

[[[-1.000000000, .6582500000e-12], [-.1316500000e-11, -1.000000000], [1., 0.], [.6582500000e-12, 1.000000000]], [2.00000000000263, [phi = -.78539816339679], 143]]

Original example:

P:=[[0,2],[1,4],[2,3.5],[4,4],[5,1],[4,0.75],[3,0]]:

MinAreaRect(P);

plots[display](plottools[polygon](P), plottools[polygon](%[1], style=line, color=green, thickness=2), scaling=constrained);

h1:=0.2: h2:=0.3:
plot([eval(eq1,h=h1), eval(eq2,h=h2)], D=0.5..1);

Slight improvement of the procedure  Location  can also work with simplex lying in spaces of higher dimension, such as a segment in the plane  R^2  or a triangle or a segment in space  R^3, etc.

Location:=proc(T::set(list), P::list)  # T - a simplex, P - a point

local m, n, lambda, Sol;

m:=nops(T); n:=nops(P);

Sol:=solve({seq(add(lambda[i]*T[i,j], i=1..m)=P[j], j=1..n), add(lambda[i], i=1..m)=1});

   if Sol<>NULL then assign(Sol);

   if convert([seq(lambda[i]>0 and lambda[i]<1, i=1..m)], `and`) then return `Inside` else

   if mul(lambda[i], i=1..m)=0 and convert([seq(lambda[i]>=0 and lambda[i]<=1, i=1..m)],         `and`)

then return `At the Border`  else return `Outside` fi; fi;  else `Outside` fi;

end proc:

Examples: 

Location({[1,0,0], [0,1,0]}, [1/2,1/2,0]);

Location({[1,0,0], [0,1,0], [0,0,1]}, [1/3,1/3,1/3]);

Location({[1,0,0], [0,1,0], [0,0,1]}, [-1,-1,3]);

Location({[1,0,0], [0,1,0], [0,0,1]}, [0,0,0]);

Location({[1,0,0], [0,1,0], [0,0,1]}, [1/2,1/2,0]);

                             Inside

                             Inside

                            Outside

                            Outside

                       At the Border

Markiyan!  Your idea about barycentric coordinates is good! I used it in the procedure  Location . This procedure solves the original problem for arbitrary dimension (for a segment, a triangle, a tetrahedron and so on).

restart;

Location:=proc(T::set(list), P::list)  # T - a simplex, P - a point

local m, n, lambda;

m:=nops(T); n:=nops(P);

solve({seq(add(lambda[i]*T[i,j], i=1..m)=P[j], j=1..n), add(lambda[i], i=1..m)=1});

assign(%);

   if convert([seq(lambda[i]>0 and lambda[i]<1, i=1..m)], `and`) then return `Inside` else

   if mul(lambda[i], i=1..m)=0 and convert([seq(lambda[i]>=0 and lambda[i]<=1, i=1..m)],    `and`) then return `At the Border` else `Outside` fi; fi;

end proc:

Examples:

Location({[0,0,0], [1,0,0], [0,1,0], [0,0,1]}, [1/4,1/4,1/4]);

Location({[0,0,0], [1,0,0], [0,1,0], [0,0,1]}, [1/3,1/3,1/3]);

Location({[0,0,0], [1,0,0], [0,1,0], [0,0,1]}, [1,1,1]);

                                          Inside
                                    At the Border
                                         Outside

Of course, for pedagogical purposes, my algorithm is incomplete and can only check whether the point of the interior of the tetrahedron.

If the point lie on the surface of the tetrahedron, for example on the face  ABC , then verify it:

1) Solve the system of  x * A + y * B + z * C = P .

2) If  x + y + z = 1  and  0 <= x <= 1  and  0 <= y <= 1  and  0 <= z <= 1  then the point  P  belongs to the face ABC .

3) If the point does not lie within the tetrahedron and does not lie on it's surface, it lies outside of the tetrahedron.

1) Find the equations of planes bounding the tetrahedron. Let they be F1, F2, F3, F4

2) Take a point lying exactly inside the tetrahedron. It be  K((x1+x2+x3+x4)/4, (y1+y2+y3+y4)/4, (z1+z2+z3+z4)/4)=K(x5, y5, z5)

3) If all the equalities hold  sign(F1(x0, y0, z0))=sign(F1(x5, y5, z5)),   sign(F2(x0, y0, z0))=sign(F2(x5, y5, z5)),  sign(F3(x0, y0, z0))=sign(F3(x5, y5, z5)),  sign(F4(x0, y0, z0))=sign(F4(x5, y5, z5)) , then the point  P(x0,y0,z0)  lies  inside tetrahedron ABCD .