It seems to me unlikely that closed-form formulas for these integrals exist.  If you give numerical values to all parameters, you can integrate numerically.  Thus (since gamma is a named constant for Maple, I'll use g instead):

> with(Statistics):
   h:= s + (v - s)/(1+g);
   phi:= unapply(PDF(Normal(mu,sigma),x),mu,sigma,x);
   f:= int(phi(h,sigma,x)*phi(g0,sigma0,g), g = -infinity .. infinity);
   M:= int(g*phi(h,sigma,x)*phi(g0,sigma0,g), g = -infinity .. infinity);
   evalf(eval(M/f, {x=1,s=2,v=1,sigma=1,sigma0=1,g0=1}));

                             1.059728453

This is a bug in Maple 12 (and 13).  It was fixed in Maple 14.  It works if you include `dsolve/numeric/data/modules` in the save command, but save to a text file rather than a ".m" file.  For example:

> save Sol, `dsolve/numeric/data/modules`,"c:/test/myfile.mpl";

I believe this can be done with JavaViewLib.  See <http://www.maplesoft.com/applications/view.aspx?SID=1503&L=G>

solve is still not very good at handling assumptions.  I think it's actually simplest to use solve here without any assumptions. 

> solve({mu=Mean(Y9),sigma^2=Variance(Y9)},{a,b});
                        2
                      mu
  {a = ln(---------------------------),
                       2     2     2
          RootOf(-sigma  - mu  + _Z )

                                           2     2     2
                     2        RootOf(-sigma  - mu  + _Z )
        b = RootOf(_Z  - 2 ln(---------------------------))}
                                          mu


> allvalues(%);
             2                  1/2 1/2
           mu          1/2    %1
  {a = ln(-----), b = 2    ln(-----)   },
            1/2                mu
          %1

                   2                   1/2 1/2
                 mu           1/2    %1
        {a = ln(-----), b = -2    ln(-----)   },
                  1/2                 mu
                %1

                     2                    1/2 1/2
                   mu          1/2      %1
        {a = ln(- -----), b = 2    ln(- -----)   },
                    1/2                  mu
                  %1

                     2                     1/2 1/2
                   mu           1/2      %1
        {a = ln(- -----), b = -2    ln(- -----)   }
                    1/2                   mu
                  %1

             2     2
  %1 := sigma  + mu

And since we want a to be real and b positive, it's clear that the solution we want is the first one.

I think the problem is that dsol doesn't compute values of diff(a(t), t$2).  However, you can use the differential equations to express this in terms of quantities that it does compute.

> att:= subs(solve(sys, {diff(a(t),t$2), diff(phi(t),t$2)}), diff(a(t),t$2));
                                    2     2
  att := - 0.02500000000 (-1. phi[t]  %1 a  + 20. a[t] %1 phi[t] a

                     2            2
         + 8000. a[t]  + 120. a[t]  %1)/(a (3. %1 + 400.))

  %1 := exp(0.05000000000 phi)

> odeplot(dsol, [t, att], t=-50 .. 50);

The implicitplot3d command is for plotting a surface defined by a single implicit equation F(t,p,s) = 0.  For example,

implicitplot3d(x^2 + y^2 + z^2 - 1, x = -1 .. 1, y = -1 .. 1, z = -1 .. 1) would give you a sphere.  You don't seem to be talking about implicit equations.  x = F(t,p,s), y = G(t,p,s), z = R(t,p,s) would give you a parametric surface for each value of s.  You could try an animation, showing how that surface changes as s goes from 0 to 1:

> animate(plot3d, [[F(t,p,s), G(t,p,s), R(t,p,s)], t=0..1, p=0..1], s = 0 .. 1);

There's nothing specifically for a solid, but you could try plotting the six parametric surfaces for t=0, t=1, p=0, p=1, s=0, s=1 (which might form the boundary of your solid):

> display([ seq(plot3d([F(t,p,s), G(t,p,s), R(t,p,s)], p=0..1, s=0..1), t=0..1),   
                seq(plot3d([F(t,p,s), G(t,p,s), R(t,p,s)], t=0..1, s=0..1), p=0..1),
                seq(plot3d([F(t,p,s), G(t,p,s), R(t,p,s)], t=0..1, p=0..1),s=0..1)]);
              

Questions of convergence for power series in two variables are more complicated than for one variable.  However, consider the line y = -x, where l(x,-x) = sqrt(1+cos(x)^2).  For this function, the closest singularities to 0 in the complex plane are at x = (+/-) Pi/2 (+/-) ln(1+sqrt(2))*i, so the radius of convergence of its Taylor series around 0 is

sqrt((Pi/2)^2 + (ln(1+sqrt(2)))^2) = approximately 1.801171980

So along this line, the bivariate Taylor series can't converge past |x| = sqrt((Pi/2)^2 + (ln(1+sqrt(2)))^2) either.

You can import data directly from an Excel file.  For example:

> with(ExcelTools):
   R:= Import("c:/myfolder/mydata.xls");

This should produce a Matrix with as many rows and columns as you have in the file.  So for example if the X and Y data form all of columns A and B respectively, 

> X := R(.., 1); Y := R(.., 2);

For your pointplot, you need the plots package, and for NonlinearFit, you need Statistics.  The worksheet you uploaded did have

> restart:with(plots):with(Statistics):

but this line doesn't seem to have been executed.  Re-execute your worksheet, e.g. using the !!! button, and the pointplot will be produced:

> pp := pointplot(X, Y, symbolsize = 25, color = blue);

     pp := PLOT(...)

Because the plot structure has been assigned to a variable, you don't see the plot here, but evaluate pp by itself and you will see it.

On the other hand, the NonlinearFit command produces an error message:

Error, (in Statistics:-NonlinearFit) sizes of independent and dependent data do not match

The complaint is correct: your X has 12 entries while Y has just 11.  Add another entry to Y or remove one from X and it should work.

> J := Int(F(x,y), x = 0 .. 1);

Let u = cos(a*y) and x = t/a.  The integral becomes

> Jt:= IntegrationTools:-Change(eval(J,cos(a*y)=u),x=t/a,t);

                  /a                    (1/2)   
                 |   /     2           \        
                 |   \1 + u  + u cos(t)/      dt
                 |                              
                /0                              
                --------------------------------
                               a                

Maple can actually evaluate this, at least under appropriate assumptions on a and u.

> value(Jt) assuming u>0, a > 0, a < 2*Pi;

  /                          /                             (1/2)\     
  |/             2\          |   /1  \   (1/2) /    u     \     |     
  |\2 u + 2 + 2 u / EllipticF|sin|- a|, 2      |----------|     |     
  |                          |   \2  /         |         2|     |     
1 |                          \                 \1 + u + u /     /     
- |-------------------------------------------------------------- - 2 
a |                                  (1/2)                            
  |                      /         2\                                 
  \                      \1 + u + u /                                 

              (1/2)          /                             (1/2)\     
  /         2\               |   /1  \   (1/2) /    u     \     |     
  \1 + u + u /      EllipticF|sin|- a|, 2      |----------|     | + 2 
                             |   \2  /         |         2|     |     
                             \                 \1 + u + u /     /     

              (1/2)          /                             (1/2)\
  /         2\               |   /1  \   (1/2) /    u     \     |
  \1 + u + u /      EllipticE|sin|- a|, 2      |----------|     |
                             |   \2  /         |         2|     |
                             \                 \1 + u + u /     /

  \
  |
  |
  |
  |
  |
  |
  |
  /

If Maple returns a difficult integral unevaluated, it's quite likely that there is no closed-form solution.  We might be able to say more if you included the actual input.  For a numerical computation, you would need to supply numerical values for all parameters, not just assume one of them is positive.  In this case it seems Maple can't answer the question because you haven't told it enough to decide whether a certain expression is between 0 and a certain other expression.  This typically arises in a definite integral where Maple can find an antiderivative, but needs to determine whether the antiderivative (which might have singularities or branch cuts) is continuous on the interval in question.  If you want to live dangerously, you can try the continuous option.

For a series in two or more variables, try mtaylor.

The usual question would be to diagonalize the matrix by writing A = P D P^(-1).  If you really want A = P D P, the requirement for P to be invertible is not needed (if P was not invertible, P D P couldn't be invertible either).  The straightforward way to solve it is like this (I'll use d instead of D, because D is reserved for the derivative operator).

> A:= <<-3,0,0>|<0,-3,0>|<6,-6,3>>;
   P:= Matrix(3,3, symbol=p);
   d:= LinearAlgebra[DiagonalMatrix]([d1,d2,d3]);
   eqs:= convert(A - P . d . P, set);
   solve(eqs);

This produces six solutions, each with 2 to 4 arbitrary parameters.

Hint: see the help pages for ?solve and/or  ?LinearAlgebra/LinearSolve .  If you use solve, you may want to use eval or subs to put the answer in vector form. 

Interesting: your 2D math input in the assignment to X looks like one large left bracket and one small right bracket, but when you convert it to 1D Maple input it becomes three left and three right brackets.  I don't know how you produced this.  Of course, what you want is

> X := Vector([0,1.01,2.01,3.01,4.02,5.02,6.03,7.03,8.03,9.04,1.11,1.21], datatype=float[8]):

which you can type in just like that, in either 1D or 2D.  Your assignment to Y is also, it seems, not what it looks like: after converting to 1D Maple input, it looks like 

> Y:= Vector(""[""[[-4.77-11,1.33E-09,2.75E-09,4.23E-09,5.79E-09,7.24E-09,8.60E-09,1.02E-08,1.17E-08,1.32E-08,1.48E-08,1.65E-08]]""],datatype=float[8]):

Again, I don't know where the quotes or the extra brackets came from.  I'd consider it a bug that the 2D math input doesn't show what is wrong here.

Your third error is, of course, the result of the fact that X and Y haven't been assigned values due to the first two errors.

Maple has a MmaTranslator package which can (sometimes) translate Mathematica code into Maple.  However, it doesn't seem to work in this case.  Here's what I might do:

> L:= [seq([tt,piecewise(tt<=0,0,tt<0.0104167,5800,0)],tt=0..2.5, 0.0005)]:
   plot(L,view=[0..0.012, DEFAULT], colour=green);

You can define the pde and boundary conditions using symbolic constants, but if you want to solve numerically you'll have to substitute numerical values for them.  However, D is a bad choice to use for a name, because it is the derivative operator (which also occurs in your boundary conditions).  I'll use d instead.  By the way, you need an additional boundary condition.  For example:

> pde:= diff(f(x, t), t) =d*(diff(f(x, t), x, x));
   ibc:= {f(x, 0) = 0.25, D[1](f)(0, t) = (1-N)*sqrt((2*k)/t)/d, f(1,t) = 0.25};
   params1:= {d = 2, N = 0.3, k = 0.4};
   sol1:= pdsolve(eval(pde,params1), eval(ibc, params1), numeric);