The symmetry group of the space that preserves the origin and integer cubic lattice consists of 48 members: 2^3=8 changes of signs of  coordinates ( reflections with respect to the coordinate planes ) and 6 permutations of the coordinates ( reflections with respect to the bisector planes). Therefore we get 8 * 6 = 48. 

The following code replicates obtained above solutions (the list  L ). Then we find the positions of the original examples   A(-13,-5,5), B(-5,11,-11), C(-3,-9,15)  and  A(-6,6,-1), B(5,-1,-3), C(2,10,7)  in this list :

p:=combinat[permute](3);  # Group of the permutations of 3 elements

q:=combinat[permute]([-1$3, 1$3],3);  # Group of the changes the signs of coordinates

L1:=map(t->seq([seq([seq(t[i][k[j]], j=1..3)], i=1..3)], k=p), L):

L2:=map(t->seq([seq([seq(t[i][j]*k[j], j=1..3)], i=1..3)], k=q), L1):  # The final list of the all solutions (nops(L2) =3567*48=171216 solutions)

P1:=[[-13,-5,5], [-5,11,-11], [-3,-9,15]]:  P2:=[[-6,6,-1],[-5,-1,-3],[2,10,7]]:

L3:=map(t->[t[1]+P1[1], t[2]+P1[1], t[3]+P1[1]], L2):

L4:=map(t->[t[1]+P2[1], t[2]+P2[1], t[3]+P2[1]], L2):

N1:=ListTools[Search](P1, L3);  N2:=ListTools[Search](P2, L4);

Next, check these two triangles:

L3[N1];  L4[N2];  Q1:=map(convert,P1,Vector);  Q2:=map(convert,P2,Vector);

СenterOfCircle(op(Q1));

Orthocenter(op(Q1));

Centroid(op(Q1));

СenterOfCircle(op(Q2));

Orthocenter(op(Q2));

Centroid(op(Q2));

If you want all the integrals are calculated immediately, write

restart;

phi[1]:=(x,y)->exp(x-y);

for jj to 2 do

phi[jj+1]:=unapply(int(phi[1](x,s)*phi[jj](s,y),s=y..x), x,y);

end do;

If you want to calculate integrals later, replace  int  by  Int .

For purely visual presentation (without assignments) you can write

restart;

phi[1]=((x,y)->exp(x-y));

for jj to 2 do

phi[jj+1]=unapply(Int(phi[1](x,s)*phi[jj](s,y),s=y..x), x,y);

end do;

The first example for standard worksheet:

%piecewise(``, x+y-2*z = 1, ``, 2*x+y-3*z = 5, ``, -x+y+z = 1);

 There is also  the procedure from Aladjev's library that works in both classic and standard (I checked in Maple 12). This procedure  additionally numbers equations of the system.

braces := proc(L::{list(relation), set(relation)})

local a, f, k;

global __vsvak38;

assign(a = "__vsv38ak:=piecewise(", f = cat(currentdir(), "/__$$$__"));

seq(assign('a'=cat(a, convert(L[k], string), ",", convert([k], string), ",")), k=1..nops(L));

writeline(f, cat(a[1 .. -2], ");")), close(f);

(proc() read f end proc)(f), fremove(f), __vsv38ak, unassign('__vsv38ak')

end proc: 

Example of use:

braces([x+y-2*z=1, 2*x+y-3*z=5, -x+y+z=1]);

 The procedure  braces  is useful if you want to refer to the individual equations of the system.

If the first version of the procedure  brace  has problems, then bring another variant which seems to work everywhere:

restart;

braces := proc (L::anything)

local a, f, k, t;

if whattype(L) = function then try assign(a = convert(L, 'list')),

assign(f = [seq(a[k], k = {seq(2*t-1, t = 1 .. (1/2)*nops(a))})]),

`if`(nops(f) = 0, NULL, f) catch "": NULL end try;

elif type(L, {list(relation), set(relation)}) then assign(a = "piecewise("); seq(assign('a' = cat(a, convert(L[k], 'string'), ",", convert([k], 'string'), ",")), k = 1 .. nops(L));

eval(parse(cat(a[1 .. -2], ")")));

else  end if;

end proc:

Rewrite this equation as follows:

2^(sin(x)^4+sin(x)^2-1)+sin(x)^2 = 2^(cos(x)^4+cos(x)^2-1)+cos(x)^2

and consider the function  t->2^(t^2+t-1)+t . This function is strictly increasing for  t>=0  and so the original equation is equivalent to the equation  sin(x)^2=cos(x)^2

restart;

f := 2^(t^2+t-1)+t:

solve({diff(f, t) > 0, t >= 0});  # The proof of increasing the function  f  for  t>=0

solve(combine(sin(x)^2 = cos(x)^2), AllSolutions);  # The all solutions

z2 := cos(x)^(n-1)*sin(x)+n*(Int(cos(x)^(n-2), x))-(Int(cos(x)^(n-2), x));

collect(z2, Int);

Just apply  simplify  command to the last line:

