Sometimes Maple seems to be a bit stupid, so help it

Int(Int((r^2+r*cos(theta))*r*(a^2-r^2)^(1/2),r = 0 .. a),theta = 0 .. 2*Pi);
value(%) assuming (0 < a, a < 2*Pi);
         2 Pi    a
        /       /
       |       |     2                     2    2 1/2
       |       |   (r  + r cos(theta)) r (a  - r )    dr dtheta
       |       |
      /       /
        0       0


                                  5
                               4 a  Pi
                               -------
                                 15

But - is that really stupid? Only in some sense, I think - but I am not sure
(though it always makes me grumbling): if you think complex, then it is not
clear to have an ordering, or if the path is a (real) line providing an order ...

So give merci to the generality and feed assumtions in case.

PS: providing input in ASCII would make it more easy to copy/paste for answering.

In addition to Duncan:

maplemint(interpol); # give it a try ...    

    This equation was used as a statement: 
      p1(x) = f(x0)+f[x0,x1](x-x0)
    This equation was used as a statement: 
      y[i] = (y[i]-y[i-1])/(x[i]-x[i-j])
    This equation was used as a statement: 
      p[k] = y[1]
    These names were used as global names, but were not declared: 
      x0, x1

For this example it seems to diverge, can not remember conditions for that method.

But what is the original task / exercise you want to answer?

Look up the help for ?evalc and ?Re, ?Im

PS: do not use [...] for those grouping, always use (...).
And do not hesitate to use *, the multiplication operator.

First I would care for correct dimensions in A.x=b, so it actually is a linear equation.

LinearSolve should find a solution - but there may be more than only just 1
(or otherwise said: why should D be invertible?)

Usually one can not see that in examples by random values, since generically
a matrix is invertible

Edited for pagan: the reason that a 'test' fails is, that  Det(A) = 0 is an (analytical)
meager set, thus meating it has prob = 0 (as far as I remember that stuff), it is
like having probability f(x) = 0 in dimension = 1, i.e. that exactly just 1 value is taken.

In this vague form certainly not: take 'any expression' and ask from which sum it comes,
Sum(expr(k), k=1 .. 1) is always true ...

Though it is hazardous: looking at sum as an analogy to int it would be asking for a DE,
variable could be the a or b. So finite differences *may* stay for DE.

Hm ... 'umbral calculus' comes to my mind (some kind of voodoo for me ...), may be you
look in that field to get ideas for a more solid statement of the problem.

 a:= (2*x-1)/sqrt(1-4*x^2); 
 b:=sqrt(1-2*x)/sqrt(1+2*x);
 
 'radsimp(a/b)': '%'=%;

                          radsimp(a/b) = -1

I think the problem comes atleast through your 'indNonNegative':

If J stands for your integrand, then it re-writes as

