Second call  f:=x->a*x  is quivalent to a procedure  but any procedure works (as you want) only if  its parameters are specified:

a:=3:

f:=proc(a,x)

a*x;

end;

f(a,x);

Your code works correctly, I have checked.
Here is a simplified version of your code (no need in any packages), which does the same:

restart;

s:={(1,3)=a,(1,4)=-a,(1,5)=-d,(1,6)=d,(2,4)=-d,(2,5)=d,(2,6)=-d,(2,7)=d,(3,5)=a,(3,6)=d,(3,7)=-d,(4,6)=-d,(4,7)=-a,(5,7)=a}:

A:=Matrix(7, s, shape=antisymmetric);

for a in [-1,1] do  for d in [-1,1] do  print('a'=a,'d'=d, A)  od; od;

You can just copy this code and run.

Addition. I understood the reason for error in the repetition of your code . You use variables  a  and  d , which are already assigned specific values. Therefore, if you want to avoid errors in repeating of the code, at the beginning of your code insert the line  a:='a': d:='d':  or  restart;  as in my version.

Here is another variant, which does not need   a:='a': d:='d':  or  restart;  under repetition:

s:={(1,3)=a,(1,4)=-a,(1,5)=-d,(1,6)=d,(2,4)=-d,(2,5)=d,(2,6)=-d,(2,7)=d,(3,5)=a,(3,6)=d,(3,7)=-d,(4,6)=-d,(4,7)=-a,(5,7)=a}:

A:=Matrix(7, s, shape=antisymmetric);

for i in [-1,1] do  for j in [-1,1] do  print(a=i, d=j, eval(A, {a=i, d=j}))  od; od;

Your system has no real solutions. For solving I assumed that angles  alpha2, .. ,alpha5  are specified in degrees (not in radians) and added 4 extra equations.  The system has become a polynomial system. Maple can find all the numerical solutions of such a system, but it did not find anything:

restart;

eq1 := (cos(beta2)-1)*w11-sin(beta2)*w12+(cos(alpha2)-1)*z11-sin(alpha2)*z12-cos(delta2) = 0:

eq2 := (cos(beta2)-1)*w12+sin(beta2)*w11+(cos(alpha2)-1)*z12+sin(alpha2)*z11-2*sin(delta2) = 0:

eq3 := (cos(beta3)-1)*w11-sin(beta3)*w12+(cos(alpha3)-1)*z11-sin(alpha3)*z12-3*cos(delta3) = 0:

eq4 := (cos(beta3)-1)*w12+sin(beta3)*w11+(cos(alpha3)-1)*z12+sin(alpha3)*z11-4*sin(delta3) = 0:

eq5 := (cos(beta4)-1)*w11-sin(beta4)*w12+(cos(alpha4)-1)*z11-sin(alpha4)*z12-5*cos(delta4) = 0:

eq6 := (cos(beta4)-1)*w12+sin(beta4)*w11+(cos(alpha4)-1)*z12+sin(alpha4)*z11-6*sin(delta4) = 0:

eq7 := (cos(beta5)-1)*w11-sin(beta5)*w12+(cos(alpha5)-1)*z11-sin(alpha5)*z12-7*cos(delta5) = 0:

eq8 := (cos(beta5)-1)*w12+sin(beta5)*w11+(cos(alpha5)-1)*z12+sin(alpha5)*z11-8*sin(delta5) = 0:

alpha2 := -20*Pi/180: alpha3 := -45*Pi/180: alpha4 := -75*Pi/180: alpha5 := -90*Pi/180: delta2 := 15.5: delta3 := -15.9829: delta4 := -13.6018: delta5 := -16.7388: P21 = .5217: P31 = 1.3421: P41 = 2.3116: P51 = 3.1780:

Sys:={seq(eq||i, i=1..8)}:

Sys1:=eval(Sys,{cos(beta2)=cb2,sin(beta2)=sb2,cos(beta3)=cb3,sin(beta3)=sb3,cos(beta4)=cb4,sin(beta4)=sb4,cos(beta5)=cb5,sin(beta5)=sb5}) union {cb2^2+sb2^2=1,cb3^2+sb3^2=1,cb4^2+sb4^2=1,cb5^2+sb5^2=1};

RealDomain[solve](Sys1);

restart;

Digits:=20:

f:=x->(1+1/x)^x;

seq((evalf(1+1/10^n))^(10^n), n=1..10);

