I corrected a few errors in your file, wrote about Mac Dude. Maple found 4 solutions of the first 9 equations. But none of them does not satisfy the additional condition  C[1] . 

For details, see the attached file:

NL_new.mw

Addition: If a system has аn expression  A , instead of an equation, Maple treats this as the equation  A=0

It seems that this is not possible to do interactively only  programmatically.

An example:

plot(x^2, x = 0 .. 2, style = point, numpoints = 20);

To create a matrix  A  does not necessarily advance to create its entries. It is easier and shorter to specify the function (i,j)->A[i,j] that assigns a certain number  A[i,j]  to ith row and jth column :

psi := (1/4)*(1-Zeta)^2*(2+Zeta), (1/4)*(1-Zeta)^2*(Zeta+1), (1/4)*(Zeta+1)^2*(2+Zeta), (1/4)*(Zeta+1)^2*(Zeta-1);

A := Matrix(2,  (i, j) -> 2*psi[i]*psi[j]);

The procedure  Root  solves the problem:

Root:=a->fsolve(-x^3+a*x^2-ln(x)=0, x=0..infinity);

Examples of use:

Root(0.5),  Root(1),  Root(10);

                       0.8130586971,  1.,  9.976890600

Now the code works:

with(plots):

 m := 0.46e-1; d := 0.42e-1; v := 60; alpha0 := 12; g := 9.81; pa := 1.205; cd := .2; n := 100; omega := 2*Pi*(1/60);                           

 p := 6*m/(Pi*d^3);

 k1 := (3/4)*cd*pa/(d*p); k2 := (3/8)*omega*pa/p;                                       

 gl1 := vx(t) = diff(x(t), t);                   

 gl2 := vy(t) = diff(y(t), t);                             

 gl3 := diff(vx(t), t) = -k1*vx(t)*(vx(t)^2+vy(t)^2)^(1/2)-k2*vy(t);

 gl4 := diff(vy(t), t) = -g-k1*vy(t)*(vx(t)^2+vy(t)^2)^(1/2)+k2*vy(t);

 init1 := x(0) = 0;

 init2 := y(0) = 0;

 init3 := vx(0) = v*cos((1/15)*Pi);                              

 init4 := vy(0) = v*sin((1/15)*Pi);

sol:=dsolve({gl1, gl2, gl3, gl4, init1, init2, init3, init4}, {vx(t), vy(t), x(t), y(t)}, type = numeric);

sol(0.1);                        

odeplot(sol, [[t,x(t)], [t, y(t)]], t = 0 .. 1, color=[red,blue]);

In the example below, to keep the order of the elements I replaced the braces on brackets, also I added a few more negative numbers:

a:=[10, 11, -13, -9, 20, 74, -10]:

map(t->`if`(t<0,0, t), a);

                          [10, 11, 0, 0, 20, 74, 0]

Addition.  Easy modification of your code also works

a:=[10, 11, -13, -9, 20, 74, -10]:

for k from 1 to nops(a)  do

  if a[k]<0 then a:=subs(a[k]=0, a) end if:

end do:

a;

                            [10, 11, 0, 0, 20, 74, 0]

The procedure  SM  solves the problem.

SM := n->Matrix(n+1, (i,j)-> `if`(j <= i, (-1)^(i+j-2)*(i+j-2)!/((i-j)!*(j-1)!), 0)):

Example of use:

SM(4);

It seems that this is exactly what you want.

Edited.

I have seen on your picture what you are trying to solve the equation  f(a, x)=x  for  x  with the parameter  a . Similar equations can be solved numerically for each parameter value  a , eg using  fsolve  command. We first prove that this equation has only the unique positive root for any real values of  a . To do this, we are expressing  the parameter  a  from the equation, going to the equivalent equation  a=g(x) . It is easy to prove that the function  g(x)=(x^3+ln(x))/x^2  takes any real values and strictly increasing for x>0 . This implies the existence and uniqueness of the root of the equation   f(a, x)=x  for any  a .

