I think basically you need to work on a frame-by-frame basis. 

display([
  seq(pointplot([[3, 2], [2, 5], [5, 2], [5, 5]],symbol=solidcircle,colour=red,transparency=A), 
        A = 0 .. 0.25, 0.05),
  seq(display([pointplot([[3, 2], [2, 5], [5, 2], [5, 5]],symbol=solidcircle,colour=red,transparency=A),
     pointplot([[2, 2], [4, 5], [3, 2], [3, 3]],symbol=solidcircle,colour=blue,transparency=A-0.3)]), 
        A = 0.3 .. 1, 0.05),
  seq(pointplot([[2, 2], [4, 5], [3, 2], [3, 3]],symbol=solidcircle,colour=blue,transparency=A-0.3), 
        A = 1.05 .. 1.3, 0.05)],
      insequence=true);

If plot3d(w, t = 0 .. 5, y = 0 .. 2) works, try something like

plot(eval(w, t=5), y = 0 .. 2);

I think you want to have the transpose of your Matrix in that PLOT3D, rather than the Matrix itself.  Try it with

PositionPlot := PLOT3D(CURVES(pos^%T), COLOR(RGB, 0, 0, 1), TEXT([1, 60, 1], 'pitch', ALIGNABOVE));

You should use LinearAlgebra, not linalg which is deprecated.  But it looks to me like the question was about a symbolic solution for arbitrary n1 and n2.  Maple can't really help much with getting the formula.  But one formula is

Abig^(-1) = << E | F>, <G | H>> where

E = (A11 - A12.A22^(-1).A21)^(-1),
F = -A11^(-1).A12.(A22 - A21.A11^(-1).A12)^(-1),
G = -A22^(-1).A21.(A11 - A12.A22^(-1).A21)^(-1),
H = (A22 - A21.A11^(-1).A12)^(-1)

assuming all matrices this takes the inverse of are invertible.

To test this out (I'll try the case n1 = 2, n2 = 3):

> with(LinearAlgebra):
   n1:= 2: n2:= 3:
   A11:= Matrix(n1,n1,symbol=a11);
   A12:= Matrix(n1,n2,symbol=a12);
   A21:= Matrix(n2,n1,symbol=a21);
   A22:= Matrix(n2,n2,symbol=a22);
   M:= <<A11 | A12>, <A21 | A22>>:
   E:= (A11 - A12 . A22^(-1) . A21)^(-1): 
   F:= -A11^(-1) . A12 . (A22 - A21 . A11^(-1) . A12)^(-1):
   G:= -A22^(-1) . A21 . (A11 - A12 . A22^(-1) . A21)^(-1):
   H:= (A22 - A21 . A11^(-1) . A12)^(-1):
   Mi:= <<E|F>, <G|H>>:
   map(normal, M.Mi);

                       [1    0    0    0    0]
                       [                     ]
                       [0    1    0    0    0]
                       [                     ]
                       [0    0    1    0    0]
                       [                     ]
                       [0    0    0    1    0]
                       [                     ]
                       [0    0    0    0    1]

pi/pi^2 should automatically simplify to 1/pi.  However, I suspect that one of these is Pi and the other is pi.

Remember that Maple is case-sensitive, and the constant 3.14159... is Pi, not pi.

@chiwhalee : I doubt that you can get Maple to do the double integral as-is, but if you reverse the order of integration it should work for suitable f (at least if f and its Fourier transform are both integrable).  For example:

> f:= x -> exp(-(x-x0)^2/m^2);
   int(int(f(x)*exp(I*p*x),x=-infinity..infinity), p=-infinity .. infinity);
   simplify(%) assuming m > 0;

                                     2
                                   x0
                           2 exp(- ---) Pi
                                    2
                                   m

The problem seems to be in evalhf, which produces an undefined result with some of the sinh factors (which produce numbers greater than 10^308 or so).  For example:

> evalhf(sinh(240*Pi)/sinh(400*Pi));

                        Float(undefined)

A simple work-around is to prevent evalhf from being used, by setting Digits higher than evalhf(Digits).

Warning: that will cause the calculation to take much longer. For example:

> Digits:= 16; 
  F:= Sum(10*sinh(n*Pi*(L-y)/L)*sin(n*Pi*x/L)*(2*(1+(-1)^(n+1)))/(sinh(n*Pi)*n*Pi),n=1..N);
  plot3d(eval(F,{L=5,N=600}),x=0..5,y=0..5,axes=box);

For example:

U:= sin(x)*sin(y);
plots:-densityplot(-U, x=0..Pi, y=0..Pi, colorstyle=HUE, scaletorange = -1 .. 0.4);

@Axel Vogt :

The justification for the decay of the tail integral can be done using integration by parts.

> with(IntegrationTools);
   J:= Int(sin(t)*t^(-b),t=T..infinity);
   Parts(J,t^(-b)) assuming b > 0;

                               infinity
                              /                 (-b)
                     (-b)    |          cos(t) t     b
             cos(T) T     -  |          -------------- dt
                             |                t
                            /
                              T

and the integral on the right can be estimated directly, as its absolute value is bounded by

> int(t^(-b)*b/t,t=T..infinity) assuming b>0,T>0;

                                 (-b)
                                T

This is an interesting suggestion, but it could cause problems if it was done automatically.

1)  GlobalOptimization is likely to take a very long time on problems involving more than 3 or 4 unknowns.