evalf(exp(1));  # Exact value for comparison

In  diff(f, x) the first argument  f  should be an expression not a function. Carl gave the best solution for your problem. If you still want to use  diff  command then  unapply  command helps you:

f:=x->x^2;

derivative:=unapply(diff(f(x),x), x);

derivative(2);

You can not write  

derivative:=x->diff(f(x), x);

derivative(2); 

because in this case Maple substitutes  x=2  before the derivative is calculated and an error occurs.

Here is a solution using a finite difference method. We use a simple formulas:  D(y)(x[i])=(y[i+1]-y[i-1])/2/h ,(D@@2)(y)(x[i])=(y[i+1]-2*y[i]+y[i-1])/h^2  and Simpson formula for calculation of the integral. This approach requires the boundary conditions. In the first example we specify the boundary conditions as the exact solution  y=exp(x)  has:

restart;

Eq:=x^2*diff(y(x), x, x)+50*x*diff(y(x), x)-35*y(x) = (1-exp(x+1))/(x+1)+(x^2+50*x-35)*exp(x)+Int(exp(x*t)*y(t), t = 0 .. 1):

n:=10:

h:=1/2/n:

Simpson:=L->h/3*(L[1]+2*add(L[3+2*k],k=0..n-2)+4*add(L[2+2*k],k=0..n-2)+L[-1]):

y[0]:=1: y[2*n]:=exp(1):

L:=[y[0],seq(y[i], i=1..2*n-1),y[2*n]]:

for i from 1 to 2*n-1 do

eq[i]:=eval(Eq,[diff(y(x),x)=(y[i+1]-y[i-1])/2/h, diff(y(x),x,x)=(y[i+1]-2*y[i]+y[i-1])/h^2,y(x)=y[i], Int(exp(x*t)*y(t), t = 0 .. 1)=evalf(Simpson([seq(exp(i/2/n*k/2/n)*L[k], k=1..2*n)])), x=i/2/n]);

od:

Sys:=convert(eq,set):

R:=solve(Sys);

plot([exp(x), [[0,y[0]],seq([i/2/n,rhs(R[i])],i=1..2*n-1),[1,y[2*n]]]], x=0..1, color=[blue,red], thickness=0, view=[0..1,0..y[2*n]]);  # Visual comparison of exact and approximate solutions (the plots merge)

In the second example we specify the other boundary conditions and get a different solution:

restart:

Eq:=x^2*diff(y(x), x, x)+50*x*diff(y(x), x)-35*y(x) = (1-exp(x+1))/(x+1)+(x^2+50*x-35)*exp(x)+Int(exp(x*t)*y(t), t = 0 .. 1):

n:=10:

h:=1/2/n:

Simpson:=L->h/3*(L[1]+2*add(L[3+2*k],k=0..n-2)+4*add(L[2+2*k],k=0..n-2)+L[-1]):

y[0]:=1: y[2*n]:=1.5:

L:=[y[0],seq(y[i], i=1..2*n-1),y[2*n]];

for i from 1 to 2*n-1 do

eq[i]:=eval(Eq,[diff(y(x),x)=(y[i+1]-y[i-1])/2/h, diff(y(x),x,x)=(y[i+1]-2*y[i]+y[i-1])/h^2,y(x)=y[i], Int(exp(x*t)*y(t), t = 0 .. 1)=evalf(Simpson([seq(exp(i/2/n*k/2/n)*L[k], k=1..2*n)])), x=i/2/n]);

od:

Sys:=convert(eq,set):

R:=solve(Sys);

plot([[0,y[0]],seq([i/2/n,rhs(R[i])],i=1..2*n-1),[1,y[2*n]]],x=0..1,color=red, thickness=0, view=[0..1,0..y[2*n]]);

Can we conclude that there is more than one solution of the original problem?

See help on  Explore  command. To plot vectors you can use  plots[arrow]  command. The imaginary unit is coded (by default) in Maple as  I  not  j .

Of course the second calculation  evalf(gg(2))=-0.3333333333  is incorrect. Probably Maple uses formal geometric series summation  b+b*q+b*q^2+..+b*q^n+..=b/(1-q)  which is correct only if  |q|<1 .

Indeed for the series  gg(2)=Sum(2^jj*2^jj, jj = 0 .. infinity)   we have  b=1 ,  q=4    and   1/(1-4)=-0.3333333333

