For any real numbers b>0, c>0  the expression sqrt(sqrt(9)*sqrt((1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c))+3+(3*b^2+6*b+3)*c^2+(8*b-6)*c)  will be strictly positive . Here is the simple solution (in fact by hand):

u:=sqrt(sqrt(9)*sqrt((1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c))+3+(3*b^2+6*b+3)*c^2+(8*b-6)*c);  # The original expression

v:=normal(u^2);  # The simplified expression under the square root

v1:=3*b^2*c^2+6*b*c^2+3^(1/2)*((3*b^2*c^2+6*b*c^2+10*b*c+3*(c-1)^2)*(b^2*c^2+2*b*c^2+2*b*c+(c-1)^2))^(1/2)+8*c*b+3*(c-1)^2;  # In  v  we substitute c^2-2*c+1=(c-1)^2  # All terms are strictly positive if  c<>1, c>0, b>0

eval(v1, c=1);  # for c=1

If the parameters  beta, gamma, R4 and C2  are not necessarily positive, vv's solution does not work. Assume  command helps in this case. Also I wrote  local gamma;  because  gamma is protected constant (Euler's constant).

restart;

local gamma;

assume(gamma>0,beta<0,R4>0,C2<0):

f:=gamma*sqrt(4)*sqrt(omega^2*C2^2*R4^2/(C2^4*R4^4*beta^2*gamma^2*omega^4+C2^2*(1+gamma^2*(beta+1)^2-2*gamma)*R4^2*omega^2+1)):

diff(f,omega):

solve(%=0, omega);

select(t->is(t>0),[%])[];

I think the reason is a low accuracy of calculations (10 digits), while the coefficients of the polynomial specified with high precision. I set  Digits:=100  and use  fsolve  command, which finds all real roots for polynomials:

Digits:=100:

P:=-12116320194738194778134937600000000*t^26+167589596741213731838990745600000000*t^24+1058345691529498270472972795904000000*t^22-4276605572538658673086219419648000000*t^20-23240154739806540070988490473472000000*t^18-5442849111209103187871341215744000000*t^16+49009931453396028716875310432256000000*t^14+74247033158233643322704589225984000000*t^12-2762178990802317464801412907008000000*t^10-25947900993773120244883450232832000000*t^8-7468990043547273070742668836864000000*t^6-567730116675454293925108383744000000*t^4+3703566799705707258760396800000000*t^2-4742330812072533924249600000000:

Sol:=fsolve(P):

seq(eval(P,t=s),s=[Sol]);  # Check

-.25965814962389295505e-49, -.39367217335124e-55, -.4842468019e-59, .1022039e-62, -.1e-68, 0., 0., -.1e-68, .1022039e-62, -.4842468019e-59, -.39367217335124e-55, -.25965814962389295505e-49

Vic, Carl gave a solution for a matrix with 2 parameters  a  and  d .  The following procedure  named  Tr  solves your problem for any number of parameters. The procedure returns the sequence of 2^N  matrices, where  N is the number of parameters :

Tr:=proc(n,s)

local A, S, N;

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

S:=indets(s);

N:=nops(S);

seq(eval(A,S=~p), p=combinat[permute]([-1$N,1$N],N));

end proc:

Examples of use.

The first one is yours:

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}:

Tr(7,s);

Tr(7,s)[1]  is the first matrix in this sequence and so on.

The second example solves the same problem if we have 3 parameters  a , d , e :

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

Tr(5,s);

This equation you can solve only numerically. Here is the least positive root:

fsolve(2*x-tan(x), x=0..Pi/2, avoid={x=0});

                                        1.165561185

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));