If a function was specified by an expression, as the questioner did, then

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

diff(f(x), x);

simplify(%);

The result is the same.

restart;

S:={$0..9}: L:=[]:

for C in (S minus {0}) do

for N in (S minus {C}) do

for R in (S minus {C,N}) do

for P in (S minus {C,N,R}) do

for K in (S minus {0,C,N,R,P}) do

for H in (S minus {C,N,R,P,K}) do

for M in (S minus {C,N,R,P,K,H}) do

for E in (S minus {C,N,R,P,K,H,M}) do

a:=C*10^3+N*10^2+R*10+P:  b:=C*10^2+P*10+P:  c:=K*10^4+H*10^3+M*10^2+E*10+R:

if a + b = c  then  L:=[op(L), [CNRP=a, CPP=b, KHMER=c]]  fi:

od: od: od: od: od: od: od: od:

op(L);

Received a unique solution.               

The first term is already linear in  f .

expression := diff(f(r), r$2) + exp(-f(r))*(1 - g(r))^2:

op(1,expression)+subs({f=f(r), g=g(r)},mtaylor(eval(op(2,expression), {f(r)=f, g(r)=g}), [f, g], 2));

Let  in the pocket  a is the number of  50c coins,  b  is the number of  20c coins,  c  is the number of  10c coins,  d  is the number of  5c coins. It is evident that  a<=1,  b<=4,  c<=9,  d<=19

M:=[]:

for a from 0 to 1 do

for b from 0 to 4 do

for c from 0 to 9 do

for d from 0 to 19 do

if 50*a+20*b+10*c+5*d>100 and convert([seq(seq(seq(seq(50*i+20*j+10*k+5*l<>100, i=0..a),j=0..b), k=0..c), l=0..d)], `and`) then M:=[op(M),[a,b,c,d]] fi:

od: od: od: od:

M;

max(seq(50*M[i,1]+20*M[i,2]+10*M[i,3]+5*M[i,4], i=1..nops(M)));

Another way with two arguments:

seq(a[4 - i], i=1 .. 3);

Carl's  method does not work for me (Maple 16).

The function  x->x^x  is defined only for  x>0 .

Your example A can be easily solved by hand, if we note that this limit is the derivative of the function  f  at  the point  x=n :

f:=x->x^x:

Limit(((n+k)^(n+k)-n^n)/k, k=0)=D(f)(n);

Direct calculation of the limit with Maple give the same result:

limit(((n+k)^(n+k)-n^n)/k, k=0);

If you need in manual solution or want to verify Carl's result you can rewrite the sum to reduce it to simple definite integral:

sum(1/(n*(1+k/n)), k = 1 .. n):

limit(%, n = infinity) = Int(1/(n+1), n = 0 .. 1);

int(1/(n+1), n = 0 .. 1);

Markiyan told you about the decision in the complex domain.
In the real domain, this equation has only two solutions  {a=0, b=0}  and  {a=1, b=1} . This is evident from the plot of your system:

sys := [b-a*sqrt(1+a^2+b^2)-a^2*(a*b-sqrt(1+a^2+b^2)) = 0, a-b*sqrt(1+a^2+b^2)-b^2*(a*b-sqrt(1+a^2+b^2)) = 0]:

plots[implicitplot](sys, a = -5 .. 5, b = -5 .. 5, color = [red, blue], numpoints = 50000);

For a rigorous solution try the command  RealDomain[solve]  or solve by hand. In the manual solution it is useful  firstly to factor the equations:

sys := [b-a*sqrt(1+a^2+b^2)-a^2*(a*b-sqrt(1+a^2+b^2)) = 0, a-b*sqrt(1+a^2+b^2)-b^2*(a*b-sqrt(1+a^2+b^2)) = 0]:

factor(sys);

Look this application  http://www.maplesoft.com/applications/view.aspx?SID=146470

restart;

for i do

a:=ithprime(i):

if a>300 then break fi:

L[i]:=a:

end do:

L:=[seq(L[k], k=1..i-1)];

L := [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293]

Output is  more compact and can be used in the future, for example:

nops(L),  L[20];
                                   62, 71

Thanks Markiyan for your latest comment! This  hinted to me the idea of ​​solution.                                                    

First, we add the two equations. The resulting equation is a consequence of the original system and contains all of its roots. Then we estimate the left and right sides:

sys := [sqrt(sin(x)^2+1/sin(x)^2)+sqrt(cos(y)^2+1/cos(y)^2) = sqrt(20*y/(x+y)), sqrt(sin(y)^2+1/sin(y)^2)+sqrt(cos(x)^2+1/cos(x)^2) = sqrt(20*x/(x+y))]:

A := lhs(sys[1])+lhs(sys[2]);   B := rhs(sys[1])+rhs(sys[2]);

M := minimize(2*(eval(op(1, A)+op(2, A), [sin(x)^2 = p, 1/sin(x)^2 = 1/p, cos(x)^2 = 1-p, 1/cos(x)^2 = 1/(1-p)])), p = 0 .. 1, location);

N := maximize(eval(B, y = u*x), u = 0 .. infinity, location);

We see that the minimum value of the left-hand side is equal to the maximum value of the right side. The maximum value of the right-hand side is attained for any  x=y  (of cause  x<>0 and y<>0). The minimum value of the left-hand side is attained for any  x  and  y  such that  sin(x)^2 = 1/2  and  cos(x)^2 = 1/2 .

