For an unknown reason, this 10-year-old thread has been refreshed. So, I thought that this would be a good place to point out that the Question can now be Answered with a stock solution built-in to sort:

sort(x, 'output'= 'permutation');
             [1, 3, 2, 4]

Another way of thinking of it is that we want the 8th roots of 1+I:

(M,A):= (abs, argument)(1+I):
plot(((M^(1/8)*~[cos,sin])~([seq]((2*Pi*k+A)/8, k= 0..8))), style= pointline);

Maple does not automatically recognize vectors (or matrices) with equal entries as being equal, but they will recognized if they are converted to lists. If at the end you do

K1:= [seq]~(x, x=~ K);

then they'll be converted to a set of 2-lists. If you do

K1:= `<,>`~([seq]~(x, x=~ K)); 

then they'll be converted to a set of 2-vectors.

Are you trying to find the orbit of L under mulplication by M?

Your major problem is that doSum doesn't use its parameter first.

I suspect that this is an issue related to installation and licensing. I've seen numerous similar reports where that was the issue. You have entered the differential equation correctly. 

Try going to menu Tools -> Check for Updates. 

Here's my GaussSeidel, coded for Maple 2019 or later. The Maple 2019 syntax makes it easy to see that the algorithm is simply adding the residuals until the residual is small. I post this because this procedure is competitive (both in efficiency and accuracy) with stock methods provided by LinearAlgebra:-LinearSolve. Also, I show a little-known trick for including a plot in a userinfo statement.
 

Download GaussSeidel.mw

See the Sockets package: ?Sockets

In your first line of NEUZMinus,

NEUZMinus:= proc(Unten, Oben, f,G,Liste,n)::real;

the ::real is what I call a return-value type assertion^[*1]. It needs to be terminated with a semicolon, which you have. The semicolon terminates the assertion itself; it's purpose is not to terminate the proc(...), which needs no terminator. However, your first line of Basenweschel

Basenwechsel:=proc(Dividend, m);

is incorrect because of the semicolon without a type assertion. It is only a coincidence that this particular procedure is syntactically accepable in Maple 2020; it would be incorrect if there were any declarations other than local on the next line, even in 2020.

[*1] There is no stock type named real; instead, real is a property. If you want to include a type assertion, you should make it realcons. Type assertions other than of procedure parameters are only checked if you set kernelopts(assertlevel= 2).

Your problem has nothing to do with the Task programming model. If you use lists and nops, then nops(L) is always a valid index of list L. The same is not always true of a 1-dimensional Array and rtable_num_elems if the Array starts at an index other than 1. The command lowerbound(A) is threadsafe and will return the first index of A.

The command coeffs has a long-standing flaw that it's awkward to match monomials to their coefficients. The following procedure corrects that:

Coeffs:= proc(f, V::{set,list}(name))
local C, t;
    C:= coeffs(f, V, t);
    table(sparse, [t]=~[C])
end proc
:
F:= 2*x+6*y+4*x^2+12*x*y-5*y^2:
C:= Coeffs(F, [x,y]):
C[x], C[x^2], C[x*y];
                            2, 4, 12

Here's a start:

EvalPts:= proc(
    a::algebraic, h::algebraic, n::posint, 
    {method::identical(left, right, midpoint):= right}
)
local x:= a-(if method='left' then h elif method='midpoint' then h/2 else 0 fi);
    [to n do x+= h od]
end proc
:
RiemannSum:= proc(
    f, intv::name=range, n::posint, 
    {
        method::identical(left, right, midpoint):= right, 
        output::identical(value, plot):= value
    }
)
local 
    a:= lhs(rhs(intv)), b:= rhs(rhs(intv)), h:= (b-a)/n,
    E:= EvalPts(a, h, n, ':-method'= method), F:= eval~(f, lhs(intv)=~ E),
    X:= [a, EvalPts(a, h, n)[]]
;
    if output='value' then return h*add(F) fi;
    plot(
        [f, [seq]([[X[k],0],[X[k],F[k]],[X[k+1],F[k]],[X[k+1],0]][], k= 1..n)],
        intv, color= [red, gray],
        caption= typeset( 
            "The ", method, " Riemann sum for ", f, " on ", intv,
            " using ", n, " evaluations is ", evalf(h*add(F)). "."
        )
    )
end proc
:
RiemannSum(1/(1+x^2), x= -2..2, 9);
                           1101459956
                           ----------
                           498612465 

evalf(%);
                          2.209050181

RiemannSum(1/(1+x^2), x= -2..2, 9, method= midpoint, output= plot);

Like this

<xl12, xl13, xl23> =~ 
   LinearAlgebra:-LinearSolve(<1,1,-1; 1,-1,1; -1,1,1>/2, <X1,X2,X3>)
; 

You almost had it. You defined RHS but you didn't use it in the plot. You used sys instead. You also have a typo in the bounds for y1.

You only need 1 loop. I can't even conceive of using two.

Edit: Upon further thought, I do see that two nested loops could be used: The outer loop checks convergence conditions, and the inner loop does the updating of the solution vector. 

You need a multiplication sign after the first x. A space may also work, but I prefer an explicit sign.