The task is easy to solve by using Banach fixed-point theorem (contraction mapping principle). See  here

Consider the function   f:=x->1/2*(x+2/x) . The segment  [1, 2]  is a complete metric space in standard metric. 

Enough to check two conditions:
1) The function  f  maps the segment  [1, 2]  into itself.
2) The mapping  f  is a contraction.

This is easily done manually or with Maple.

Checking with Maple:

f:=x->1/2*(x+2/x):

is(f(x)>=1 and f(x)<=2) assuming x>=1,x<=2;  # Checking of first condition

abs('f'(a)-'f'(b))<=maximize(abs(diff(f(x),x)), x=1..2)*abs(a-b);  # Checking of second condition

The limit of the sequence  {u(n+1)=1/2*(u(n)+2/u(n)), u(0)=1}  is the fixed point of mapping  f  ;

solve({x=f(x), x>=1, x<=2});

Procedure  P  solves your problem for any number of couples  (m is the number of couples).

N:=proc(m::posint)

local P, C1, C, P1;

P:=(2*m)!;  # Total number of variants

C1:=2*binomial(2*m-1,1)*(2*m-2)!;  # Number of variants in which the members of the first couple are close

C:=k->binomial(2*m-k,k)*2^k*k!*(2*m-2*k)!;  # C(k) is the number of variants in which the members of the first  k couples are close

P1:=binomial(m,1)*C1+add(binomial(m,k)*(-1)^(k-1)*C(k), k=2..m); # Number of variants in which the members at least one couple are close (the exclusion-inclusion principle)

P-P1;

end proc:

Examples:

seq(N(m), m=1..20);  # The last number is the answer to your problem

Addition: independent check by the brute force method for  m=3 :

N:=0:  # {a,b}, {c,d}, {e,f} are three couples

for i in combinat[permute]([a,b,c,d,e,f]) do

if {i[1],i[2]}<>{a,b} and {i[1],i[2]}<>{c,d} and {i[1],i[2]}<>{e,f} and

{i[2],i[3]}<>{a,b} and {i[2],i[3]}<>{c,d} and {i[2],i[3]}<>{e,f} and

{i[3],i[4]}<>{a,b} and {i[3],i[4]}<>{c,d} and {i[3],i[4]}<>{e,f} and

{i[4],i[5]}<>{a,b} and {i[4],i[5]}<>{c,d} and {i[4],i[5]}<>{e,f} and

{i[5],i[6]}<>{a,b} and {i[5],i[6]}<>{c,d} and {i[5],i[6]}<>{e,f} then N:=N+1 fi:

od:

N;

                                                             240

assume(a::realcons, b::realcons):

eq:=I*(a+3*b)+2*b+5*a = 3+2*I:

solve({Re(lhs(eq))=Re(rhs(eq)), Im(lhs(eq))=Im(rhs(eq))});

                              {a = 5/13, b = 7/13}

If you want to use your rule to specific real numbers, then  write  n::realcons .

In general case write  n::anything .

Examples:

applyrule( n::realcons*F(K,L) = F(n*K,n*L), Pi*F(K,L) );

applyrule( n::anything*F(K,L) = F(n*K,n*L), k*F(K,L) );

                                   F(Pi*K, Pi*L)

                                    F(k*K, k*L)

Example 1. 

Find the least positive root of the equation  tan(x)+a*x=0  for  a  in the list  [0.1, 0.2, ...,0.9, 1]  (total 10 numbers).

restart;

map(fsolve, [seq(tan(x)+a*x=0, a=[seq(0.1*n, n=1..10)])], x=Pi/2..3*Pi/2);

assign(seq(x||i=%[i], i=1..10));

x1, x3;

   [2.862772588, 2.653662400, 2.498399006, 2.380644485, 2.288929728, 2.215707303, 2.155983633, 2.106375674, 2.064528178, 2.028757838]

                                                     2.862772588, 2.498399006

Example 2.

 Find the first  10  positive root of the equation  tan(x)+0.1*x=0 . As the parameter we mean the range in which we look for the root.

restart;

map2(fsolve, tan(x)+0.1*x=0, [seq(x=-Pi/2+Pi*n..Pi/2+Pi*n, n=1..10)]);

assign(seq(x||i=%[i], i=1..10));

x2, x5;

[2.862772588, 5.760557933, 8.708313831, 11.70267808, 14.73347234, 17.79083538, 20.86723817, 23.95736763, 27.05755029, 30.16523655]

                                                      5.760557933, 14.73347234

We have used the obvious fact that in each range  x=-Pi/2+Pi*n .. Pi/2+Pi*n  (n is integer)  the equation  tan(x)+a*x=0  (a>0)    has a unique root.

Default Maple operates with angles expressed in radians:

convert((58+1/60*28+1/60^2*18)*degrees, radians);

                                 (35083/108000)*Pi

Sys1 := {diff(r(t),t) = r(t)^2*sin(theta(t)), diff(theta(t),t) = -r(t)^2*(-2*cos(theta(t))^2+1)}

DEtools[phaseportrait](Sys1, [r(t),theta(t)],t=-1..1,[[r(0)=1,theta(0)=0]]);

See  LinearAlgebra[KroneckerProduct]  command.

restart;

Sys:={diff(x(t),t)=x(t)*y(t)-x(t)^2*y(t)+y(t)^3, diff(y(t),t)=y(t)^2+x(t)^3-x(t)*y(t)^2}:

trans := {x(t) = r(t)*cos(theta(t)), y(t) = r(t)*sin(theta(t))}:

eval(Sys, trans);

solve(%, {diff(r(t),t),diff(theta(t),t)});

simplify((-1)^k)  assuming k::odd;

                       -1

2  and  3  are now symbols, not integers:

expand(exp(`2`*a+b)+exp(`3`*b+b));

In Maple 12:

T:=unapply(sum( sum(binomial(m,k+p), k=1..m-p)*x^p, p=0..m-1), x,m);

 T(x,3);

 simplify(%);

or

T:=(x,m)->sum( sum(binomial(m,k+p),k=1..m-p)*x^p,p=0..m-1);

eval(T(x, m), m=3);

simplify(%);

See information about  OrthogonalExpansions  package on  http://www.maplesoft.com/applications/view.aspx?SID=33406 

You can use a simple syntax for writing a triple integral. The example can be easily solved by direct calculation in Cartesian coordinates:

int(y, [y=0 .. sqrt(1-x^2-z^2), z=-sqrt(1-x^2) .. sqrt(1-x^2), x=-1 .. 1]); 

Solution of your example:

DEtools[convertsys]({diff(y(x), x, x)+2*y(x)+y(x)^2 = 0}, {y(0) = y0, (D(y))(0) = y1}, {y(x)}, x);