I was wrong about "a lot of time." Still, the capabilities of modern computers boggle my mind ! Less than 7 minutes the program went through  21^6 = 85766121 variants and  found all solutions. The first 10 and last 10 solutions displayed.

restart;

ts := time(): N := 0:

for c from -10 to 10 do
for d from -10 to 10 do
for e from -10 to 10 do
for f from -10 to 10 do
for g from -10 to 10 do
for h from -10 to 10 do
a := c+2*d+4*e+6*f+10*g+22*h-16: b := -2*c-3*d-5*e-7*f-11*g-23*h+20:
if a >= -10 and a <= 10 and b >= -10 and b <= 10 then N := N+1:
Sol[N] := [a, b, c, d, e, f, g, h] end if:
end do: end do: end do: end do: end do: end do:

N; seq(Sol[i], i = 1 .. 10); seq(Sol[i], i = N-9 .. N);
time()-ts;

Of a given system  a  and  b  can be easily expressed  through the other variables, and then in the loop we find one solution:

restart;

for c from -10 to 10 do

for d from -10 to 10 do

for e from -10 to 10 do

for f from -10 to 10 do

for g from -10 to 10 do

for h from -10 to 10 do

a:=c+2*d+4*e+6*f+10*g+22*h-16: b:=-2*c-3*d-5*e-7*f-11*g-23*h+20:

if a>=-10 and a<=10 and b>=-10 and b<=10 then Sol:=[a,b,c,d,e,f,g,h]: break: fi:

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

Sol;

                                        [-10, -10, 10, 10, 10, 1, -7, 0]

Finding all solutions will require a lot of time.

Use  mtaylor  command. In your example,  a  can not be 0, because it is in the denominator:

mtaylor((1+1/(2*a^2)+(1/2)*b^2-1/(2*a))^a, [a=1,b=0], 5);

From the plot  

plot([ln(x+1),x^2-1], x=-2..3, -3..3);  

we can see that there are exactly two points of intersection of the graphs of functions. Solving numerically the corresponding equations, we finally obtain

RealRange(Open(-1),fsolve(ln(x+1)=x^2-1, x=-infinity..0)),RealRange(fsolve(ln(x+1)=x^2-1, x=0..infinity), infinity) ;

This can be done as follows (example):

A := plot(sqrt(x), x = 0 .. 4, tickmarks = [0, default]):
B := seq(plottools[line]([k, 0], [k, 0.4e-1]), k = 1 .. 4):
C := seq(plottools[line]([k, 0], [k, 0.2e-1]), k = 0 .. 4, 0.1):
T := plots[textplot]([seq([k, 0.3e-1, k], k = 0.5 .. 4, 0.5)], align = above): plots[display](A, B, C, T, scaling = constrained);

It is easy to derive a formula for the  n-th  member in this sequence:

f:=n->ceil((-1+sqrt(1+8*n))/2);

Examples:

f(2013);   f(10^10);

                 63
             141421

PS.  About  ceil  function see  ?trunc  and  http://en.wikipedia.org/wiki/Floor_and_ceiling_functions

Basic concepts, motivation, examples and references, see here

http://en.wikipedia.org/wiki/Group_representation

Use  curve=a*x  option.

Example:

with(CurveFitting):

LeastSquares([1, 2, 3], [1, 5/2, 3], x);  #  without the option

LeastSquares([1, 2, 3], [1, 5/2, 3], x, curve=a*x);  #  with the option

                                         1/6+x

                                       15/14*x

If the triangle has integer sides and the coordinates of its vertices are integers, it is obvious that the length of each side (if it is not parallel to the axes) has the hypotenuse in a Pythagorean triangle. But not exists  Pythagorean  triangle with hypotenuse to be equal   28:

L:=[]:

for b from 2 to 39 do

for a from 1 to b-1 do

for c from 3 to 40 do

if a^2+b^2=c^2 then L:=[op(L),[a,b,c]] fi:

od: od: od:

L;  #  All Pythagorean triangles with the hypotenuse  <=40

[[3, 4, 5], [6, 8, 10], [5, 12, 13], [9, 12, 15], [8, 15, 17], [12, 16, 20], [15, 20, 25], [20, 21, 29], [7, 24, 25], [10, 24, 26], [18, 24, 30], [21, 28, 35], [16, 30, 34], [24, 32, 40], [12, 35, 37], [15, 36, 39]]

Therefore, we can consider  A(0, 0) ,  B(28, 0),  C(x,y) :

solve({x^2+y^2=17^2, (28-x)^2+y^2=39^2, y>0});

                                    {x = -8, y = 15}

