Procedure  PartialSum calculates partial sums of the series   a1+a2+...+an+...

Formal parameters:  N  is the number of  terms, an is an expression for nth term of the series.

PartialSum:=proc(N, an)
local psum, f, k;

psum:=0;
f:=unapply(an, n);
for k from 1 to N do
psum:=psum+f(k);
od;
psum;

end proc:

Example of use:

PartialSum(10, 1/n);

                               7381/2520

Addition:  Of course,  PartialSum  procedure is simply a correction of  OP's attempt. In actual practice, you can simply write 

 add(an, n=1..N);  # for numeric  N

sum(an, n=1..N);  # for symbolic  N

Examples:

add(1/2^n, n=1..10);

sum(1/2^n, n=1..N);

eval(%, N=10);

Markiyan Hirnyk  In fact, you've got the original parameterization  x=9z^2+3z+367, y=6z^2 +z+244  and nothing new.

If  x, y  are integers, we will prove that there must be  z  is an integer (this is obvious)  or  z  is  an integer+1/3

Let   x, y  are integers. It is easy to get  3*y-x=9*z^2+365 . It follows that 3*z must be an integer. Let  3*z=n ,  z=n/3 , n is integer.  It is clear that  x  be  an integer for all  z=n/3 , n is integer. It remains to find conditions on n, y to be an integer.  To do this, we have to solve in integers  the equation  6*(n/3)^2+n/3=m :

isolve(6*(n/3)^2+n/3=m);

             {m = 6*_Z1^2-_Z1, n = -3*_Z1}, {m = 6*_Z1^2-5*_Z1+1, n = 1-3*_Z1}

From the second solution we have  z=n/3=-_Z1+1/3

minimize(sin(x+y), x=-1..1, y=-1..1);

                                    -1

With  location  option:

minimize(sin(x+y), x=-1..1, y=-1..1, location);

               -1, {[{x = 1-1/2*Pi, y = -1}, -1]}

Your cone depends on two parameters  R  and  Theta, so we can set a procedure with these parameters, and use it by setting specific parameters.

Two ways:

Cone := (R, Theta) -> plot3d(eval([r*cos(phi)*sin(theta), r*sin(phi)*sin(theta), r*cos(theta)], theta = Theta), r = 0 .. R, phi = 0 .. 2*Pi, scaling = constrained) :

# or

Cone1 := (R, Theta) ->plot3d(eval([r, phi, theta], theta = Theta), r = 0 .. R, phi = 0 .. 2*Pi, coords = spherical, scaling = constrained) :

Examples of use:

Cone(1, (1/6)*Pi);

Cone1(2, (1/4)*Pi); 

Addition: The third way - to use  plottools[cone]  command:

Cone2 := (R, Theta)->plots[display](plottools[cone]([0, 0, 0], R, R*cot(Theta)), axes = box, scaling = constrained):

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]