2) NonlinearFit would perform differently for different users (one who has GlobalOptimization and another who doesn't), and the reason for the difference would not be obvious.

Perhaps it would be better to provide it as an option.

Yes, it is possible.  For example:

> sys:= {diff(u(x,t),t$2) = diff(u(x,t),x$2) + v(x,t),
             diff(v(x,t), t$2) = diff(v(x,t),x$2) - u(x,t)};
   IBC:= {v(x,0) = 0, u(x,0) = 1, v(0,t) = 0, v(5,t) = 0, u(0,t) = 1, u(5,t) = 1};
   pdsolve(sys,IBC,numeric);

My suspicion is that your system may in some way be singular at t=0, so pdsolve can't solve for the derivatives it needs to.  But it's hard to say without looking at the actual equations and initial/boundary conditions.

In a PDE system for pdsolve, I think all dependent variables must be functions of the same independent variables.  So instead of vs(t), you'd need vs(z,t).  Of course it will end up not depending on z if its differential equation and initial condition don't depend on z.  The following seems to work (I set Fr, k, etc to 1).

> PDEs:= subs({Fr = 1, k = 1, m = 1, Apop = 1, Csat = 1, Vsist = 1, vmed = 1},
             { diff(Cf(z,t),t)+vmed*diff(Cf(z,t),z)=k*Apop*(Csat-Cf(z,t))/Vsist,
                diff(mpop(z,t),t)+vs(z,t)*(diff(mpop(z,t),z))=-k*Apop*(Csat-Cf(z,t)),
                diff(vs(z,t),t)=Fr/m}); 
   IBCs:={Cf(z,0)=0,Cf(0,t)=0,mpop(z,0)=0,mpop(0,t)=0,vs(z,0)=0};
   Sol := pdsolve(PDEs, IBCs, numeric, time=t, range=0..1);

Your [x,y] values [3,9], [5,10], [1,-2], [-8,3] and [7,5] don't seem to be organized in any particular way, and 5 (the number of points) is prime so you can't have a regular grid.  Well, you can (by hand) make them into triangles:

> a:= <3,5,1,-8,7>: 
   b:= <9,10,-2,3,5>:
   c:= <10, 5, 2, 7, 8>:
   Pt:= [seq([a[i],b[i],c[i]], i=1..5)];
   PP:= [[Pt[1],Pt[2],Pt[5]],[Pt[1],Pt[3],Pt[5]],[Pt[1],Pt[3],Pt[4]]];
   plots:-polygonplot3d(PP, style=patchcontour, axes=box, labels=['a','b','c'],
     orientation=[-40,70]);

Is that what you had in mind?

I don't think you mean a contour plot in the sense Maple uses that word.  You seem to have points lying (approximately) on a curve.

> plots[spacecurve]([seq([a[i,1],b[i,1],c[i,1]],i=1..12)], axes=frame, labels=["a","b","c"]);

You can add arbitrary text in arbitrary places using textplot in the plots package, although it's a fair amount of work to get everything placed right.