K:=piecewise((1/3)*x11+(1/3)*x12+(1/3)*x13-(1/3)*x21-(1/3)*x22-(1/3)*x23 < 0, 0, 0 <= (1/3)*x11+(1/3)*x12+(1/3)*x13-(1/3)*x21-(1/3)*x22-(1/3)*x23, 1)*piecewise((1/3)*x11+(1/3)*x12+(1/3)*x13-(1/3)*x31-(1/3)*x32-(1/3)*x33 < 0, 0, 0 <= (1/3)*x11+(1/3)*x12+(1/3)*x13-(1/3)*x31-(1/3)*x32-(1/3)*x33, 1);

       /{          x11   x12   x13   x21   x22   x23     \
       |{ 0        --- + --- + --- - --- - --- - --- < 0 |
       |{           3     3     3     3     3     3      |
  K := |{                                                |
       |{               x11   x12   x13   x21   x22   x23|
       |{ 1        0 <= --- + --- + --- - --- - --- - ---|
       \{                3     3     3     3     3     3 /

        /{          x11   x12   x13   x31   x32   x33     \
        |{ 0        --- + --- + --- - --- - --- - --- < 0 |
        |{           3     3     3     3     3     3      |
        |{                                                |
        |{               x11   x12   x13   x31   x32   x33|
        |{ 1        0 <= --- + --- + --- - --- - --- - ---|
        \{                3     3     3     3     3     3 /

which is a product (difficult to recognize and see it here) and one can achieve that by

theAssumptions:=NULL;
for t in {x11, x12, x13, x21, x22, x23, x31, x32, x33} do
  theAssumptions:=theAssumptions, 0 < t, t < 1;
end do:
theAssumptions:

K:=simplify(J) assuming theAssumptions;

You integrate over the unit hypercube, may be you have to transform boundaries
fitting the linear conditions above.

Otherwise any integrator has to check those conditions in every point of evaluation
and I am not even sure, whether the above is continous over the hypercube.

PS: only in rare cases you can copy from a *.mw in a reasonable way.

First I have to state that Maple 12 behave quite buggy for that *.mw for
the graphical interface: already for those expressions it clutters up the
display and I have to scroll up and down to correct it.

PS: Preben Alsholm has answered as well while typing this, I do not want
to update my writing and there may be similar thoughts about the task ...


Now for your task (and in a 3rd step you also may check what hirnyk said),
I take the notations from the uploaded sheets (to long to repeat it here),
especially I use upper case Int to avoid evaluations if not desired and
use Digits:=14 [or better write Digits:=round(evalhf(Digits)) ].

1. 
Looking at the result so far for 0 < t only constants come out and they
are either zero or of opposite sign, stemming from JL or JR Before.

Looking at those they *only* differ in sign, Maple shows that through

  'JLBefore + JRBefore': '%'= combine(%);

                           JLBefore + JRBefore = 0

For 0 < t those are zero (by your piecewise definition) and for t < 0 it
is an integral over the whole Reals, not involving t (giving a constant).

One can evaluate it very fast in 'double precision' using the NAG routines.

  'JLBefore assuming t < 0'; #lprint(%); # to copy for real work sheets ... 
  combine(%);                # to have all under the integral sign

  # now append the appropriate command enforcing NAG methods, see the help
  Int(op(%), method = _d01amc);
  evalf(%);

                              0.12075268864975

The answer is *immediate*.

2. 
For 0 < t only 2 of your integrals depend on t, the other 4 are constants,
either 0 (the 1st and the 3rd) or similar to above:

  'JLAfter + JRAfter': '%'= combine(%);

                            JLAfter + JRAfter = 0

The value also can be computed very fast like shown above.

3. 
All the above can be done before you even start to plot and there are only
two integrals to be plot, depending on t, 0 < t.

One of them is your JLPos, for which I want to go a similar way as above.

  JLPos assuming 0 < t:
  combine(%);

That gives me 2 integrals, where the 2nd is of form Real(Int(...)).

Starring at the integrand it is clear, that it is real valued (check it!)
and thus there is no need to use it - which allows me to combine those into
just 1 integral:

  JLPos assuming 0 < t:
  combine(%): 

  eval(%, Re = 'u -> u'): # replace Re by identity
  combine(%):             # everything under one integral
  simplify(%, size):      # to have a more simple integrand
  
  JLP:=Int(op(%), method = _d01ajc);

The last command again is just a way to use NAG procedures again (without
re-typing it): this integral is over a finite range, so we need the method
for that, see the help.

Then you can plot it quickly

  plot(JLP, t=0..3);



Certainly you can combine all that into a revised version of your solution
(but try to do it in 'worksheet' mode).

I think there are several issues and I have not gone through all.

Mainly I replaced evalf(int( ...)) by Int and let that do the plot,
except in those cases, where it is seen to give an constant.

I would suggest to care for real and imaginary on the integrands,
not on the integrals.

But I did not do: your sheet is in document mode and that is the
worst case for others to work with, one can not even copy or export
in a save way.

If you really are forced to work with the new *.mw instead of the
*.mws sheet, then at least switch to 'worksheet' mode.

QD_TH_TI_103_revise.zip

You need to know more about you parameters.

data:=[A=4.00E-03, B=-6.00E-07, C=0, W=2.6];
data:=map(convert, data, rational);

# final situation for t[n]=(W-1)/(A+B*t[n-1]+C*t[n-1]^2);
T=(W-1)/(A+B*T+C*T^2);

eval(%, data);
S:=[solve(%, T)]; #nops(%);
evalf(%);
                                 1/2                   1/2
                          2000 19               2000 19
          S := [10000/3 + ----------, 10000/3 - ----------]
                              3                     3

                      [6239.265962, 427.400704]

Now some iterations starting with t[0]=0 (like in the xls sheet)
idicate, that you have to take S[2].

Can one do that with some 'empty default plot' as well?

What I mean: the plot carrying the color in the given codes must have a range
and that should fit to the actual plot (more or less).

It would convenient if having such independently - is that possible?

W is the equation of an affine hyperplane (perhaps over the Reals or Rationals in dim=28).

For that you can find the orthogonal vector, which gives the full space.

From that conversely you can find the hyperplane as the orthogonal complement
plus the shift off the origin.

Also note that you can do this 'exact' and do not need numerical procedures, since you
actually have an (equivalent) system over the Rationals with integer coefficients:

  1000*W: 
  convert(%, rational);

  -112 - 6 x(4) + 2 x(18) + 4 x(5) + 4 x(22) + 17 x(23) + 12 x(20)
         - x(21) + x(19) + 4 x(24) + 8 x(10) + 10 x(27) + 16 x(28)
         + x(6) + 9 x(25) + 6 x(26) + x(17) + x(7) + 4 x(9) + 6 x(3)
         + 17 x(12) + 2 x(13) + 13 x(14) + 4 x(2) + 3 x(1) + 12 x(15)
         + 6 x(16) + 6 x(11)

May be that is the question (of some homework)?

................