It looks like what you want is to control the accuracy of the result, rather than the precision used in the computations.  In general this might not be easy to do (Mathematica, I am told, tries to do this, but is not always successful).  There is no simple way to do it in Maple in general. 

You should heed the warning.  Without the restrictions of positivity, the problem is indeed unbounded.
For example, from any feasible solution you can add -2, 0 and 1 to x, y, z respectively and have another solution with a lower objective value.  The solution Maple returns is feasible, but of course not optimal.  I imagine it is the last solution Maple obtains before discovering that the problem is unbounded.

When you put in a delay condition, you don't have a differential equation any more, you have a delay-differential equation.  As far as I know, there are no facilities in Maple specifically for dealing with these, but there are various approaches you can use.  For example, see
http://www.mapleprimes.com/questions/39942-Differential-Equations-With-Delay

> simplify(convert(z, polar));

Only one answer, but note that in Maple w^(1/4) is the "principal branch" fourth root of w, not a multivalued fourth root.  The other branches are, of course, this one multiplied by fourth roots of unity.

I assume you're talking about the Student[Calculus1] package.  I think you may find it a useful tool (yes, it does have step-by-step facilities for limits, differentiation and integration: see ?Student/Calculus1/SingleStepOverview ), but it will not teach you calculus by itself.  More likely, the main thing you need is a textbook.  You might try e.g. Strang's Calculus or Keshet's Math 102 and Math 103 notes, all freely available on-line.  

Hints:

1) The answer should be 1 + (number of diagonals) + (number of interior intersections of diagonals).

2) If the vertices are numbered 1 to n (either clockwise or counterclockwise around the polygon), the diagonal from vertex i to vertex j (1 <= i < j <= n) has an interior intersection with all and only the diagonals where one vertex number is strictly between i and j and the other is less than i or greater than j.

i) The Binomial distribution says the probability of exactly s successes in k Bernoulli trials, with probability p of success in each one, is binomial(k,s)*p^s*(1-p)^(k-s).  The case he's looking at is s=1, p=1/6.  For a fair coin you'd have p=1/2.

ii) If f(k) = k*p*(1-p)^(k-1), f(k+1)/f(k) = (k+1)*(1-p)/k.  The maximum (for integer k) will be the first time this <= 1,
namely k = ceil((1-p)/p).

You didn't attach any code. But here's a simple way to do it.

> with(LinearAlgebra):
   M:= [seq(RandomMatrix(4,4), i=1..8)];
   for p in combinat[permute](8) do
      M[p[1]] . M[p[2]] . M[p[3]] . M[p[4]] . M[p[5]] . M[p[6]] . M[p[7]] . M[p[8]]
   end do:

This took over 555 seconds of CPU time on my computer.

A slightly more sophisticated way would be to build up a tree of products.  This would reduce the number of matrix
multiplications needed from 7*8! = 282240 to 8!/6! + ... + 8!/1! = 69272

> for i from 1 to 8 do P[[i]]:= M[i] end do:
for r from 2 to 8 do
  for p in combinat[permute](8, r) do
     P[p]:= P[p[1..r-1]] . M[p[r]]
  end do
end do:

This took just under 40 seconds.

Ah, I think I see the problem.  It is not dT/dx that is continuous across the boundary, but the heat flux lambda*dT/dx.  I think one way to tackle this is by a piecewise transformation of the position variable x, so that a displacement by dx in region j corresponds to a change of dX = dx/lambda in the new variable X.  Then dT/dX is continuous across the boundary. 

See e.g. <http://www.mapleprimes.com/questions/35560-External-Image-As-Background-In-2d-Plot-#comment44271>
for how to get a .jpg image as a 2D plot.  Scale and translate as necessary using plottools[scale] and plottools[transform], and combine this with your own plot using plots[display].

I think you just need to define the PDE using piecewise.

PDE := diff(T(x, t), t) = a*(diff(T(x, t), x, x));
a:= piecewise(x<0.2, 1e-6, 5e-7);
L := .3; h := 10; lambda := 3; `epsilon;` := .85; sigma := 5.6697*10^(-8);IBC := {T(0, t) = 1000, T(x, 0) = 298, (D[1](T))(L, t) = -h*(T(L, t)-298)/lambda+`epsilon;`*sigma(T(L, t)^4-298^4)};
pds := pdsolve(PDE, IBC, numeric);

Your n-tuples [q[1],q[2],...,q[n]] run over the n-fold Cartesian product of [1,2,3].  So you could do something like this:

> with(combinat):
   T:= cartprod([[1,2,3]$n]):
   while  not T[finished] do
     L:= T[nextvalue]();
     s[op(L)]:= convert(L,`+`)
   end do:

Positive semidefinite matrices form a closed convex subset of the space of all n x n real symmetric matrices.  Thus in terms of the 10 parameters x[1] ... x[10] we have a closed convex subset C of R^10 consisting of those x for which A is psd.  If the n x n matrix B is not psd, then there is an n x n psd matrix C such that Trace(B.C) < 0 (in fact C can be of the form V . V^%T for some column vector V, but I don't think that's important here).  You might look up "semidefinite programming".  I think there may be effective numerical methods for tackling this problem, but I don't know if there's anything specific implemented within Maple.

Well, you could use piecewise, with one alternative being undefined:

> f:= x -> piecewise(x>=-Pi and x <= Pi, abs(x), undefined);

The bottom line is: it's best to use 1D Maple input rather than 2D math input.