##########
# Solve a linear system of equations
# using the SOR (successive over-relaxation) method
##########

# Set the matrix system Ax=b:

N:=2:

A:= array([[5,3],[-1,4]]):

b:= array([8,3]):


# Set initial guess and 
# convergence tolerance

x[1]:= 0.:  x[2]:=0.:

eps:= 0.001;
err:= 99.;

# Set the relaxation parameter:

w:= 0.6;

# Iterate to solution

for k from 1 while err > eps do
	old1:= x[1];
	old2:= x[2];
	for i from 1 to N do
	mysum:=0:
		for j from 1 to N do
		if (i<>j) then mysum:= mysum - A[i,j]*x[j] end if;          # <> means NOT EQUAL TO
		end do;	
# And relax...
	x[i]:= (1-w)*x[i] + w*(mysum+b[i])/A[i,i];

	end do;
print("k=",k);
print("x1=",x[1]);
print("x2=",x[2]);
#
# Compute  a measure of the error
err:= sqrt( (x[1]-old1)^2 + (x[2]-old2)^2 );
end do;