A: matrix(
    [2,4,-4,0],
    [1,5,-5,-3],
    [2,3,1,3],
    [1,4,-2,2]
);
B: matrix(
    [12],
    [18],
    [8],
    [8]
);
X: matrix(
    [0],
    [0],
    [0],
    [0]
);
X0: matrix(
    [x1],
    [x2],
    [x3],
    [x4]
);
n: length(A);
decimales:4;
prec:10^decimales;

for k:1 thru n-1 do (
    for i:k+1 thru n do (
        factor:A[i,k]/A[k,k],
         B[i]:B[i]-factor*B[k],
        for j:k thru n do (
            A[i,j]:A[i,j]-factor*A[k,j]
        )
       
    )
);
X[n]:B[n]/A[n,n];
for i:n-1 thru 1 step -1 do (
    sum:B[i],
    for j:i+1 thru n do (
        sum:sum-A[i,j]*X[j],
    X[i]:sum/A[i,i]   
    )
);
print("Sistema = ",float(round(A*prec)/prec),X0,"=",float(B))$
print("Solución = ",float(X))$