PS: Using Pythagorean triangles, you can significantly increase the search range, e.g., to  x=-1000..1000, y=-1000..1000 

Of course, your example may be easier to solve by hand, but if  with Maple, you can use the usual substitution in  Green's formula (from Wiki)

eval((diff(M(x,y),x)-diff(L(x,y),y)), {L(x,y)=diff(f(x,y),y), M(x,y)=-diff(f(x,y),x)});

simplify(%,{diff(f(x,y),x,x)+diff(f(x,y),y,y)=0});

int(%, [y=y1(x)..y2(x), x=a..b]); 

knn:=Matrix(10,1):

for m from 1 to 10 do

P:=fsolve({y=1/kn,y=-tan(kn)},{kn=(m-1)*Pi..m*Pi,y=0..5});

knn[m,1]:=rhs(P[1]);

end do:

knn[3,1];

knn;

The solution closest to brian's decision:

P := {3, 4, 5, 6, 7, 8, 28, 30, 35}:
Total:=simplify(convert(P, `*`)^(1/3));
A := combinat:-choose(P, 3):
B:=select(a->is(convert(a,`*`)=Total), A);
C:=combinat:-choose(B, 3):
op(select(c->is(`union`(op(c))=P), C));

There is the unique solution:

S:={3,4,5,6,7,8,28,30,35}: L:={}:

T:=combinat[choose](S, 3):

for t in T do

S1:=S minus t:

U:=combinat[choose](S1, 3):

for u in U do

S2:=S1 minus u:

if convert(t, `*`)=convert(u, `*`) and convert(u, `*`)=convert(S2, `*`) then

L:={op(L), {t,u,S2}} fi:

od: od:

L;

                                  {{3, 8, 35}, {4, 7, 30}, {5, 6, 28}}

If you do not use any special properties of the dodecahedron, then the solution is not very complex, but rather cumbersome. 

We will look for the vertex coordinates of a dodecahedron inscribed in the unit sphere. We assume that the uppermost and lowermost vertex lies on the axis Oz, the vertex adjacent to the top lies in the plane xOz. First we look for the edge length of the dodecahedron (see figure):

restart;
AB:=2*a*convert(sin(54*Pi/180), radical);
BK:=AB/sqrt(3):
OK:=sqrt(1-BK^2):
EK:=sqrt(a^2-BK^2):
solve(OK+EK=1);
a:=op(select(t->is(t>0), [%]));  # the length of the edge

Then we find the coordinates of the vertex  A  and vertex  B , and  C  using the rotation around the axis Oz:

c:=a/2:  # Cosine of the angle OEA
with(geom3d):
point(E, 0, 0, 1):
point(A,simplify(a*sqrt(1-c^2)), 0, simplify(1-a*c)):
coordinates(A);
rotation(B, A, 2*Pi/3, line(l,[0,0,t],t)):
rotation(C, A, 4*Pi/3, line(l,[0,0,t],t)):
coordinates(B);
coordinates(C);

Next, we find the coordinates of the vertices  F and G, belonging to the face  AFGBE, using a rotation about an axis passing through the origin and center of this face:

EO1:=simplify(a/2/convert(cos(54*Pi/180), radical)):
midpoint(L, A, B):
point(O1, convert(simplify(convert(coordinates(E), Vector)+convert(coordinates(L)-coordinates(E), Vector)*EO1/distance(E,L)), list)):
dsegment(OO1, [point(o, 0,0,0), O1]):
rotation(F, A, 72*Pi/180, line(l1,[o,OO1])):
rationalize(simplify(convert(coordinates(F), radical)));
point(F, %):
rotation(G, A, 144*Pi/180, line(l1,[o,OO1])):
expand(rationalize(simplify(convert(coordinates(G), radical))));
point(G, %):
Pol1:=plottools[polygon]([coordinates(A),coordinates(F),coordinates(G),coordinates(B),coordinates(E)], style=line, thickness=2, color=red):
s1:=plottools[line](coordinates(A),coordinates(B), linestyle=3, color=black):
s2:=plottools[line](coordinates(E),coordinates(O1), linestyle=3, color=black):
s3:=plottools[line](coordinates(o),coordinates(O1), linestyle=3, color=black):
T1:=plots[textplot3d]([[0,0,0,"O"], [op(coordinates(A)),"A"], [op(coordinates(B)),"B"], [op(coordinates(E)),"E"], [op(coordinates(F)),"F"], [op(coordinates(G)),"G"], [op(coordinates(L)),"L"], [op(coordinates(O1)),"O1"]], align=left, color=black, font=[TIMES,ROMAN,16]):
plots[display](Pol1, s1, s2, s3, T1, axes=normal, scaling=constrained, tickmarks=[0, 0, 0], view=[-1.4..1.4,-1.4..1.4,-0.4..1.4]);