JQDFE := simplify(eval(J, [S = Lambda/(beta-mu[A]), T = 0, H = Lambda*(beta/(mu+mu[A])-1)/(beta-mu[A]), C = 0, C1 = 0, C2 = 0, C1M = 0, C2M = 0]));

The following code finds all the such triangles in specified ranges. Without loss of generality, to reduce the volume of calculations, we take one vertex at the origin, and take one vertex in the first octant with non-decreasing coordinates. It is obvious that by translations and reflections with respect to the coordinate planes and relatively bisector planes, any triangle can be reduced to this form.

restart;

CP:=proc(A::Vector, B::Vector, C::Vector)

uses LinearAlgebra;

convert(CrossProduct(B-A, C-A), list);

end proc:

СenterOfCircle:=proc(A::Vector, B::Vector, C::Vector)

local P1, P2, P3, X, x, y, z;

uses LinearAlgebra;

X:=<x,y,z>;

P1:=Determinant(Matrix([X-A, B-A, C-A]))=0;

P2:=DotProduct((X-(A+B)/2),(B-A),conjugate=false)=0;

P3:=DotProduct((X-(A+C)/2),(C-A),conjugate=false)=0;

solve({P1,P2,P3}, {x,y,z});

assign(%);

[x,y,z];

end proc:

Orthocenter:=proc(A::Vector, B::Vector, C::Vector)

local P1, P2, P3, X, x, y, z;

uses LinearAlgebra;

X:=<x,y,z>;

P1:=Determinant(Matrix([X-A, B-A, C-A]))=0;

P2:=DotProduct((X-A),(B-C),conjugate=false)=0;

P3:=DotProduct((X-B),(C-A),conjugate=false)=0;

solve({P1,P2,P3}, {x,y,z});

assign(%);

[x,y,z];

end proc:

Centroid:=proc(A::Vector, B::Vector, C::Vector)

uses LinearAlgebra;

convert((A+B+C)/3, list);

end proc:

A, B, C:=<0,0,0>, <x1,y1,z1>, <x2,y2,z2>:

T:=CP(A,B,C):

CC:=СenterOfCircle(A,B,C):

Oc:=Orthocenter(A,B,C):

Ct:=Centroid(A,B,C):

k:=0:

for x1 from 0 to 20 do

for y1 from x1 to 20 do

for z1 from y1 to 20 do

for x2 from -10 to 10 do

for y2 from -10 to 10 do

for z2 from -10 to 10 do

A, B, C:=<0,0,0>, <x1,y1,z1>, <x2,y2,z2>:

if (T[1]<>0 or T[2]<>0 or T[3]<>0) then if (CC[1]::integer and CC[2]::integer and CC[3]::integer) and (Oc[1]::integer and Oc[2]::integer and Oc[3]::integer) and (Ct[1]::integer and Ct[2]::integer and Ct[3]::integer) then k:=k+1: L[k]:=[A,B,C]: fi: fi:

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

L:=convert(L,list):

nops(L); # Number of solutions

seq(L[1+350*(i-1)], i=1..10); # 10 solutions from the total 3567 ones

This code searches 3567 solutions for about 8 minutes. 

Because  abs(x+y)=sqrt((x+y)^2)  then

restart;

Eqs:= [sqrt((x-1)^2+(y-5)^2)+(1/2)*sqrt((x+y)^2) = 3*sqrt(2), sqrt(abs(x+2)) = (2-y)]:

solve(Eqs);

                                                       {x = -2, y = 2}

Such examples are useful to check by the visualization:

plots[implicitplot]([sqrt((x-1)^2+(y-5)^2)+(1/2)*abs(x+y) = 3*sqrt(2), sqrt(abs(x+2)) = 2-y], x = -7 .. 3, y = -3 .. 6, color = [red, blue], thickness = 2, numpoints = 10000, scaling = constrained);

1. Your system  for solving by  fsolve  is composed wrong. First, you define a set of equations  eq[i]  based on  ode1 , then under the same names  other equations based on  ode2.

2. To solve a system by  fsolve  the number of equations must equal the number of unknowns.

May be as follows:

solve({x=z,y=z, z=t}, {x,y,z});  # t is the parameter

             {x = t, y = t, z = t}

If not, please give the full text of the problem.

To be more precise, the planets revolve not in circles, but in ellipses, and the Sun is at one focus. So first we take the polar equation of the ellipse (the Sun is in the pole), and then all the constructions using parametric equations, based on the original polar equations. 

For example, all parameters are arbitrary and have no relation to the actual planets. The plane of rotation of the red planet taken at an angle of 18 degrees with respect to the plane of the blue planet. Rotation period of the blue planet taken 2 times more than the red planet. Initial phases are also different.

As a result of the code, you can see two animations. The first animation shows the motion of the two planets in their orbits, the second one also shows the planes of the orbits:

