\documentclass[a4paper]{article}
\usepackage{pgfplots} % fuer plots

\usepackage{amsmath}
\usepackage{amssymb}

\def\showfourxfour{%
	\[
	\def\entry##1,##2{%
		\pgfplotsmatrixselect##1,##2\of\pgfplotsmatrix\to\elem
		\pgfmathprintnumber\elem
	}%
	\begin{bmatrix}
		\entry0,0 & \entry0,1 & \entry0,2 & \entry0,3 \\
		\entry1,0 & \entry1,1 & \entry1,2 & \entry1,3 \\
		\entry2,0 & \entry2,1 & \entry2,2 & \entry2,3 \\
		\entry3,0 & \entry3,1 & \entry3,2 & \entry3,3 
	\end{bmatrix}
	\]
}%
\def\showvectorfour#1{%
	\def\showvectorfourneedsclose{0}%
	\ifmmode\else \[ \def\showvectorfourneedsclose{1}\fi
	\def\entry##1{%
		\pgfplotsarrayselect##1\of#1\to\elem
		\pgfmathprintnumber\elem
	}%
	\begin{bmatrix}
		\entry0 & \entry1 & \entry2 & \entry3 \\
	\end{bmatrix}
	\if1\showvectorfourneedsclose \]\fi
}%

\begin{document}
\parindent=0pt
\parskip=\baselineskip

{
\pgfplotsmatrixnewempty\pgfplotsmatrix
\pgfplotsmatrixresize\pgfplotsmatrix44

\pgfplotsmatrixset 0,0\of\pgfplotsmatrix\to{1}
\pgfplotsmatrixset 0,1\of\pgfplotsmatrix\to{2}
\pgfplotsmatrixset 0,2\of\pgfplotsmatrix\to{3}
\pgfplotsmatrixset 0,3\of\pgfplotsmatrix\to{4}

\pgfplotsmatrixset 1,0\of\pgfplotsmatrix\to{5}
\pgfplotsmatrixset 1,1\of\pgfplotsmatrix\to{6}
\pgfplotsmatrixset 1,2\of\pgfplotsmatrix\to{7}
\pgfplotsmatrixset 1,3\of\pgfplotsmatrix\to{8}

\pgfplotsmatrixset 2,0\of\pgfplotsmatrix\to{-1}
\pgfplotsmatrixset 2,1\of\pgfplotsmatrix\to{4}
\pgfplotsmatrixset 2,2\of\pgfplotsmatrix\to{3}
\pgfplotsmatrixset 2,3\of\pgfplotsmatrix\to{2}

\pgfplotsmatrixset 3,0\of\pgfplotsmatrix\to{5}
\pgfplotsmatrixset 3,1\of\pgfplotsmatrix\to{6}
\pgfplotsmatrixset 3,2\of\pgfplotsmatrix\to{2}
\pgfplotsmatrixset 3,3\of\pgfplotsmatrix\to{-1}


INPUT MATRIX:

\tracingmacros=2 \tracingcommands=2

\showfourxfour

LU result (encoded):
\pgfplotsmatrixLUdecomp\pgfplotsmatrix\perm\pgfplotsperm\sign\pgfplotspermsign

\showfourxfour

with permsign \pgfplotspermsign\ and perm

\showvectorfour\pgfplotsperm

Solve $Ax = b$:
\pgfplotsarraynew\pgfplotsrhs{1,5,-1,5}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixLUbacksubst\pgfplotsmatrix\perm\pgfplotsperm\inout\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{2,6,4,6}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixLUbacksubst\pgfplotsmatrix\perm\pgfplotsperm\inout\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{3,7,3,2}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixLUbacksubst\pgfplotsmatrix\perm\pgfplotsperm\inout\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{4,8,2,-1}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixLUbacksubst\pgfplotsmatrix\perm\pgfplotsperm\inout\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

}


{
Once more again, this time using {\tt \string\pgfplotsmatrixsolveLEQS}:
\pgfplotsmatrixnewempty\pgfplotsmatrix
\pgfplotsmatrixresize\pgfplotsmatrix44

\pgfplotsmatrixset 0,0\of\pgfplotsmatrix\to{1}
\pgfplotsmatrixset 0,1\of\pgfplotsmatrix\to{2}
\pgfplotsmatrixset 0,2\of\pgfplotsmatrix\to{3}
\pgfplotsmatrixset 0,3\of\pgfplotsmatrix\to{4}

\pgfplotsmatrixset 1,0\of\pgfplotsmatrix\to{5}
\pgfplotsmatrixset 1,1\of\pgfplotsmatrix\to{6}
\pgfplotsmatrixset 1,2\of\pgfplotsmatrix\to{7}
\pgfplotsmatrixset 1,3\of\pgfplotsmatrix\to{8}

\pgfplotsmatrixset 2,0\of\pgfplotsmatrix\to{-1}
\pgfplotsmatrixset 2,1\of\pgfplotsmatrix\to{4}
\pgfplotsmatrixset 2,2\of\pgfplotsmatrix\to{3}
\pgfplotsmatrixset 2,3\of\pgfplotsmatrix\to{2}

\pgfplotsmatrixset 3,0\of\pgfplotsmatrix\to{5}
\pgfplotsmatrixset 3,1\of\pgfplotsmatrix\to{6}
\pgfplotsmatrixset 3,2\of\pgfplotsmatrix\to{2}
\pgfplotsmatrixset 3,3\of\pgfplotsmatrix\to{-1}


INPUT MATRIX:

\tracingmacros=2 \tracingcommands=2

\showfourxfour

Solve $Ax = b$:
\pgfplotsarraynew\pgfplotsrhs{1,5,-1,5}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixsolveLEQS\pgfplotsmatrix=\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{2,6,4,6}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixsolveLEQS\pgfplotsmatrix=\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{3,7,3,2}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixsolveLEQS\pgfplotsmatrix=\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]

\pgfplotsarraynew\pgfplotsrhs{4,8,2,-1}
\[ b= \showvectorfour\pgfplotsrhs
\leadsto A^{-1} b =
\pgfplotsmatrixsolveLEQS\pgfplotsmatrix=\pgfplotsrhs
\showvectorfour\pgfplotsrhs \]
}
\end{document}