Then we find the coordinates of all vertices yet two faces using  rotation of the face  AFGBE  around the axis Oz:

rotation(M, F, 2*Pi/3, line(l,[0,0,t],t)):
rotation(N, G, 2*Pi/3, line(l,[0,0,t],t)):
rotation(P, F, 4*Pi/3, line(l,[0,0,t],t)):
rotation(Q, G, 4*Pi/3, line(l,[0,0,t],t)):
coordinates(M);
expand(rationalize(coordinates(N)));
coordinates(P);
expand(rationalize(coordinates(Q)));
Pol2:=plottools[polygon]([coordinates(B),coordinates(M),coordinates(N),coordinates(C),coordinates(E)], style=line, thickness=2, color=red):
Pol3:=plottools[polygon]([coordinates(C),coordinates(P),coordinates(Q),coordinates(A),coordinates(E)], style=line, thickness=2, color=red):
T2:=plots[textplot3d]([[0,0,0,"O"], [op(coordinates(A)),"A"], [op(coordinates(B)),"B"], [op(coordinates(C)),"C"], [op(coordinates(E)),"E"], [op(coordinates(F)),"F"], [op(coordinates(G)),"G"], [op(coordinates(M)),"M"], [op(coordinates(N)),"N"], [op(coordinates(P)),"P"], [op(coordinates(Q)),"Q"]], align=left, color=black, font=[TIMES,ROMAN,16]):
plots[display](Pol1, Pol2, Pol3, T2, axes=normal, scaling=constrained, tickmarks=[0, 0, 0], view=[-1.4..1.4,-1.4..1.4,-0.4..1.4]);

Finally, we find all the other vertices of the dodecahedron, using the reflection about the plane xOy and rotation around Oz:

X:=plots[display](Pol1, Pol2, Pol3):
Y:=plottools[rotate](plottools[reflect](X, [[0,0,0], [1,0,0], [0,1,0]]), Pi/3, [[0,0,0], [0,0,1]]):
rotation(P1, reflection(P0, P, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(Q1, reflection(Q0, Q, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(F1, reflection(F0, F, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(G1, reflection(G0, G, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(M1, reflection(M0, M, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(N1, reflection(N0, N, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(A1, reflection(A0, A, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(B1, reflection(B0, B, plane(xOy, z=0, [x,y,z])), Pi/3, l):
rotation(C1, reflection(C0, C, plane(xOy, z=0, [x,y,z])), Pi/3, l):
reflection(E1, E, plane(xOy, z=0, [x,y,z])):
T4:=plots[textplot3d]([[op(coordinates(A1)),"A1"], [op(coordinates(B1)),"B1"], [op(coordinates(C1)),"C1"], [op(coordinates(E1)),"E1"], [op(coordinates(F1)),"F1"], [op(coordinates(G1)),"G1"], [op(coordinates(M1)),"M1"], [op(coordinates(N1)),"N1"], [op(coordinates(P1)),"P1"], [op(coordinates(Q1)),"Q1"]], align=left, color=black, font=[TIMES,ROMAN,14]):
P1Q:=plottools[line](coordinates(P1),coordinates(Q), color=red, thickness=2):
Q1F:=plottools[line](coordinates(Q1),coordinates(F), color=red, thickness=2):
F1G:=plottools[line](coordinates(F1),coordinates(G), color=red, thickness=2):
G1M:=plottools[line](coordinates(G1),coordinates(M), color=red, thickness=2):
M1N:=plottools[line](coordinates(M1),coordinates(N), color=red, thickness=2):
N1P:=plottools[line](coordinates(N1),coordinates(P), color=red, thickness=2):
T5:=plots[textplot3d]([[1.7,0,0,"x"], [0,1.7,0,"y"], [0,0,1.7,"z"]], align=above, color=black, font=[TIMES,ROMAN,16]):
plots[display](X, Y, T2, T5, T4, P1Q, Q1F, F1G, G1M, M1N, N1P, axes=normal, scaling=constrained, tickmarks=[0, 0, 0], view=[-1.7..1.7,-1.7..1.7,-1.7..1.7]);

 dodecahedron.mws

aa:=``(delta[1]*n[1])/``(n[f1]*delta[1]);

 If you want to reduce the fraction, then

expand(aa);