Use mathematical' 
        s function "Solve";

    A = {{1, 3, 1, -2}, {0, 2, -1, -3}, {0, 1, 2, -2}, {0, 
            0, 1, 3}};

    Print["A=", 
        A // MatrixForm];
  

    H = {{2}, {1}, {2}, {0}};

    X = {{x1}, {x2}, {x3}, {x4}};

    Xp = Solve[
            A . X == H, {x1, x2, x3, 
              x4}][[1]];

    x1 = x1 /. Xp; x2 = x2 /. Xp; x3 = x3 /. Xp; 
    x4 = x4 /. Xp;

    Xp = {{x1}, {x2}, {x3}, {x4}};

    Print["Xp=", Xp // MatrixForm];

    Print["A.Xp=", A . Xp // MatrixForm];

Use Cramer' 
        s rule to find the particular solution;
A = {{1, 3, 1, -2}, {0, 2, -1, -3}, {0, 1, 
            2, -2}, {0, 0, 1, 3}};
Print["A=", A // MatrixForm];
A1 = {{2, 3, 1, -2}, {1, 2, -1, -3}, {2, 1, 
            2, -2}, {0, 0, 1, 3}};
A2 = {{1, 2, 1, -2}, {0, 1, -1, -3}, {0, 2, 
            2, -2}, {0, 0, 1, 3}};
A3 = {{1, 3, 2, -2}, {0, 2, 1, -3}, {0, 1, 2, -2}, {0, 0, 
            0, 3}};
A4 = {{1, 3, 1, 2}, {0, 2, -1, 1}, {0, 1, 2, 2}, {0, 0, 1, 
            0}};

Print["A=", A // MatrixForm];
Print["A1=", A1 // MatrixForm];
Print["A2=", A2 // MatrixForm];
Print["A3=", A3 // MatrixForm];
Print["A4=", A4 // MatrixForm];
x1 = Det[A1]/Det[A]; x2 = Det[A2]/Det[A]; x3 = Det[A3]/Det[A]; 
    x4 = Det[A4]/Det[A];
Xp = {{x1}, {x2}, {x3}, {x4}};
Print["Xp=", 
        Xp // MatrixForm];

Now,  solve the Homogeneous part, A . X = 0;

Print["Det[A]=", Det[A]];

X = {{xx1}, {xx2}, {xx3}, {xx4}};
sol = Solve[A . X  ==  0];
Print["sol=", sol] // MatrixForm;

The only sol for the homogeneous part is the trivial solution, Xh = 0; Hence the complete solution is X=Xp=(({{-(7/16)}, {(1/2)}, {(9/16)}, {-(3/16)}})).