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.

Tuples:=proc(n, b)

local L, It;

L:=[seq([k], k=0..b)];

if n=1 then return L fi;

It:=proc(M)

[seq(seq([k, op(M[i])], k=0..b), i=1..nops(M))];

end proc;

(It@@(n-1))(L);

end proc;

Example:

Tuples(4, 2);

[[0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 0, 0], [0, 1, 0, 0], [1, 1, 0, 0], [2, 1, 0, 0], [0, 2, 0, 0], [1, 2, 0, 0], [2, 2, 0, 0], [0, 0, 1, 0], [1, 0, 1, 0], [2, 0, 1, 0], [0, 1, 1, 0], [1, 1, 1, 0], [2, 1, 1, 0], [0, 2, 1, 0], [1, 2, 1, 0], [2, 2, 1, 0], [0, 0, 2, 0], [1, 0, 2, 0], [2, 0, 2, 0], [0, 1, 2, 0], [1, 1, 2, 0], [2, 1, 2, 0], [0, 2, 2, 0], [1, 2, 2, 0], [2, 2, 2, 0], [0, 0, 0, 1], [1, 0, 0, 1], [2, 0, 0, 1], [0, 1, 0, 1], [1, 1, 0, 1], [2, 1, 0, 1], [0, 2, 0, 1], [1, 2, 0, 1], [2, 2, 0, 1], [0, 0, 1, 1], [1, 0, 1, 1], [2, 0, 1, 1], [0, 1, 1, 1], [1, 1, 1, 1], [2, 1, 1, 1], [0, 2, 1, 1], [1, 2, 1, 1], [2, 2, 1, 1], [0, 0, 2, 1], [1, 0, 2, 1], [2, 0, 2, 1], [0, 1, 2, 1], [1, 1, 2, 1], [2, 1, 2, 1], [0, 2, 2, 1], [1, 2, 2, 1], [2, 2, 2, 1], [0, 0, 0, 2], [1, 0, 0, 2], [2, 0, 0, 2], [0, 1, 0, 2], [1, 1, 0, 2], [2, 1, 0, 2], [0, 2, 0, 2], [1, 2, 0, 2], [2, 2, 0, 2], [0, 0, 1, 2], [1, 0, 1, 2], [2, 0, 1, 2], [0, 1, 1, 2], [1, 1, 1, 2], [2, 1, 1, 2], [0, 2, 1, 2], [1, 2, 1, 2], [2, 2, 1, 2], [0, 0, 2, 2], [1, 0, 2, 2], [2, 0, 2, 2], [0, 1, 2, 2], [1, 1, 2, 2], [2, 1, 2, 2], [0, 2, 2, 2], [1, 2, 2, 2], [2, 2, 2, 2]]

Of cause, 25^2+24^2 . You can see it:

with(plottools):

A:=curve([[0,0],[200,0],[200,200],[0,200],[0,0]], thickness=2, color=black):

B:=seq(seq(disk([4+8*i,4+8*j], 1.5, color=green), j=0..24), i=0..24):

C:=seq(seq(disk([8+8*i,8+8*j], 1.5, color=yellow), j=0..23), i=0..23):

plots[display](A, B, C, axes=none);  # All in one

plots[display](A, B, axes=none);  # Only green trees

plots[display](A, C, axes=none);  # Only yellow trees

Two squares: the first square of green trees (with a side of 25 trees), and the second one of yellow trees (with the side of 24  tree):

I usually work in the classic interface. The text of the code can without any problems be copied and pasted into a text editor of mapleprimes. But if you do so, then the front of each line of code check mark appears. So at first I copy the code into Word, and then from Word into the text editor.

If I work in the standard interface, to copy the code I first select it, and then by context menu convert to 1-D Math, then cope without any problems .

for C from 2 to 10 do

s[C] := lhs(op(allvalues(solve({K > 0, K*(K-1) > 6*C-2})))):

end do:

L := [seq(floor(s[i]+1), i = 2 .. 10)];

                                      L := [4, 5, 6, 6, 7, 7, 8, 8, 9]

If you make a change  x=sqrt(lambda) , lambda>=0  it can be clearly seen from the graphs

plot([tan(x), x], x=-Pi..10*Pi, -5..35);

that in each range  x = (1/2)*Pi+Pi*k .. (1/2)*Pi+Pi*(k+1), k>=-1 there is a single root.

Finding  the first 10 roots:

seq(fsolve(tan(x) = x, x = (1/2)*Pi+Pi*k .. (1/2)*Pi+Pi*(k+1))^2, k = -1 .. 8);

0., 20.19072856, 59.67951595, 118.8998692, 197.8578111, 296.5544121, 414.9899843, 553.1646459, 711.0784498, 888.7314224

PS. This is interesting: my list does not coincide with Carl's one.

An interesting problem!

It can be solved in different ways. In my view, the most short way is to use a double integral with the change of variables. The change  u=y/x, v=(a-x)/y  maps the original region on the rectangle.

