surfdata is only for points organized in a more-or-less regular grid, with n points in each of m columns.  Looking at your points, they don't seem to be organized in that way at all.  You can use pointplot3d to plot the actual points.  Or you can try the Surplot package, which you can download from the links at

<http://www.mapleprimes.com/questions/80862-Plot-A-3d-Surface-With-Data-Points>

I only see two equations in what you wrote:

-28 + x(1) + x(2) + x(3) + x(4) + x(5) + x(6) + x(7) + x(8)

   + x(9) + x(10) + x(11) + x(12) + x(13) + x(14) + x(15) + x(16)

   + x(17) + x(18) + x(19) + x(20) + x(21) + x(22) + x(23)

   + x(24) + x(25) + x(26) + x(27) + x(28) = 0

and

0.003000000000 x(1) + 0.004000000000 x(2) + 0.006000000000 x(3)

   - 0.006000000000 x(4) + 0.004000000000 x(5)

   + 0.001000000000 x(6) + 0.001000000000 x(7)

   + 0.004000000000 x(9) + 0.008000000000 x(10)

   + 0.006000000000 x(11) + 0.01700000000 x(12)

   + 0.002000000000 x(13) + 0.01300000000 x(14)

   + 0.01200000000 x(15) + 0.006000000000 x(16)

   + 0.001000000000 x(17) + 0.002000000000 x(18)

   + 0.001000000000 x(19) + 0.01200000000 x(20)

   - 0.001000000000 x(21) + 0.004000000000 x(22)

   + 0.01700000000 x(23) + 0.004000000000 x(24)

   + 0.009000000000 x(25) + 0.006000000000 x(26)

   + 0.01000000000 x(27) + 0.01600000000 x(28) - 0.1120000000 = 0.

You seem to have 28 variables x(1) to x(28) and only two equations, so of course you won't have unique solutions.  But I suspect you aren't telling us the whole story.  What is the actual problem?

You might also try DEplot in the DEtools package.  For example:

> with(DEtools): de := diff(y(x),x,x) + 4*y(x) = sin(2*x);
   DEplot( de, y(x), x=0..3*Pi, y=-3..3,[seq([y(0)=0,D(y)(0)=i],i=-3..3)],linecolour=[red,orange,yellow,green,blue,cyan,magenta],numpoints=100);

I don't know why you want to use powcreate for this.  That's a command for working with formal power series (where the coefficients are algebraic expressions, not lists).  And the equations defining the coefficients must be algebraic, not procedures.  Such constructions as subsop(1 = NULL, u(n-1)) don't work the way you expect them to when u(n-1) is a symbolic term: it will just produce u(), because the first operand of u(n-1) is n-1.

Well, I guess you could say you take the formal power series in two variables

f(x,y) = sum(sum(binomial(n,k)*x^k*y^n, k=0..n), n = 0 .. infinity) = sum((1 + x)^n y^n, n=0 .. infinity).

Something like this:

> with(powseries):
   powcreate(u(n)=u(n-1)*(1+x),u(0)=1):
   expand(tpsform(u,y,7));

1+(1+x)*y+(1+2*x+x^2)*y^2+(1+3*x+3*x^2+x^3)*y^3+(1+4*x+6*x^2+4*x^3+x^4)*y^4+(1+5*x+10*x^2+10*x^3+5*x^4+x^5)*y^5+(1+6*x+15*x^2+20*x^3+15*x^4+6*x^5+x^6)*y^6+O(y^7)

I think you want dualaxisplot in the plots package.  This was introduced in Maple 12, and requires the Standard interface.

> dualaxisplot(aa,bb);

According to the general theory of differential equations, a second-order ODE such as this with initial conditions specifying y(0) and D(y)(0) will have a unique solution on some interval containing 0.  Are you trying to determine a value of the parameter a such that y(2) = 0?