We find:

allvalues(solve({x = y, cos(x)^2 = 1/2, sin(x)^2 = 1/2}, AllSolutions));

The resulting solutions can be written more compactly:

x = Pi/4 + Pi*n/2,  y = Pi/4 + Pi*n/2,  n is integer

Since  the equation  A=B  is  a consequence of the original system, the solutions should be checked:

simplify(eval(sys, {x = (1/4)*Pi+(1/2)*n*Pi, y = (1/4)*Pi+(1/2)*n*Pi, 1/cos(x)^2 = 2, cos(x)^2 = 1/2, 1/cos(y)^2 = 2, cos(y)^2 = 1/2, 1/sin(x)^2 = 2, sin(x)^2 = 1/2, 1/sin(y)^2 = 2, sin(y)^2 = 1/2})) assuming n::integer;

                      [sqrt(5)*sqrt(2) = sqrt(5)*sqrt(2), sqrt(5)*sqrt(2) = sqrt(5)*sqrt(2)]

We prove that if the number  n  contains a lot of numbers (more than 6) it will not double.

At first, the idea of ​​proof by example. Let  n  has 7 digits. Then inequality  1000000 <= n <= 9999999  holds. We have 2000000<=2* n<=19999998 . Then  2*n[in base 7] <= 7^7+add(6*7^k, k=0..6) = 1647085 < 2000000

The general case reduces to the proof the inequality 

7^m+sum(6*7^k, k=0..m-1)<2*10^(m-1);   # m>=7

                        2*7^m - 1 < 2*10^(m-1)

Is sufficient to prove the stronger inequality  2*7^m < 2*10^(m-1)  for  m>=7  . It is equivalent to  10<(10/7)^m

isolve(10<(10/7)^m);
about(_NN1);

Checking for  n<1000000

t:=time():

N:=0:

for a0 from 0 to 1 do

for a1 from 0 to 6 do

for a2 from 0 to 6 do

for a3 from 0 to 6 do

for a4 from 0 to 6 do

for a5 from 0 to 6 do

for a6 from 0 by 2 to 6 do

a:=add(a||i*7^(6-i), i=0..6): b:=add(a||i*10^(6-i), i=0..6):

if 2*a=b then N:=N+1: L[N]:=a fi:

od: od: od: od: od: od: od:

[seq(L[i], i=1..N)];

time()-t;

                               [0, 51, 102, 105, 153, 156, 207, 210, 258, 261, 312, 315]
                                                                       2.188

Carl. your code is compact and elegant, but it works too slowly. Can you explain why?

Instead of  f(2)  write  eval(f(2), a = 2)  or change the procedure:

f := proc (b)

if b <= 0 then 0

elif b <= evalf(e*sin((1/2)*B)) then eval(eq1, a = b)

elif b <= evalf(2*e*sin((1/2)*B)) then eval(eq2, a = b)

else 0 end if

end proc;

Example:

f(2);

      1.358226754

restart;

g := proc(i)

if i = 1 then a else 0 fi

end proc:

h :=f->sum('g(i)', i=1 .. f):

h(3);

                 a

The use of the package  DirectSearch  is not proof in the mathematical sense.

The problem reduces to the optimization of a function of two variables because we can assume that  a=x, b=y, c=1  with restrictions

x^2<=y^2+1, y^2<=x^2+1, 1<=x^2+y^2

I took nonstrict inequalities, as maximum and minimum may be achieved within the domain or on its boundary.

Plot of the domain

plots[implicitplot]([x^2+y^2=1, y^2-x^2=1, x^2-y^2=1], x=0..4, y=0..4, color=black, thickness=2);

The domain of the function  is restricted by 3 smooth lines.

Find the function:

Expr:=proc(a,b,c)

local p, S, R;

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

S:=sqrt(p*(p-a)*(p-b)*(p-c));

R:=a*b*c/4/S;

simplify(R*p/(2*a*R+b*c));

end proc:

f:=unapply(Expr(x,y,1), x,y); 

Find the critical points of  f  within the domain:

solve({diff(f(x,y),x),diff(f(x,y),y)});

simplify(eval(f(x,y), %)); 

Find maximum and minimum of  f  on the boundaries of the domain: 

A:=subs([x=cos(t), y=sin(t)], f(x,y)):

simplify([maximize(A, t=0..Pi/2, location)]); evalf(%);

simplify([minimize(A, t=0..Pi/2, location)]); evalf(%);

B:=subs([x=cosh(t), y=sinh(t)], f(x,y)):

simplify([maximize(B, t=0..infinity,location)]); evalf(%);

simplify([minimize(B, t=0..infinity, location)]); evalf(%);

C:=subs([y=cosh(t), x=sinh(t)], f(x,y)):

simplify([maximize(C, t=0..infinity,location)]); evalf(%);

simplify([minimize(C, t=0..infinity, location)]); evalf(%);

The inequality  2/5 <= R*p/(2*a*R+b*c) < 1/2  is prooved. The lower limit  2/5  is reached for the triangle with sides  6/5*C, C, C, where C is arbitrary positive. The upper limit  1/2  is not achieved for any acute triangle.