You can search for a numerical solution using  numeric  option. To do this, all the parameters and initial conditions must be specified. With the resulting procedure, you can find solutions in some points, build plots, etc.

sys := {diff(x(t), t) = -k*x(t)*sqrt(x(t)^2+y(t)^2)/m, diff(y(t), t) = -k*y(t)*sqrt(x(t)^2+y(t)^2)/m - g}:

sol := dsolve({eval(op(sys), {g = 9.8, k = 1, m = 2}), x(0) = 0, y(0) = 0}, numeric);

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 2, color = [red, blue], thickness = 3);

If you mean the 50 roots of this equation, then here are the first 50 positive roots:

seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49);

 Addition. You can easily plot these points with the following code (for clarity, plotted the first 10 roots in the form of small red circles):

L:=[seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49)]:

plots[display](plot(0.5*lambda*tan(lambda)-1, lambda=0..Pi/2+9*Pi+0.01, -10..10, color=blue, discont), plot(map(t->[t,0], L[1..10]), style=point, symbol=solidcircle, symbolsize=12, color=red, scaling=constrained));

Complex function of a complex argument may be regarded as a mapping from R^2 into R^2. On the left - the space of arguments,  on the right - the space of their images:

R:=evalc(1/(1-(x+I*y)));

R1:=op(1,R),coeff(R,I);

f:=unapply([R1],x,y);

F:=plottools[transform](f):

A:=plot([seq([convert(<cos(Pi/4*i),-sin(Pi/4*i);sin(Pi/4*i),cos(Pi/4*i)>.<t,0>,list)[], t=-0.8..0.8], i=0..7), seq([r*cos(t),r*sin(t),t=0..2*Pi], r=0.2..0.8,0.2)], color=[red,gold,brown,yellow,blue,green,violet,cyan], thickness=3):

plots[display](<A | F(A)>, scaling=constrained);

 Edited.

If we go to polar coordinates, Maple can easily find the solution in explicit form, even simpler than Mathematica does:

restart;

T:={x(t)=r(t)*cos(phi(t)), y(t)=r(t)*sin(phi(t))}:  # The transformation to polar coordinates

sys:={diff(x(t),t) = -k/m*x(t)*sqrt(x(t)^2+y(t)^2), diff(y(t),t) = -k/m*y(t)*sqrt(x(t)^2+y(t)^2)}:  # The original system

subs(csgn(r(t))=1, simplify(eval(sys,T)));  # The transformation of original system to polar coordinates

sys1:=solve(%,{diff(r(t),t),diff(phi(t),t)});  # The original system in polar coordinates in standard form

dsolve(sys1, {r(t),phi(t)});  # The general solution in polar coordinates

dsolve(sys1 union {r(0)=1,phi(0)=Pi/3},{r(t),phi(t)});  # The solutin with some initial conditions

eval(T, %);  # The same solution in Cartesian coordinates

In fact, the point moves along a straight line toward the origin.

ConvIntoIneq  procedure solves the problem:

ConvIntoIneq:=proc(Rel)

local t, R, x, y;

t:=indets(Rel)[];

R:=solve(Rel);

x, y:=op(1,R), op(2,R);

if x=-infinity then if type(y,realcons) then return cat(t, "<=", y) else

cat(t, "<", op(1,y)) fi else

if type(x,realcons) then return cat(t, ">=", x) else

cat(t, ">", op(1,x)) fi; fi;

end proc:

Examples of use:

ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100) ;

 Addition. Unfortunately, if we want to get a "nice" output, the name of variable can appear on the right side of the inequality:

map(parse, [ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100)])[] ;

I added 2 options into your code: axes=normal and orientation=[15,50]  and used user lighting :

 Addition. The plot was saved in Paint as png - file.

Another variant. Additionally I added  color=khaki ,  style=surface , numpoints=20000  and changed  user lighting  :

plot([[0,0.36],[5,0.38],[10,0.3],[15,0.18],[20,0.96],[25,0.18],[30,0.16],[35,0.96]], style=line, thickness=5, color=red, view=[0..35,0..1], axis[1] = [gridlines = [7, thickness = 1, subticks = false, color = grey]], axis[2] = [gridlines = [6, thickness = 1, subticks = false, color = grey]], tickmarks=[[0=1,5=5,10=10,15=50,20=250,25=1000,30=10000,35=50000], spacing(0.2,0)], size=[650,350], font=[COURIER,ROMAN,14]);

 Edited.