In fact, there are many more solutions than indicated by vv, for example

restart;
sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0];
fsolve(eval(sys,[a=10,b=10.4]), {x=0..infinity,y=0..infinity,z=0..infinity});
;

                      sys := [x^2-x*y+y^2-1 = 0, -a^2+y^2-y*z+z^2 = 0, -b^2+x^2-x*z+z^2 = 0]
                                    {x =0.1989730189, y = 1.084528287, z = 10.49805888}

Further, we strictly show that there are positive solutions  (x, y, z)  for arbitrary a=b>=1. From considerations of continuity, this implies that for any number  a>=1 and a sufficiently small positive number c  for  (a, b=a+c)  there are also a positive solutions  (x, y, z) .

restart;
sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0];
Sol:=simplify(solve(eval(sys,[b=a]), [x,y,z], explicit));
a in solve(a>=1 and `or`(seq((eval(x,k)>0 and eval(y,k)>0 and eval(z,k)>0), k=%)));

 The result:                                          a in RealRange(1, Open(infinity))

Numerical experiments show that all solutions seem to be in the strip  1<=a<=b<a+1/2 . This does not mean that every point of this strip is a solution. This is only a necessary condition.

restart;
a := 1: b:= 2: l:= 1: alpha[1]:= 1: alpha[2]:= 3:
syspde:= [diff(u(x, t), t)-a+u(x, t)-u(x, t)^2*v(x, t)-alpha[1]*(diff(u(x, t), x$2)) = 0, diff(v(x, t), t)-b+u(x, t)^2*v(x, t)-alpha[2]*(diff(v(x, t), x$2)) = 0]:
Bcs:= [u(x,0)=1,v(x,0)=1,D[1](u)(-l, t) = 0,D[1](u)(l, t) = 0,D[1](v)(-l, t) = 0,D[1](v)(l, t) = 0]:
sol:=pdsolve(syspde, Bcs, numeric);

sol:-value(u(x, t)+diff(u(x, t), x)+diff(u(x, t), t), t=0.7)(0.5);

                                          sol:=module() ... end module
               [x = 0.5, t = 0.7, u(x, t)+diff(u(x, t), x)+diff(u(x, t), t) = 2.34305224942915]

restart;	
s := Sum(a*X[n]+b, n=1..N);
expand(value(s));
op(2,%)%+op(1,%); 

s := Sum(X[n]+Y[n], n=1..N);
(expand@value)(s);

I use Maple 2018.2

Do you mean this:

Example:

for n from 1 to 10 do
a[n]:=n^2;
od:

seq(a[n], n=1..10);
                                    1, 4, 9, 16, 25, 36, 49, 64, 81, 100

The reason for the lack of a solution can be not only a comma. You solve a linear system with 9 unknowns, but in the solve command you specify only 6 unknowns. Could well be a situation that the system is consistent, but the unknowns, that you specified, are not expressed in terms of the remaining unknowns. Here is a simple example:

solve({x-2*y+3*z=0,2*x-4*y-z=1},{x,y,z});
solve({x-2*y+3*z=0,2*x-4*y-z=1},{x,y});
solve({x-2*y+3*z=0,2*x-4*y-z=1},{x,z});

                    {x = 2*y+3/7,  y = y,  z = -1/7}
                                        NULL
                          {x = 2*y+3/7, z = -1/7}

In the first case, Maple returns a general solution. We see that there are infinitely many solutions (a straight line of the solutions). But in the second case Maple does not return anything (NULL), because {x,y} cannot be expressed through  z (the line is perpendicular to  Oz-axis). The third case works again.

The other two answers consider only simple examples in which one can easily find a parameterization or one variable is explicitly expressed through another. Here is a more complicated case in which these methods do not work. Here, we first plot the curve given by the implicit equation, then extract the data using  plottools:-data  command and finally do the animation using these data:

This animation is quite slow, since it contains many (about 500) points necessary for a good quality graph. By default, 10 frames are played back in one second. If you want faster playback, then open the file, click on the animation plot and in the animation panel replace the number 10 with the number 30, for example. Then the animation will play 3 times faster.

Download animcurve.mw

Try this way:

restart;
alias(X=inttrans:-laplace);

Example of use:

x:=t^2+sin(t): 
inttrans:-laplace(x, s, t);
X(x,s,t);
 

I fixed the error (which Carl indicated), removed one condition (Maple writes that one condition is superfluous), and also reduced  T (if  T>0.2  then Maple returns an error). Now the code is working:

restart;
hBar:= 1:m:= 1:Fu:= 0.2:Fv:= 0.1:
pdeu := diff(u(x,t),t)+u(x,t)/m*(diff(u(x,t),x)) = Fu;
pdev := diff(v(x, then Maplet),t)+u(x,t)/m*(diff(v(x,t),x))-hBar*(diff(u(x,t),x$2))/(2*m)+v(x,t)*(diff(u(x,t),x))/m = Fv;
ICu:={u(x,0) = 0.1*sin(2*Pi*x)};
ICv:={v(x,0) = 0.2*sin(Pi*x)};
IC := ICu union ICv;
BCu := {u(0,t) = 0.5*(1-cos(2*Pi*t))};
BCv := {v(0,t) = 0.5*sin(2*Pi*t),v(1,t)=-0.5*sin(2*Pi*t)};
BC := BCu union BCv;
pdu := pdsolve({pdeu,pdev},{IC[],BC[]},numeric, time = t,range = 0..1,spacestep = 1/66,timestep = .1);
T := 0.2; p1 := pdu:-plot3d(u,t=0..T,numpoints = 2000,x=0.0..2, shading = zhue,orientation=[-146,54,0], title = print("Figure 1",u(x, t), numeric));

