It seems to be attempting to permute the colour components in a plot structure.  For example, Cperm([3,1,2]) would take a plot structure using RGB (red-green-blue) colour specifications and make red into green, green into blue, blue into red.  However, it doesn't work correctly, because the substitutions are done sequentially (i.e. after one object's colour has been changed, it might be changed again if there is another object whose original colour matched this object's new colour).

 Hmmm, I tried to do this using evalindets or subsindets, but that suffers from the same problem, which it seems to me is a bug.  For example:

> evalindets([A,B], 
   proc(t) 
    if t=A then B+1 elif t=B then C+1 elif t=C then A+1 end if 
   end proc);

[C+2, C+1]

No, a..b x c..d means the rows are numbered from a to b and the columns are numbered from c to d.  If you always want the numbering to start with 1, you might use a Matrix rather than an Array: thus Matrix(3,1) would be a Matrix with three rows (numbered 1, 2, 3) and one column.

> solve(b*(b-4) < -4*a);

Why would you use evala?  That's for algebraic number fields. You might try normal or factor.

Try translate and rotate from the plottools package.  For example,

> with(plottools): with(plots):
  display([D0, rotate(translate(D0, 2, 3),  Pi/6, [2,3])]);

Try this:

> U:= subs(pds:-value(output=listprocedure), u(x,t));

 Then use U(x,t) to get the value of u(x,t).

Apologies for starting a new topic, but somehow the original one seems to be "locked", so replies are not possible. 

That imaginary part is on the order of 10^(-12).  It's not particularly significant.

In essence, here's what's happening.  I'll use somewhat simpler numbers, and change the variable to x.

You have an indefinite integral

> J1:= int(1/(sin(x+1/2)+1/10),x);

This antiderivative is perfectly good (where it is continuous), but it happens to have a nonzero imaginary part on some intervals.  Then

> v:= allvalues(solve(J1 - eval(J1, x=0) - 1, x));

> evalf(v);

.8473702770+.3178782144e-48*I

The tiny imaginary part here is not significant.  It is the result of roundoff error.  Even though the final result ought to be real, some of the immediate expressions are not, e.g.

> evalf(arctanh(1/33*(tan(1/4)+10)*11^(1/2)));

That's perfectly normal and no problem.  Somewhat more serious is the fact that v is not at all a solution of the equation you're trying to solve (presumably due to problems with branches of arctan and arctanh).

> evalf(eval(J1 - eval(J1, x=0) - 1, x = v));

On the other hand, it is a solution of the equation you probably meant to solve:

> evalf(Int(1/(sin(x+1/2)+1/10),x=0..v));

> with(plots):
  listplot([seq(C[i][1,1],i=1..10)], style=point);

pdsolve/numeric only deals with equations involving two independent variables (one time variable and one space variable), so I'm not sure what you would want the 3-d animation to represent.  Well, perhaps you have two dependent variables.

> Sys:= {diff(u(x,t),t) = v(x,t) + diff(u(x,t),x,x),
       diff(v(x,t),t) = -u(x,t) + diff(v(x,t),x,x)};
  IBC:= {u(x,0)=x-x^2, v(x,0)=1+x-x^2, u(0,t)=sin(t),v(0,t)=0,u(1,t)=0,v(1,t)=cos(t)};
  Sol:= pdsolve(Sys,IBC,numeric);

  Q:=Sol:-animate([x,u,v],x=0..1,t=0..2*Pi,axes=box);  
  

A look at the structure of Q shows that it consists of 16 frames, each of which contains a list of three Arrays of size 1..21 x 1..2.  The second component of each Array contains the values of x, u, v respectively at a particular time value.
 

> for i from 1 to 16 do
    L:= map(op,op([1,i],Q));
    A:= Array(1..21,1..3,datatype=float[8]);
    for j from 1 to 3 do A[1..21,j]:= L[j][1..21,2] od;
    frame[i]:= PLOT3D(CURVES(A),COLOR(RGB,0,0,0),TITLE(t=evalf((i-1)*2/15,3)*Pi))
  end do:
  plots[display]([seq](frame[i],i=1..16), insequence=true, axes=box, labels=[x,u,v]);

I tried making this into a system, with D[1](u) another dependent variable v.  However, there seem to be severe numerical problems with v at the boundaries.

> Sys:= {diff(u(x,t),t) = v(x,t)+diff(v(x,t),x), diff(u(x,t),x)=v(x,t)};
  IBC:= {u(x,0)=1,u(0,t)=0,u(1,t)=0};
  Sol:= pdsolve(Sys,IBC,numeric,time=t,timestep=0.01,spacestep=0.001);
  Sol:-plot(v, x=0, t=0..1,style=point);

Executing this in Maple 13, the first error I get is "invalid matrix/vector product", and the offending piece of code has a red dotted line around it.

Note the dot between - and tan.  That's what's causing the error.

Another mistake I see is that you have int(eqw1, bX .. cX, numeric) rather than int(eqw1, x = bX .. cX, numeric).

I hope you used ln rather than LN. 
NLPSolve does not know or care whether the function is convex.  It finds a local maximum or minimum.  If the objective and feasible region are convex, a local minimum is guaranteed to be a global minimum; otherwise, the solution NLPSolve finds is not guaranteed to be a global minimum.
However, you say you want to maximize your objective.  For a local maximum to be guaranteed to be global, you want the objective function to be concave, not convex.

Actually I get 105 solutions.

