Carl, your way from a mathematical point of view is incorrect. If a finite set of points on the curve lies in a plane, it does not mean that the entire curve lies in this plane.

Moreover  geom3d[AreCoplanar]  command can cause errors in evident examples. Look 

with(geom3d):
AreCoplanar(point(A1,[1,0,0]), point(A2,[2,0,0]), point(A3,[3,0,0]), point(A4,[1,1,0]));

    Error, (in geom3d:-plane) the points may not be AreCollinear

For a curve given by the parametric equations with smooth functions, we can use  Student[VectorCalculus][Torsion]  command.  The quote from Wiki "A plane curve with non-vanishing curvature has zero torsion at all points. Conversely, if the torsion of a regular curve with non-vanishing curvature is identically zero, then this curve belongs to a fixed plane."

Example:

A:=<cos(phi),0,-sin(phi); 0,1,0; sin(phi),0,cos(phi)>:

V:=<cos(t),sin(t),0>:

V1:=A.V;  # Circle  V  rotated through the angle phi around the axis Oy

Student[VectorCalculus][Torsion](V1, t);

After  assign(a, c);  the following two assignments  

assign(a, d);  and   a:=d;  are different things. Compare:

assign(a, c);

a := d;  b := 2;

a, b, c;

    a := d

    b := 2

    d, 2, c

plots[polarplot](1/(theta-(1/3)*Pi), theta = 1.5 .. 8);

Edited: this is a spiral. The range for the theta should be taken such that  theta-alpha<>0

The simplest way to make the bolding (input or output, or all together) is to select the desired thing by the mouse and press  B .

Addition: if you want to make the bolding throughout the worksheet, then instead of a mouse,  first select all in edit menu.

collect(algsubs(1/(T[1]*V[2])=U, B), U);

Only in standard interface:

plots[polarplot]((1/6)*Pi, r = 0 .. 3, ordering = [angular, radial]);

More short code:

L:=select(t->`*`(op(t))<>0 and nops(convert(t, set))=9, [seq(seq(seq([x, y, z, 1+x, 2+y, -3+z,-2+x, 3+y, -1+z ], z=1..9), y=1..9), x=1..9)]);

nops(L);

map(t->{a=t[1..3], b=t[4..6], c=t[7..9]}, L);

1) In Maple  e^x  should be coded as  exp(x)

2) I think that in the general case, your problem can be solved only numerically.

Example:

G1 := (lambda, k, t) -> Int(exp((lambda - k)*s - exp(- k*s)), s = -infinity .. t);

evalf(G1(4, -3, 10));

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