Your equation with initial conditions has no solution, because at  x = 0  the left side of the equation is non-negative, and the right side is negative.

shift:=proc(n, L)

if n=1 then return [L[2],L[3],L[1]] elif

n=2 then return [L[3],L[1],L[2]] elif

n=3 then L fi;

end proc:

Examples of use:

shift(1, [a,b,c]),  shift(2, [a,b,c]),  shift(3, [a,b,c]);

                     [b, c, a], [c, a, b], [a, b, c]

I reduced the range of integration over the variable  y , otherwise under the sign of the square root  a negative number is.

evalf(Int((1-exp(-5.5/cos(x)))*sin(x), [x = 0 .. arctan(300*cos(y)+sqrt(12.25-90000*sin(y)^2)), y = 0 .. arcsin(7/600)], method = _Gquad));

                                                             0.01162184070

See help on  evalf/Int  

Example:

plots[display](plot(x^2, x = -1 .. 2, color = red, thickness = 2), plots[textplot]([1, 2, y = x^2], font = [TIMES, ROMAN, 18]));

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: