display(fig1, plot(eval(Model, MinSol[2]), x = 0 .. 1.5, color = red, thickness = 0));

You can generate a random function of a certain class of functions depending on a finite number of parameters.

Example: a random function of the class  a*sin(b*x+c)  in which  [a,b,c]  is a random list of floats with the uniform distribution in the range  0..10 :

RandomTools[Generate](list(float(range=0...10, method=uniform), 3));

eval(a*sin(b*x+c), [a,b,c]=~%);

Should be  D(x)(0)  instead of  D(x) (0)

A space between  D(x)  and  (0)  Maple interprets as a multiplication.

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 .