> ode:=diff(y(x),x$2) -( a*y(x) - y(x)^3 + y(x)^5 ) = 0;
   S:=dsolve({ode,y(0)=1.31, D(y)(0)=0},y(x),numeric,parameters=[a],output=listprocedure);
   Y:= subs(S,y(x));
   YY:= proc(t) Y(parameters=[a=t]); Y(2) end proc;
   plot(YY, -10 .. 10);

> fsolve(YY, -7 .. -5);

                     -6.127404760

> fsolve(YY, -2 .. 0);

                      -1.388821449

Assuming S has an asymptotic series in inverse powers of Pec_i, we can find the coefficients.

>  A:=Theta0/C=(exp(Pec_i*k*S/(C*(R*S+1-S)))-1)*(H+1/(1-exp(-Pec_i*(1-S)/(R*S+1-S))));
    assume(C>0,H>0,R>0,k>0,Theta0>0);
    A1:= eval(lhs(A)-rhs(A),S=add(c[j]/Pec_i^j,j=1..6)+O(1/Pec_i^7));
    A2:= asympt(A1, Pec_i, 7);

Note that the result has terms in 1/Pec_i^j for j from 0 to 6, plus O(1/Pec_i^7).  To extract the coefficients, it is convenient to substitute Pec_i = 1/t and make the result a polynomial.

>  A3:= convert(subs(Pec_i=1/t, A2), polynom);
    eqs:= {coeffs(A3,t)};
    sol:= solve(eqs,{seq(c[k],k=1..6)});
    factor(sol);

Changing that sign in your solution really means changing q to -q.  Note that there is no q in the differential equation, just in one of the initial conditions.  Indeed your changed solution satisfies D(T)(0) = q/k_m rather than -q/k_m.

In general, solving equations (e.g. using solve) would not make much use, if any, of assumptions.  You might use the assumptions to show that there exists a solution in some region, but Maple wouldn't use them.  What exactly are you trying to do?

No, sorry, Maple does not understand English (or any other human language).  Translating a "word problem" into mathematics is still a job for humans. 

Basically it looks like Maple is trying a simplification involving a square root under certain assumptions (e.g. sqrt(x^2) simplifies to x if x >= 0).  It could be that an attempt to check whether some expression is nonnegative leads to an expression so large that Maple doesn't have enough memory for it.  However, there may be some way around the problem.  I can't tell without seeing the actual Maple commands that led to this.  You might upload a worksheet file.

Most integrals of complicated algebraic functions have no closed-form solution.  This seems to be an example.  integrals.com also fails to find a solution.  Do you have any reason to think there is one?

Do you mean something like this?

Gama:= (i,k) -> ((x[k]-x[i])*xdot[k]+(y[k]-y[i])*ydot[k])*dtheta[i]
    /((x[k]-x[i])^2+(y[k]-y[i])^2)^2; 
G:= i -> add(Gama(i,k), k = {$1..100} minus {i});

There are so many things wrong with your

 sum(1-w)^n-1,k=3..5);

it's hard to know where to begin.

1) Who said anything about a sum?

2) Be careful with parentheses.  What are you taking a power of? 
Also, ^ has priority over -, so a^b-c means (a^b)-c, not a^(b-c).

3) What is that k doing?

1)

> coeftayl((x+y+z)^4, [x,y,z]=[0,0,0], [1,1,2]);

            12

2)  I suspect what you meant was

(1+I*3)^2/(4-I*6)^7 

(note that sqrt(-1) is I, not i, and you want multiplication signs).

In this case

evalf[5](Re((1+I*3)^2/(4-I*6)^7));

does work, but I would advise you _not_ to do it this way: evalf[5] doesn't mean show the result to 5 significant digits, it means do the calculation with 5 significant digits, and in some cases this can lead to severe roundoff-error problems.  Better is to first use evalf with the default 10 digits, and only then round to 5 digits.  So:

evalf(Re(1+I*3)^2/(4-I*6)^7)):
evalf[5](%);