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 .

Your equation  5 = floor(4969/x)+floor(208/x)  can be solved exactly by using  simple estimates with a further solution in Maple. If  x <= 800  or  x> = 1000, then obviously no solution because in the first case    floor(4969/x)>=6  and in the second case  floor(4969/x)<5  and  floor(208/x)=0 . But if  800<x<1000  
then floor(208/x)=0  and the equation   5 = floor(4969/x)+floor(208/x)  is equivalent to equation 5 = floor(4969/x) .

solve((5 <= 4969/x and 4969/x < 6);  #  Exact solution 

                     RealRange(Open(4969/6),4969/5)

evalf(%);

                RealRange(Open(828.1666667), 993.8000000)

You can use my procedure  Picture .  See  http://www.mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-

Example (shaded circular segment: center = [2,1], radius= 3, from  -(1/3)*Pi  to Pi radians):

L := [[[2+3*cos(t), 1+3*sin(t)], t = -(1/3)*Pi .. Pi], [[-1, 1], [7/2, 1-3*sqrt(3)*(1/2)]]]:
plots[display](Picture(L, color = green, [color = black, thickness = 2]), scaling = constrained, view = [-1 .. 6, -2 .. 5]);

Another example:

L := [[3, t = -(1/2)*Pi .. (1/2)*Pi, polar], [[-4+5*cos(t), 5*sin(t)], t = arctan(3/4) .. -arctan(3/4)]]:
plots[display](Picture(L, color = yellow, [color = brown, thickness = 3]), scaling = constrained, view = [-1 .. 4, -4 .. 4]);

Your equation corresponds to the geometric image - the surface of the second order. If this surface is bounded on at least one dimension, the procedure  Isolve  finds all integer points on the surface.

Formal arguments of procedure:  Equation  is the equation of the surface,  Range  is the range of variation of a variable in the format  variable = a .. b. The easiest way to check boundness of the surface is the construction of it's plot. If the surface is unbounded, the procedure finds all the solutions in a specified range.

restart;

Isolve:=proc(Equation, Range)

local L, i, Sol;

L:=[]:

   for i from lhs(rhs(Range)) to rhs(rhs(Range)) do

      Sol:=isolve(eval(Equation, lhs(Range)=i));

      L:=[op(L),seq([op(Sol[k]),z=i], k=1..nops([Sol]))];

   od:

L;

end proc:

Your example:

plots[implicitplot3d](x^2 +y^2 +2*z^2 +x*y -y*z=3, x=-3..3, y=-3..3, z=-3..3, axes=normal, numpoints=10000);

+

We have an ellipsoid.

Isolve(x^2 +y^2 +2*z^2 +x*y -y*z=3, z=-2..2);

[[x = -1, y = 0, z = -1], [x = 1, y = -2, z = -1], [x = 1, y = 0, z = -1], [x = -2, y = 1, z = 0], [x = -1, y = -1, z = 0], [x = -1, y = 2, z = 0], [x = 1, y = -2, z = 0], [x = 1, y = 1, z = 0], [x = 2, y = -1, z = 0], [x = -1, y = 0, z = 1], [x = -1, y = 2, z = 1], [x = 1, y = 0, z = 1]]

Finding of the all integer points on unbounded surface  x^2+y^2-z^2=1  (hyperboloid) in the range  z=-100..100 :

L:=Isolve(x^2 +y^2 -z^2=1, z=-100..100):
nops(L);    # The total number of solutions
 L[1000..1020];     # A few solutions

                                                                                  2892

[[x = 31, y = 22, z = -38], [x = 34, y = -17, z = -38], [x = 34, y = 17, z = -38], [x = 38, y = -1, z = -38], [x = 38, y = 1, z = -38], [x = -37, y = -1, z = -37], [x = -37, y = 1, z = -37], [x = -29, y = -23, z = -37], [x = -29, y = 23, z = -37], [x = -23, y = -29, z = -37], [x = -23, y = 29, z = -37], [x = -1, y = -37, z = -37], [x = -1, y = 37, z = -37], [x = 1, y = -37, z = -37], [x = 1, y = 37, z = -37], [x = 23, y = -29, z = -37], [x = 23, y = 29, z = -37], [x = 29, y = -23, z = -37], [x = 29, y = 23, z = -37], [x = 37, y = -1, z = -37], [x = 37, y = 1, z = -37]]

I think that finding  all  integer points on the unbounded surface of the second order, in general case,  is not easy problem. For example,  Mathematica does not cope with this task.

We use the formula from wiki  http://en.wikipedia.org/wiki/Surface_of_revolution

sol:=solve(eval(x^2/20+y^2/13+z^2/20=1, z=0),  x);

2*Pi*int(sol[1]*sqrt(1+diff(sol[1], y)^2), y=-sqrt(13)..sqrt(13));

S:=combine(%);

evalf(S);

You need  an other  software, for example  

http://umsolver.com/usa/mathematics/free-algebra-equation-solver/ 

interface(rtablesize=infinity):

MyHandler := proc(operator,operands,default_value)

   NumericStatus( division_by_zero = false );

   return infinity;

end proc: 

NumericEventHandler(division_by_zero=MyHandler):

Matrix([[Rad, Degree, sin(phi), cos(phi), tan(phi)], [``,``,``,``,``], seq([seq(j, j=[i, i*180/Pi, convert(sin(i), radical), convert(cos(i), radical), convert(tan(i), radical)])], i=[seq(k*Pi/12, k=0..24)])]);

Fragment of output:

Acer's procedure  MyHandler  was used.

Exact periodic extension of function  g  to the whole real axis can be constructed as follows:

g := t->piecewise(0 <= t and t <= 0.2e-1, 122.63, 0.2e-1 < t and t < .287, 0, .287 <= t and t <= .307, 122.63, .307 < t and t < .335, 0):

R := t->(t/0.335-floor(t/0.335))*0.335:

G := unapply((g@R)(t), t);

Example:

plot(G(t), t = -2*.335 .. 3*.335);