eqs:=diff(Th(z, t), t) = 7.1428*(diff(Th(z, t), z))-1397941.885*(279-Tw(z, t))-0.2160487e-1*(diff(Th(z, t), z, z)),

diff(Tc(z, t), t) = -7.1428*(diff(Tc(z, t), z))+1298990.852*(Tw(z, t)-291)+0.189366e-1*(diff(Tc(z, t), z, z)), 

diff(Tw(z, t), t) = 3.3024901*(Th(z, t)-2*Tw(z, t)+Tc(z, t))+8.0029*10^(-4)*(diff(Tw(z, t), z, z)): 

bc:=Th(z, 0) = 296, (D[1](Th))(0, t) = 0, Th(1, t) = 296, 

Tc(z, 0) = 275, (D[1](Tc))(1, t) = 0, Tc(0, t) = 275,

Tw(z, 0) = 0, (D[1](Tw))(1, t) = 0,Tw(0, t) = 0: 

sol:=pdsolve({eqs},{bc}, numeric):

Using the resulting  module  sol , you can plot functions one or two variables.

Examples:

p1:=sol:-plot(Th(z,t), t=0..1, z=0.1, numpoints=1000):

p2:=sol:-plot(Th(z,t), z=0..1, t=0.3, numpoints=1000):

p3:=sol:-plot3d(Th(z,t), z=0..1, t=0..1, numpoints=5000):

plots[display](p1, axes=normal);

plots[display](p2, axes=normal);

plots[display](p3, axes=normal);

Rit, replace the last two lines of your code by the lines

f := subs([seq(originvarslist[i] = varslistm[i], i = 1 .. nops(originvarslist))], f);

f := subs([seq(varslistm[i] = varslist[i], i = 1 .. nops(originvarslist))], f);

Certainly Acer's ways are much shorter.

``(16)/``(4);

expand(%);

or

`16`/`4`;

parse(numer(%))/parse(denom(%));

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)]);