Die Vorwärtselimination ist beendet, das Rücksubstituieren kann gestartet werden. Die dazugehörige Prozedur soll "Rückwärtseinsetzen" heißen. Man benötigt zwei Laufschleifen, eine äußere i- und eine innere j-Schleife. Die i-Schleife durchläuft rückwärts die Werte von n bis eins.
Um den folgenden Schritt besser verstehen zu können, ruft man sich am besten noch einmal das gestaffelte Gleichungssystem aus I.2.1 ins Gedächtnis:

a₁₁x₁ + a₁₂x₂ +...+ a_1nx_n = a_1,n+1
a₂₂x₂ +...+ a_2nx_n = a_2,n+1
..................................
a_nnx_n = a_n,n+1

Löst man nun jeweils die i-te Gleichung nach x_i auf, so erhält man folgendes:

x₁ = (a_1,n+1 - (a₁₂x₂ +...+ a_1nx_n)):a₁₁
x₂ = (a_2,n+1 - (a₂₃x₃ +...+ a_2nx_n)):a₂₂
........................................
x_n = (a_n,n+1 - 0):a_nn

Die j-Schleife berechnet nun für die Zeile i die Summe in der innersten Klammer (oben kursiv). Die j-Schleife durchläuft die Werte i+1 bis n.
Man könnte denken, es gäbe ein Problem, wenn i=n ist, weil dann j=n+1 sein muß und somit auch x_n+1 herangezogen würde, das natürlich nicht existiert. Doch löst sich dieses Problem in Wohlgefallen auf, da die j-Schleife in diesem Fall überhaupt nicht mehr durchlaufen wird, da dann der Startwert j=n+1 größer als der Endwert n ist.
Im weiteren Verlauf der i-Schleife wird nun noch geprüft, ob a_nn=0 und gegebenenfalls wieder ein "Division by Zero" Fehler abgefangen.
Letztendlich wird die Lösung x_i berechnet.
Als Programmtext sieht dies wie folgt aus:

PROCEDURE RückwärtsEinsetzen;                                         
VAR i, j, NullPos: INTEGER;                                           
    Term: REAL;                                                       
BEGIN                                                                 
  FOR i:= n DOWNTO 1 DO                                               
    BEGIN                                                             
      Term:= 0;                                                       
      FOR j:= i + 1 TO n DO                                           
        Term:= Term + Koeffizient[i, j] * Lösung[j];                  
                                                                      
      IF Koeffizient[i, i] = 0 THEN Abbruch;                          
      Lösung[i] := (Koeffizient[i, n + 1] - Term) / Koeffizient[i, i];
    END;                                                              
END;                                                                  

Jetzt steht die simultane Lösung in dem Array "Lösung". Dieses kann man sich am Ende noch mit einer Prozedur "DruckeLösung" ausgeben lassen. Das komplette Programm liegt im Anhang als Quelltext vor.

1.2 Programmbeschreibung

Die Programmbeschreibung wird in die verschiedenen Unterpunkte Variablen und Konstanten, Prozeduren und Funktionen und Hauptprogramm unterteilt.
Das unten beschriebene Programm gibt zusätzlich auf Wunsch noch Rechenzwischenschritte aus, in der Programmbeschreibung wird darauf allerdings nicht eingegangen. Anzumerken ist noch, daß diese Zwischenschrittausgabe sehr einfach gehalten wurde, d.h. es kann zu Ausgaben kommen, die mathematisch nicht korrekt sind (z.b. zwei aufeinanderfolgende Verknüpfungszeichen).
Die einzelnen Schritte werden in der Reihenfolge beschrieben, in der sie auch im Programmtext zu finden sind.