This is a very simple Pascal program I wrote in mid 1990s, which solves a sparse NxN linear system using the Succesive Overrelaxation Method (SOR), which in turn is based on Gauss Seidel Method.
Both methods are iterative.

N is initially set at maximum …64.
You need the W factor for the SOR to complete successful. The theory assures us that it must be W=1.062.

PROGRAM SOR;
uses crt;
VAR
 A: array [1..64,1..64] of real;
 B,X,XT: array [1..64] of real;
 I,J,K,L,M,N,ROW1,COL1,REP,MAX : integer;
 SUM,W: real;
 LABEL ALPHA,BETA;
 
begin
 
 clrscr;
 
 SUM:=0;
 
 WRITELN ('PLEASE ENTER N:');
 READLN ( N);
 
 CLRSCR;
 
 WRITELN ('PLEASE, ENTER MATRIX A NxN, LINE BY LINE:');
 FOR I:=1 TO N DO
  FOR J:=1 TO N DO
  BEGIN
   READ (A[I,J]);
  END;   
 
 CLRSCR;
 
 WRITELN ('PLEASE, ENTER MATRIX B Nx1 :');
 FOR I:=1 TO N DO
  BEGIN 
   READ (B[I]);
  END; 
 
 CLRSCR;
 
 WRITELN ('PLEASE, ENTER YOUR VECTOR X:');
 FOR I:=1 TO N DO
  BEGIN
   READ (X[I]);
 END;
 
 CLRSCR;
 
 WRITELN ('PLEASE ENTER NUMBER OF REPEATS:');
 READ (MAX);
 
 CLRSCR;
 
 WRITELN ('PLEASE ENTER W:');
 READ ( W);
 
{ MAIN GAUSS SEIDEL ALGORITHM }
 
 FOR REP:=1 TO MAX DO
  BEGIN
   FOR I:=1 TO N DO
    BEGIN
     SUM:=0;
     IF I=1 THEN GOTO ALPHA;
     FOR J:=1 TO I-1 DO
      BEGIN
       SUM:=SUM+A[I,J]*X[J];
      END;
     ALPHA: IF I=N THEN GOTO BETA;
     FOR K:=I+1 TO N DO
      BEGIN
       SUM:=SUM+A[I,K]*X[K];
      END;
     BETA:X[I]:=((W/A[I,I])*(B[I]-SUM))-(W-1)*X[I];
    END;
  END;
 
 CLRSCR;
 
 WRITELN ('X=');
 FOR I:=1 TO N DO
  BEGIN
  WRITELN (X[I]);
  END;
 
 WRITELN ('any key to continue!');
 READLN;READLN;
 
END.