num=5 (* Rozmiar siatki *)

n=(num-2)^2 (* Rozmiar macierzy *)

generatorMacierzy[num_]:=Module[{diagN},
diagN=ConstantArray[Join[Table[-1,num-3],{0}],num-2]//Flatten//Most;
SparseArray[{
Band[{1,1}]->4,
Band[{1,2}]->diagN,
Band[{2,1}]->diagN,
Band[{1,num-1}]->-1,
Band[{num-1,1}]->-1},{(num-2)^2,(num-2)^2}]
];

m=generatorMacierzy[num] (* Macierz 4x4 *)

MatrixForm[m]

Normal[m] (* To można przekleić do kodu w C *)

{LU,permutations,condition}=LUDecomposition[m]

Import["lu.m"]

LU//MatrixForm

U=UpperTriangularize[LU];MatrixForm[U]

L=LowerTriangularize[LU,-1]+IdentityMatrix[(num-2)^2];MatrixForm[L]

L.U-m[[permutations]]

LUlapack//MatrixForm

Ulapack=UpperTriangularize[LUlapack]

Llapack=LowerTriangularize[LUlapack,-1]+IdentityMatrix[(num-2)^2]

Llapack.Ulapack-m[[pivot]]

%/$MachineEpsilon//MatrixForm