f:=(a,x)->exp(x^2*(a-x)):

Eq:=f(a,x)=x:

a=solve(Eq, a);  # The equivalent equation  to  the equation  f(a, x)=x

g:=unapply(rhs(%), x);  # The function defined by the right side of previous equation

plot(g(x), x=0..10, -10..10, thickness=2, labels=[x,'g'(x)]);

Root:=a->fsolve(exp(x^2*(a-x))=x, x=0..infinity):  # Solution of the equation  f(a, x)=x  for any  a

Examples of use:

Root(-2),  Root(0),  Root(2),  Root(5);

                      0.5142801644,  0.7047094903,  1.819187100,  4.934442571

Replace  style = line  in your code with  linestyle=solid  or  style=surface

Do you mean step by step solution from the definition of the derivative?

An example (finding the derivative of the function  x->x^2 ):

f:=x->x^2:

Student[Calculus1][ShowSolution](Limit((f(x+dx)-f(x))/dx, dx=0));

 Addition.  Unfortunately, in more complicated situations, for example for piecewise functions, it does not work:

f:=x->piecewise(x=0, 0, x^2*sin(1/x)):

Student[Calculus1][ShowSolution](Limit((f(x+dx)-f(x))/dx, dx=0));

Nothing was returned.

Calculation by  diff  command is working properly:

 diff(f(x), x);

The procedure  P  generates a list  [x,y]  of 2 integers in the range  r  such that  x<y

P:=proc(r)  # r is a range

local x, y;

uses RandomTools;

x:=Generate(integer(range=lhs(r)..rhs(r)-1));

y:=Generate(integer(range=x+1..rhs(r)));

[x,y];

end proc:

Example of use:

seq(P(1..7), i=1..10);

                        [2, 7], [5, 6], [1, 6], [4, 6], [6, 7], [4, 5], [3, 6], [1, 6], [3, 7], [1, 3]

There is a unique solution to this problem. In the code below the procedure  NumbersGame  from  here  was used:

L:=[$1..9]:

[seq(op(combinat[choose](8,i)), i=1..8)]:

M:=map(t->[parse(cat(op(L[1..t[1]]))),seq(parse(cat(op(L[t[i]+1..t[i+1]]))),i=1..nops(t)-1),parse(cat(op(L[t[-1]+1..-1])))],%):

Sol:=[]:

for k in M do

S:=NumbersGame(2016, k, ["+","-","*"]):

if S<>`No solutions` then Sol:=[op(Sol),S] fi;

od:

op(Sol);

                                     1 * 2 * 34 * 5 * 6 - 7 - 8 - 9 = 2016

Addition. Also there is a unique solution to this problem for numbers  1 .. 8

                                         1 * 2 + 345 * 6 - 7 * 8 = 2016

See help on  Statistics[LinearFit]  command.

expr:=(-(1/2)*ib-(1/2)*ia+ic)*vc+(-(1/2)*ia-(1/2)*ic+ib)*vb+(-(1/2)*ic-(1/2)*ib+ia)*va;

A:=algsubs(ia+ib+ic=0, expr);

subsop([3,2]=simplify(op([3,2], A), {ia+ib+ic=0}), A);

The procedure  SM  solves your problem for both numeric and symbolic coefficients. Formal parameters: ex  - your expression,  U - the list of variables that determine the rows of the matrix, A - the list of variables that determine the columns of the matrix:

SM := proc (ex, U, A)

Matrix(map(t->[seq(coeff(t, A[i]), i = 1 .. nops(A))], map2(coeff, ex, U)));

end proc:

Examples of use:

ex1 := (3*a[1]+5*a[2])*U[1]+(-6*a[1]+2*a[2])*U[2];

ex2 := (A1*a[1]+B1*a[2])*U[1]+(A2*a[1]+B2*a[2])*U[2];

SM(ex1, [U[1], U[2]], [a[1], a[2]]);

SM(ex2, [U[1], U[2]], [a[1], a[2]]);