I think that your system with the boundary conditions can be solved only numerically:

Digits:=20:

P:=phi(x):

Q:=psi(x):

eq1:=a11*diff(P,x,x,x,x)+a22*diff(P,x,x)+a33*P+a44*diff(Q,x,x)+a55*Q:

eq2:=a44*diff(P,x,x)+a55*P+a66*diff(Q,x,x)+a77*Q:

BV:={phi(a)=sigma1, phi(-a)=sigma1, D(phi)(a)=0, D(phi)(-a)=0, psi(a)=sigma2, psi(-a)=sigma2}:

a11:=6.36463*10^(-10):

a22:=-1.22734*10^(-9):

a33:=3.48604*10^(-10):

a44:=2.94881*10^(-11):

a55:=-5.24135*10^(-11):

a66:=-1.03829*10^(-9):

a77:=4.86344*10^(-10):

a:=1.62338:

sigma1:=1.93251*10^7:

sigma2:=9.99998*10^7:

sol:=dsolve({eq1, eq2, op(BV)}, numeric);

plots[odeplot](sol, [[x,phi(x)], [x,psi(x)]], x=-a..a, color=[red,blue], thickness=2);

Brian,  here is a  more programmatic solution that is written on the basis of your solution. This solution can easily be rewritten as a procedure that is suitable for the solution of any problem of this type.

restart:

with(combinat):

T:=[0.8, 0.7, 0.4]:

L:= [T[1]$3, T[2]$12,T[3]$15]:

P:=[seq(op(choose(L,m)), m=1..3)]: k:=0:

for i in P do

if `+`(op(i)) < 1.44

then k:=k+1; S[k]:=i; fi;

od:

S:=convert(S, list);

dcs:=add(x[i], i=1..nops(S));

for i to 3 do

eq[i]:=0:

for j to nops(S) do

for k in S[j] do

eq[i]:=eq[i]+`if`(k=T[i], x[j], 0);

od: od: od:

for i to 3 do

eq||i=eq[i];

od;

Optimization[LPSolve](dcs, [eq[1]=3, eq[2]=12, eq[3]=15], assume=[integer, nonnegative]);

My solution above gives the correct answer, but it is not quite correct, because do not cover all variants for the location of files on disks. Here is the corrected solution.

There are 8 different variants:

1.  {0.8} will be  x  disks.

2.  {0.4, 0.8} will be  y  disks.

3.  {0.7} will be  z  disks.

4.  {0.7, 0.7} will be  u  disks.

5.  {0.4, 0.7} will be  v  disks.

6.  {0.4} will be  w  disks.

7.  {0.4, 0.4} will be  s  disks.

8.  (0.4, 0.4, 0.4} will be  t  disks.

Optimization[LPSolve](x+y+z+u+v+w+s+t, {x+y = 3, z+2*u+v = 12, y+v+w+2*s+3*t = 15}, assume = {integer, nonnegative});

                                         [13, [s = 0, t = 4, u = 6, v = 0, w = 0, x = 0, y = 3, z = 0]]

Only 4 variants are possible: {0.8, 0.4} will be x disks, {0.7, 0.7} will be y disks, {0.7, 0.4} will be z disks, {0.4, 0.4, 0.4} will be u disks.

Optimization[LPSolve](x+y+z+u, {x >= 3, 2*y+z >= 12, x+z+3*u >= 15, assume = {integer, nonnegative});

                                                     [13, [u = 4, x = 3, y = 6, z = 0]];

To paint each orbit the desired color, you must first make a list of the colors, and then in a loop to invoke it:

Colors:=[color1, ... , color9]:

for j to 9 do
 pic[j]:=spacecurve([subs(E[j]=E, x[j]),subs(E[j]=E,y[j]),subs(E[j]=E,z[j])],
                    E=0..2*Pi, color=Colors[j]):
od:

Addition:  After setting the colors, I recommend in the final images  separately to represent the orbits of  the first 4 planets, else on the common picture you can not see them, because their orbits are very different from  the orbits of more distant planets:

display([op(pic1[1..4]), pic1[10], op(pic1[11..14])],scaling=CONSTRAINED);

CR := (a, b, c, d)->(a-c)*(b-d)/((a-d)*(b-c));

e1 := CR(b, a, c, d);

isolate(CR(a, b, c, d) = lambda, a-c);

subs(%, e1);

t := [seq(0+i*(2*evalf(Pi)*(1/10)), i = 0 .. 10)]:

x := [3, 4, 2*evalf(Pi), 7.83]:

y := [2.5, 4.3, 6, 2*evalf(Pi)+2]:

z := [-2, .3, 1.5, evalf(Pi)]:

S := [seq([seq([x[i]*cos(t[j]), x[i]*sin(t[j]), y[i]*cos(t[j]-z[i])], i = 1 .. 4)], j = 1 .. 11)]:

plots[surfdata](S);

Maple easily solves your system with symbolic parameters. For specific solutions and plotting I arbitrarily gave numerical values of the parameters.

restart;

dsolve({diff(V(t)*C(t), t) = G - K*C(t), diff(V(t), t) = alpha - beta, V(0)=V0, C(0)=C0}, {V(t),C(t)});

assign(%);

Parameters:={G=10.5, K=2.1, alpha=6.3, beta=4.4, V0=2.5, C0=0.1}:

f1:=eval(C(t), Parameters);  f2:=eval(V(t), Parameters);

plot([f1, f2], t=0..5, color=[red, blue], thickness=2, legend=['C'(t), 'V'(t)]);

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.