Use  Physics:-diff  command for this:

restart;
alias(q = q(x, t), p = p(x, t));                             
H := lambda*p*q+conjugate(lambda)*conjugate(p)*conjugate(q)+(1/2*(p^2+conjugate(q)^2))*(conjugate(p)^2+q^2);                                   
Physics:-diff(H, q);

If you use the kernel  diff command, then the second argument in it should be a name and not a function, for example  should be  diff(sin(x^2), x)  rather than  diff(sin(x^2), x^2)

In my opinion, a more expressive result is obtained by a significantly simpler and shorter code. Why to color the lines if the surface is painted and why to rename the imaginary unit?

restart:
f:= (x,y)-> Re(sqrt(x+I*y)):
A:=plot3d([x, y, f(x,y)],x= -1..1, y= -1..1,grid= [25$2], color=f(x,y), labels= ['x', 'y',  Re(sqrt(x+I*y))], labelfont=[times,14]):
B:=plottools:-reflect(A,[[0,0,0],[1,0,0],[0,1,0]]):
 plots:-display(A,B, axes=frame, orientation=[40,80], lightmodel=light1);

   
 Such coloring emphasizes the mutual opposite of these surfaces to regard to the plane z=0.      

Here is a simple solution without DEtools package. Since you did not specify specific initial conditions, I took arbitrary  x1(0) = x10, x2(0) = x20

restart;
x:=t-><x1(t),x2(t)>;
A:=<1,-2; 4,-5>;
b:=<3,7>;
eq:=diff(x(t),t)=A.x(t)+b;
S:=Equate(op(eq));
ics:=<x10,x20>;
seq(x(0)=~ics);
Sol:=dsolve({S[],seq(x(0)=~ics)}, {seq(x(t))});

Edit. The following extra code converts the Sol (see the code above) into matrix/vector form:

Y1:=coeff(rhs(Sol[1]),exp(-t)); Y2:=coeff(rhs(Sol[1]),exp(-3*t));
Sol1:=simplify(Sol,{Y1=X1,Y2=X2});
Sol2:=(rhs=lhs)~(Sol1);
LinearAlgebra:-GenerateMatrix(Sol2,[X1,X2]);
A1,B:=%[1], x(t)-%[2];
x(t)=A1%.<Y1,Y2>+B;  # The final result

is(convert(Equate(op(value(%))),set)=Sol);  # Check

The final result (although it’s certainly easier to get it out of  Sol  manually):

DE.mw

de:=diff(y(x),x,x) - y(x) = -4*sin(x)^3 + 9*sin(x);
dsolve({de, y(0)=y(2*Pi)}, y(x));  

We see that the solution depends on one arbitrary constant  _C1 .

Here is another way to solve the original problem. We do not add a specific extra condition, but  remove the condition  y(L)=0 , leaving one condition  y(0)=0 . We get a general solution with one arbitrary constant  _С1  (then we re-designate it  as  C ). Then we impose the condition  y(L)=0  and thus get the general solution to the original problem. As an example, we find the first particular solution obtained by Rouben Rostamian, imposing additional conditions  n=10, D(y)(0)=1 :

restart; 	
de := diff(y(x),x,x) + a*y(x) = 0;
Sol:=dsolve({de,y(0)=0}, y(x));
y:=unapply(eval(eval(y(x),Sol),_C1=C), x);
a:=solve(y(L)=0,a, allsolutions);
about(_Z1);
a:=eval(a,_Z1=n);
y(x);
y:=unapply(simplify(%), x) assuming L>0,n::posint; # This is the general solution to de
eval(de,x=0), eval(de,x=L);  # Check

n:=10: # Example of  Rouben Rostamian
C0:=solve(D(y)(0)=1, C);
eval(y(x),C=C0);

 
df.mw

If you mean polynomial interpolation, then use  CurveFitting:-PolynomialInterpolation  command. Before this, replace each curly brace with a square one, because in Maple braces creates a set, and we need lists.

data := [[13, -2], [12, -1], [11, 0.0], [10, 1], [9, 2], [8, 3], [7, 
   4], [6, 5], [5, 6], [4, 7], [3, 8], [2, 9], [1, 10]]:
CurveFitting:-PolynomialInterpolation(data, x);

Output:                         11.-1.*x

We see that all these points lie on one straight line  y = 11 - x . Visualization confirms the result:

with(plots):
display(pointplot(data), plot(11-x, x=1..13), view=[0..13,-2..10]);

It is impossible to understand from OP's code whether he is trying to calculate the gradient of the function of two variables (this is a vector) or the mixed partial derivative with respect to variables  x1  and   x2  (this is a scalar). In both cases there is no need to call any packages (it will be shorter and faster).

restart;
x:=<x1,x2>:
W:=<<w11, w12>|<w21, w22>>;
Cal:=(W.x)^%T.(W.x);
<seq(diff(Cal,t),t=x)>;  # Gradient
diff(Cal,x1,x2); # Partial derivative