FWIW, here's my rational roots procedure:

RationalRoots:= proc(f::And(polynom(integer), satisfies(f-> nops(indets(f,name))=1)))
uses NT= numtheory;
local x:= indets(f,name)[], p, q;
     select(
          r-> eval(f, x= r) = 0,      
          {0, seq(seq([p/q, -p/q][], p= NT:-divisors(tcoeff(f))), q= NT:-divisors(lcoeff(f)))}
     )
end proc:

f:= x-> piecewise(x<0, x^2-2*x+1, x>0, 1-x, undefined);

Why do you say "Piecewise function in Maple are continuous function"?

FWIW, here's my somewhat minimalist Jacobi method procedure. I simply used the Wikipedia article "Jacobi method" and modified the algorithm to use Maple's vector-mode operations. I use two convergence criteria: the magnitudes of both the error and the residual need to be less than epsilon.

Jacobi:= proc(A::Matrix, b::Vector, x0::Vector, epsilon::positive, m::posint)  
uses LA= LinearAlgebra;  
local
     k,
     x:= Vector(x0),
     p:= Vector(op(1,x)),
     D:= Matrix(A, shape= diagonal),
     R:= A-D,
     B:= Vector(b)
;  
     D:= D^(-1);  B:= D.B;  R:= D.R;
     for k to m while LA:-Norm(A.x-b) > epsilon or LA:-Norm(x-p) > epsilon do
          p:= x;
          x:= B-R.p      
     end do;
     if k > m then WARNING("Max iterations reached.") end if;    
     x  
end proc:

Example of use:

n:= 3:
A:= LinearAlgebra:-RandomMatrix(n):
#Force A to be strictly diagonally dominant:
A:= A^+.A:
A:= Matrix(op(1,A), (i,j)-> `if`(i=j, n, 1)*A[i,j]):
b:= LinearAlgebra:-RandomVector(n):
x:= Jacobi(evalf(A), evalf(b), Vector(n, fill= 1.), 1e-6, 99);

A.x-b;

Here's how to create a random real symmetric float matrix with determinant 2:

macro(LA= LinearAlgebra):
n:= 9:
A:= LA:-RandomMatrix(n, datatype= float[8], shape= symmetric):
det:= LA:-Determinant(A):
A:= signum(det)*(2/abs(det))^(1/n)*A;

The major problem with your procedure is that f isn't defined. You've made f local; it should be a parameter.

A second problem is that a parameter can't be indexed. So change p[0] to p0.

A third problem is that a parameter can't (usually) be assigned to. You attempt to assign to p[0]. If you want to update it, then you need to copy it to a local variable.

omega:= (y::function(name))-> y*diff(y, op(y))/op(y):

omega(y(x));

(In the future, please post small worksheets inline as well as attaching them.)

You had many errors, mostly syntax. I don't have time right now to detail all the syntax errors, but I corrected them all. You should carefully compare what I wrote with what you wrote. Many of the changes I made were simply my personal preference syntax rather than correction of errors.

The glaring logic error that you made is

seq(EnergySpacingGUE(1,1), i= 1..10)

The value of n must be smaller than N, because if there are N values, then there are only N-1 spacings. I changed this to

seq(EnergySpacingGUE(i,1), i= 2..11)

but I'm only guessing your intention!

I guess that you realize that a histogram of only ten data points is ridiculous, and that this was only a proof-of-concept case.

restart:
UseHardwareFloats:= true:
GUE:= proc(N::posint)
uses LA= LinearAlgebra, RT= RandomTools;
local
     R:= RT:-Generate(distribution(Normal(0,1)), makeproc),
     G:= 'Matrix((N,N), R, datatype= float[8])' $ 2,
     C:= G[1] + I*~G[2],
     H:= (C + C^*)/~2
;
     convert(LA:-Eigenvalues(H)/~sqrt(4.*N), list)[]
end proc:

EnergySpacingGUE:= proc(N::posint, n::posint)
local
     Eigenvsorted:= sort(Re~([GUE(N)])),
     Spacing:= abs~(Eigenvsorted[..-2] -~ Eigenvsorted[2..])
;
     Spacing[n]*N/`+`(Spacing[])
end proc:
      
data:= [seq(EnergySpacingGUE(i,1), i= 2..10)];

plots:-display([
     plot(32/Pi^2*x*exp(-Pi/4*x^2), x= 0..2),
     Statistics:-Histogram(data, binwidth= 0.05, range= 0..2)
]);

The linestyle option will accept an integer in 0..7 or a list of such integers.

In the vast majority of cases, when you want to sort a list of complex numbers, you want to sort them by absolute value (aka magnitude). Since Maple 18, option key to sort provides an easier, more-intuitive, and more-efficient way to do this than using a comparison procedure as in Kitonum's Answer. Let C be a list or vector of complex numbers. Then do

sort(C, key= evalf@abs);

If the list C is already in floating-point form, that can be shortened to

sort(C, key= abs);

If you want a lexicographic sort, as in Preben's Answer, this can be done using key by

sort(C, key= evalf);

If C is already floating-point, then that can be simply

sort(C);

That's because lexicographic order is the default order for numeric complex constants.

In the place where you've entered the piecewise expression, in the second row, second column, you've entered two inequalities separated by a comma. This isn't valid, yet the despicable 2-D input doesn't catch the error. Instead, it interprets it as a three-row piecewise with nonsense in the second row, first column. Instead of a comma between conditions, you may use the keyword and.

But, you haven't specifed any value for when |x| >= 10; you might as well let the value be 0. It that case, there's no need to even mention the 10. So, the entire piecewise can be entered as

piecewise(abs(x) < 5, -1, 0)

That's the 1-D input. If you want the 2-D input, I'll leave it to you or someone else to translate it.

You have two errors. The first is that you need to call GramSchmidt with option normalized. The second is that the line ord:= Basis(...) simply overwrites the results of GramSchmidt. Here's an example of constructing an orthonormal basis starting with the eigenvectors.

restart:
macro(LA= LinearAlgebra):
n:= 5:
A:= LA:-RandomMatrix(n, datatype= float[8]):

#The 2nd return value of Eigenvectors is a matrix whose columns are the eigenvectors.
E:= LA:-Eigenvectors(A)[2]:

#GramSchmidt expects a list of vectors as input, so you need to split E into columns.
#Option 'normalized' is required if you want normalized vectors as output.
G:= LA:-GramSchmidt(['E[..,k]' $ k= 1..n], 'normalized'):

#Form the orthonormalized vectors into a basis matrix.
B:= `<|>`(seq(G[k], k= 1..n)):

#B^* is the conjugate transpose of B. This line is simply verifying that
#B is orthonormal (upto floating-point rounding error).
simplify(fnormal~(B.B^*), zero);

The following nonsensical subexpression occurs seven times in your formula for G:

[7.6613, 8.0102, 8.1790]*[1]

What is it supposed to mean?

The simplest solution is to simply not use gamma as a variable. The next simplest solution is to begin your code with the line

local gamma;

This latter solution will only work in relatively recent versions of Maple.

If either of these solutions is unacceptable, there are still other alternatives.

Zero-point and negative-point users have never appeared on the Users' list in the three years that I've been here. That's why I give zero-point users a vote up every chance I get (except for spammers, of course). Since only Questions, Posts, and Answers can be voted up---and Mcloone has never made any of those---I never had a chance to give him a vote up. Since only Questions, Posts, and Answers get indexed, you won't find him by searching for words from his Comments and Replies.

The following worksheet shows that there are two solutions for N which produce equally minimal residuals and which both lead to n = 14.

Download SolveForN.mw