restart;

solve({u=y/x, v=(a-x)/y}, {x,y}):

assign(%):

int(Student[MultivariateCalculus][Jacobian]([x,y], [u,v], output = determinant), [u=1/2..1, v=1..3]) assuming a>0; 

                                                            7/120*a^2

A := plots[implicitplot](max(2-r, r-5, 3*Pi*(1/4)-theta, theta-5*Pi*(1/4)) = 0, r = 0 .. 6, theta = 0 .. 2*Pi, coords = polar, axiscoordinates = polar, gridrefine = 3):

B := plottools[polygon]([[2*cos(3*Pi*(1/4)), 2*sin(3*Pi*(1/4))], [5*cos(3*Pi*(1/4)), 5*sin(3*Pi*(1/4))], seq([5*cos(3*Pi*(1/4)+(1/200)*Pi*i), 5*sin(3*Pi*(1/4)+(1/200)*Pi*i)], i = 1 .. 100), [2*cos(5*Pi*(1/4)), 2*sin(5*Pi*(1/4))], seq([2*cos(5*Pi*(1/4)-(1/100)*Pi*i), 2*sin(5*Pi*(1/4)-(1/100)*Pi*i)], i = 1 .. 49)], color = green):

plots[display](A, B);

To automate the plotting of complicated plane figures, and to calculate their areas and perimeters, you can see my work

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

Of cause, the problem can easily be solved by brute force method:

N:=0:

for n from 1 to 2013 do

if (irem(n,3)=0 and irem(n+1,4)=0) or (irem(n,4)=0 and irem(n+1,3)=0) then

N:=N+1; fi;

od:

N; 

                     336

In fact, the problem can be solved analytically for any range. The decision is based on the following arguments:

1. If the number  n  is divisible by 3, and the number  n+1  is divisible by 4, then the general formula for all such numbers is obtained as the solution of Diophantine equation  3*k+1=4*m . Similarly, if  n  is divisible by 4, and  n+1  is divisible by 3.

2. If  a .. b  is any  real range (a<=b), the number of integer points in this range is  floor(b) - ceil(a)+1 .

P:=proc(N1, N2)  # N1 and N2 specify the range N1..N2

local sol1, sol2;

isolve(3*k+1=4*m);

sol1:=solve({3*rhs(%[1])>=N1, 3*rhs(%[1])<=N2-1});

isolve(4*k+1=3*m);

sol2:=solve({4*rhs(%[1])>=N1, 4*rhs(%[1])<=N2-1});

floor(rhs(sol1[2]))-ceil(lhs(sol1[1]))+floor(rhs(sol2[2]))-ceil(lhs(sol2[1]))+2;

end proc: 

Examples:

P(1, 2014);

P(10^10, 10^20);

                      336

      16666666665000000000

If you want to use a symbol  V , and that there was no contradiction with the already defined vector  V , you can write 

V:=Vector(5):

for i to 20 do 
...
V||i =...
...
od:

Example:

V:=Vector(5);

for i to 20 do

V||i:=i^2;

od:

V||9;

There are infinitely many solutions of the form  a=i^2,  b=2*i*j,  c=j^2  (i, j  are integers),  that is the consequence of the identity

i^2*4^n+2*i*j*6^n+j^2*9^n = (i*2^n+j*3^n)^2

It remains to prove that there are no other solutions.

a := [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:

b := convert([ 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011], fraction):

c := convert([ -0.88, -8.87, -0.86, -0.82, -0.77, -0.71, -0.66, -0.62, -0.57, -0.54, -0.89, -0.88, -0.85, -0.81, -0.76, -0.71, -0.66, -0.61, -0.57, -0.53], fraction):

P1:=CurveFitting[PolynomialInterpolation](a[1..10], c[1..10], x):

P2:=CurveFitting[PolynomialInterpolation](a[11..20], c[11..20], x):

F:=unapply(expand(P1*(y-b[11])/(b[1]-b[11])+P2*(y-b[1])/(b[11]-b[1])), x, y); 

The polynomial  F  is exact on your data:

is([seq(F(a[i],b[i]), i=1..20)] = c);

                          true

Apply a function  tan  to both sides of the equation and then after of simplifications  use   isolve  command. 

The only answer   {x = 0, y = 4}

From conditions of the problem all the angles are easy to find, and from the law of sines  all sides of the triangle can be expressed by  a :

is(cos(Pi/7)^2+cos(2*Pi/7)^2+cos(4*Pi/7)^2=5/4);

is(sin(Pi/7)/sin(2*Pi/7)/a+sin(Pi/7)/sin(4*Pi/7)/a=1/a);

                                               true

                                               true

You can find out the nature of the roots of a cubic equation, without solving the equation itself, through its coefficients in terms of the discriminant and the resultant (Delta  and  Delta[0]). See  http://en.wikipedia.org/wiki/Cubic_function