restart;
r1 := (3*(1-0.4^2))/(1-0.4*cos(t)):
r2 := (2*(1-0.5^2))/(1-0.5*cos(2*t-(1/4)*Pi)):
R := <cos(phi), 0, -sin(phi); 0, 1, 0; sin(phi), 0, cos(phi)>:
phi := (1/10)*Pi:
L1 := [r1*cos(t), r1*sin(t), 0]:
L2 := convert(R.`<,>`(r2*cos(2*t-(1/4)*Pi), r2*sin(2*t-(1/4)*Pi), 0), list): T1 := plots[spacecurve](L1, t = 0 .. 2*Pi, linestyle = 3, color = blue, thickness = 2):
T2 := plots[spacecurve](L2, t = 0 .. 2*Pi, linestyle = 3, color = red, thickness = 2):
Sun := plottools[sphere]([0, 0, 0], 0.2, style = surface, color = yellow):
S1 := plot3d([u*cos(t), u*sin(t), 0], u = 0 .. r1, t = 0 .. 2*Pi, style = surface, color = "LightBlue", transparency = 0.5):
S2 := plot3d(convert(R.<v*cos(2*t-(1/4)*Pi), v*sin(2*t-(1/4)*Pi), 0>, list), v = 0 .. r2, t = 0 .. 2*Pi, style = surface, color = pink):
P1 := proc (x, y) plots[display](plottools[sphere](x, y, style = surface, color = blue)) end proc:
P2 := proc (x, y) plots[display](plottools[sphere](x, y, style = surface, color = red)) end proc:
Planet1 := plots[animate](P1, [L1, 0.15], t = 0 .. 2*Pi, frames = 50):
Planet2 := plots[animate](P2, [L2, 0.1], t = 0 .. 2*Pi, frames = 50):

plots[display](T1, T2, Sun, Planet1, Planet2, scaling = constrained, axes = framed, orientation = [65, 75], lightmodel = light1);
plots[display](T1, T2, S1, S2, Sun, Planet1, Planet2, scaling = constrained, axes = framed, orientation = [65, 75], lightmodel = light1);

evalf(map(log10, [25,5,1,10,4,20]));

map(t->10^t, %);

[1.397940008, .6989700041, 0., 1., .6020599914, 1.301029996]

[24.99999996, 4.999999997, 1., 10., 4.000000001, 20.00000002]

Maybe you mean the plane domain  sqrt(x)<y<=1/sqrt(2) ?

Plotting in Maple 12 classic:

restart;

A:=plot(sqrt(x), x=-0.2..0.7, color=blue, thickness=2):

B:=plot(1/sqrt(2), x=-0.2..0.7, color=red, thickness=2):

A1:=plot(sqrt(x), x=0..1/2, color=white, filled=true):

B1:=plot(1/sqrt(2), x=0..1/2, color=yellow, filled=true):

B2:=plot(1/sqrt(2), x=0..1/2, color=yellow, thickness=3):

C:=plot([1/2, t, t=0..1/sqrt(2)], linestyle=3, color=black):

C1:=plots[textplot]([1/2,-0.02,1/2], font=[TIMES,ROMAN,16]):

plots[display](B2, A, B, A1, B1, C, C1, scaling=constrained);

Of course, manual axes are very laborious for a single application. But it is easy to write a procedure in which to provide any options. 

restart;

export_plot_options := axes = none, size = [850, 850]:
points := [seq([seq(exp(-(x^2+y^2)*(1/100)), x = -10 .. 10.0)], y = -10 .. 10.0)]:

A := plots:-listdensityplot(points, export_plot_options):
B := plots[arrow]([.43-2, .43], [27, 0], width = 0.4e-1, head_length = [0.4e-1, relative], head_width = [15, relative]):
C := plots[arrow]([.43, .43-2], [0, 27], width = 0.4e-1, head_length = [0.4e-1, relative], head_width = [15, relative]):
Ticks := seq(plottools[line]([2*i, .43], [2*i, .43-.3], thickness = 3), i = 0 .. 11), seq(plottools[line]([.43-.3, 2*i], [.43, 2*i], thickness = 3), i = 0 .. 11):
SubTicks := seq(plottools[line]([i, .43], [i, .43-.1], thickness = 3), i = -1 .. 23), seq(plottools[line]([.43-.1, i], [.43, i], thickness = 3), i = -1 .. 23):
Text := plots[textplot]([seq([2*i, .43-.8, 2*i], i = 1 .. 11), [25, .43-.8, x], seq([.43-1, 2*i, 2*i], i = 1 .. 11), [.43-1, 25, y]], font = [TIMES, roman, 30]):

plots[display](A, B, C, Ticks, SubTicks, Text);

I fixed a few mistakes in your code.. Symbolic solution is very cumbersome, so I show only the approximate one and its plot:

restart;

eq:=a11*diff(phi(x),x,x,x,x)+a22*diff(phi(x),x,x)+a33*phi(x):

inc:=phi(a)=sigma1, phi(-a)=sigma1, D(phi)(a)=0, D(phi)(-a)=0:

a11:=2.731e-10:

a22:=-1.651e-9:

a33:=3.09027e-10:

a:=35.714:

sigma1:=200e6:

dsolve({eq, inc});

evalf(%);

assign(%);

plot(phi(x), x=-a-1..a+1, thickness=2);