> with(Optimization): 
  cons:= {seq(add(A[i][j]*x[i],i=1..nops(A))=1,j=1..8)}:
  avoidcons:= {}:
  sols:= NULL;
  do
    try
      sol:= Maximize(0, cons union avoidcons, assume=binary)[2];
    catch:
      break
    end try;
    sols:= sols,select(t -> subs(sol,x[t])=1,{$1..28});
    avoidcons:= avoidcons union
       {add(subs(sol, x[i])*(1-2*x[i])+x[i], i=1..28)>=1};
  end do:  
  sort([sols],(a,b) -> add((a[i]-b[i])*30^(4-i),i=1..4)<0);

[{1, 14, 23, 28}, {1, 14, 24, 27}, {1, 14, 25, 26}, {1, 15, 20, 28}, {1, 15, 21, 27}, {1, 15, 22, 26}, 
{1, 16, 19, 28}, {1, 16, 21, 25}, {1, 16, 22, 24}, {1, 17, 19, 27}, {1, 17, 20, 25}, {1, 17, 22, 23}, 
{1, 18, 19, 26}, {1, 18, 20, 24}, {1, 18, 21, 23}, {2, 9, 23, 28}, {2, 9, 24, 27}, {2, 9, 25, 26}, 
{2, 10, 20, 28}, {2, 10, 21, 27}, {2, 10, 22, 26}, {2, 11, 19, 28}, {2, 11, 21, 25}, {2, 11, 22, 24}, 
{2, 12, 19, 27}, {2, 12, 20, 25}, {2, 12, 22, 23}, {2, 13, 19, 26}, {2, 13, 20, 24}, {2, 13, 21, 23}, 
{3, 8, 23, 28}, {3, 8, 24, 27}, {3, 8, 25, 26}, {3, 10, 16, 28}, {3, 10, 17, 27}, {3, 10, 18, 26}, 
{3, 11, 15, 28}, {3, 11, 17, 25}, {3, 11, 18, 24}, {3, 12, 15, 27}, {3, 12, 16, 25}, {3, 12, 18, 23}, 
{3, 13, 15, 26}, {3, 13, 16, 24}, {3, 13, 17, 23}, {4, 8, 20, 28}, {4, 8, 21, 27}, {4, 8, 22, 26}, 
{4, 9, 16, 28}, {4, 9, 17, 27}, {4, 9, 18, 26}, {4, 11, 14, 28}, {4, 11, 17, 22}, {4, 11, 18, 21}, 
{4, 12, 14, 27}, {4, 12, 16, 22}, {4, 12, 18, 20}, {4, 13, 14, 26}, {4, 13, 16, 21}, {4, 13, 17, 20}, 
{5, 8, 19, 28}, {5, 8, 21, 25}, {5, 8, 22, 24}, {5, 9, 15, 28}, {5, 9, 17, 25}, {5, 9, 18, 24}, 
{5, 10, 14, 28}, {5, 10, 17, 22}, {5, 10, 18, 21}, {5, 12, 14, 25}, {5, 12, 15, 22}, {5, 12, 18, 19}, 
{5, 13, 14, 24}, {5, 13, 15, 21}, {5, 13, 17, 19}, {6, 8, 19, 27}, {6, 8, 20, 25}, {6, 8, 22, 23}, 
{6, 9, 15, 27}, {6, 9, 16, 25}, {6, 9, 18, 23}, {6, 10, 14, 27}, {6, 10, 16, 22}, {6, 10, 18, 20}, 
{6, 11, 14, 25}, {6, 11, 15, 22}, {6, 11, 18, 19}, {6, 13, 14, 23}, {6, 13, 15, 20}, {6, 13, 16, 19}, 
{7, 8, 19, 26}, {7, 8, 20, 24}, {7, 8, 21, 23}, {7, 9, 15, 26}, {7, 9, 16, 24}, {7, 9, 17, 23}, 
{7, 10, 14, 26}, {7, 10, 16, 21}, {7, 10, 17, 20}, {7, 11, 14, 24}, {7, 11, 15, 21}, {7, 11, 17, 19}, 
{7, 12, 14, 23}, {7, 12, 15, 20}, {7, 12, 16, 19}]

Actually you have one boundary condition Sigma(0,t) = 10 and one initial condition.Sigma(R,0) = 0.  Since your pde is second order in R, you need a second boundary condition.  For example, this would work:

> PDE := diff(Sigma(R, t), t) = 3*(diff(R^.5*(diff(Sigma(R, t)*R^.5, R)), R))/R;
  IBC := {Sigma(0, t) = 10, Sigma(10, t) = 0, Sigma(R,0) = 0};
  pds := pdsolve(PDE, IBC, numeric, Sigma(R, t), time = t, range = 0 .. 10);

Note that to match your PDE and IBC, it must be Sigma(R, t), not Sigma(t, R).  You don't need to specify indepvars.

pdsolve is for problems in finite spatial domains.  It won't handle boundary conditions at infinity.  However, you might try a transformation that takes an infinite interval into a finite interval.  For example:

> with(PDEtools):
    simplify(dchange(R=tan(r), PDE, [r,t]) assuming r > -Pi/2, r < Pi/2);

Dirac deltas can't be done per se in pdsolve.  You might try an "approximate delta": something that is nonzero only on a small interval near 0 and has integral 1.  This may be tricky, however, because you'll need a very fine spatial resolution (a small value of the spacestep option).

If P is your plot,

> op([1,1],indets(P,specfunc(anything,CURVES)));