Die linearen Gleichungssysteme aus Gl. 4.23 resultieren in spärlich besetzten Matrizen. Zur Lösung dieser Systeme kann wahlweise ein iterativer Gleichungslöser oder ein direktes Verfahren verwendet werden. Beide Methoden werden von der Entwicklungsumgebung als Bibliotheksfunktion zur Verfügung gestellt [Hei93].
Der direkte Gleichungslöser arbeitet nach dem Verfahren der Gauß-Elimination [Lee68] und liefert im Gegensatz zu iterativen Methoden auch bei schlecht konditionieren Gleichungssystemen zuverlässige Lösungen. Der benötigte Speicher und Rechenaufwand bei dieser Methode ist jedoch speziell bei dreidimensionalen Problemen substantiell größer als bei iterativen Methoden, sodaß das direkte Verfahren hauptsächlich zu Testzwecken eingesetzt wird.
Die im Zuge der Lösung der Diffusionsgleichung auftretenden Gleichungssysteme lassen sich aufgrund ihrer Hauptdiagonaldominanz sehr effizient mit iterativen Methoden, wie der Methode der sukzessiven Überrelaxation (successive over-relaxation SOR), der konjugierten Gradienten (conjugate gradient CG) oder der bi-konjugierten Gradienten (bi-conjugate gradient Bi-CG) lösen.
Der verwendete iterative Gleichungslöser arbeitet nach dem stabilisierten Verfahren der bi-konjugierten Gradienten (BiCGStab [vdV93]) mit einer inkompletten L-U Faktorisierung [Saa90] als Präkonditionierverfahren und einer automatischen Skalierung der Gleichungen [Fis94c]. Im Zuge der inkompletten L-U-Faktorisierung werden zusätzliche Nebendiagonaleinträge in der Systemmatrix generiert, deren Anzahl durch eine sog. Dropping-Strategie begrenzt wird. Alle Einträge, die aufgrund ihres geringen Zahlenwertes keine wesentliche Information beinhalten, werden dabei vernachlässigt. Eine genaue Beschreibung der Implementierung ist in [Fis94a] und [Hei93] zu finden.
Die praktische Anwendung dieses Verfahrens auf Diffusionsprobleme hat gezeigt, daß der Aufwand zur Lösung einer nichtlinearen Iteration in etwa proportional mit der Anzahl der freien Variablen steigt.