@Markiyan Hirnyk 

I think that the code for the projection and the regions are enough for this problem.
Any other solid people can start from here.

@Markiyan Hirnyk 

Yes, the MMA expression seems to be the same (the 13 regions), but that triangle is a nonsense. I am not sure about the boundary points.
The MMA projection  is:

@Markiyan Hirnyk 

A projection of a triangle is a triangle. As can be seen from the first plot (with code) this is false. So, MMA is totally wrong.

@Markiyan Hirnyk 

The MMA plot does not seem to be correct.
Can you plot the (x1,x2) projection?

The 13 regions obtained with Maple are:

Region1
{442/7 < x2, x1 < 80, x2 < 64, 381/2-(7/4)*x2 < x1}
Region2
{432/7 < x2, x1 < 78, x2 < 62, -(7/4)*x2+186 < x1}
Region3
{414/7 < x2, x1 < 78, x2 < 60, -(7/4)*x2+363/2 < x1}
Region4
{438/7 < x2, x1 < 75, x2 < 442/7, -(7/4)*x2+369/2 < x1}, {442/7 <= x2, x2 <= 444/7, 74 < x1, x1 < 75}, {x1 = 74, 442/7 < x2, x2 < 64}, {74 < x1, 444/7 < x2, x1 < -(7/4)*x2+186, x2 < 64}, {x1 = -(7/4)*x2+186, 444/7 < x2, x2 < 64}
Region5
{426/7 < x2, x1 < 75, x2 < 430/7, -(7/4)*x2+363/2 < x1}, {430/7 <= x2, 74 < x1, x1 < 75, x2 < 62}, {x1 = 74, 430/7 < x2, x2 < 62}
Region6
{410/7 < x2, x1 < 76, x2 < 414/7, -(7/4)*x2+357/2 < x1}, {414/7 <= x2, 75 < x1, x1 < 76, x2 < 60}, {x1 = 75, 414/7 < x2, x2 < 60}
Region7
{404/7 < x2, x1 < 76, x2 < 58, -(7/4)*x2+177 < x1}
Region8
{408/7 < x2, x1 < 75, x2 < 412/7, -(7/4)*x2+177 < x1}, {412/7 <= x2, x2 <= 414/7, 74 < x1, x1 < 75}, {x1 = 74, x2 <= 418/7, 412/7 < x2}, {74 < x1, 414/7 < x2, x1 < -(7/4)*x2+357/2, x2 < 418/7}, {x1 = -(7/4)*x2+357/2, 414/7 < x2, x2 < 418/7}
Region9
{x2 = 60, 147/2 < x1, x1 < 74}, {60 < x2, x1 < 74, x2 < 426/7, -(7/4)*x2+357/2 < x1}, {426/7 <= x2, x2 <= 430/7, 72 < x1, x1 < 74}, {x1 = 72, 426/7 < x2, x2 < 62}, {72 < x1, 430/7 < x2, x1 < -(7/4)*x2+363/2, x2 < 62}, {x1 = -(7/4)*x2+363/2, 430/7 < x2, x2 < 62}
Region10
{58 <= x2, 72 <= x1, x2 <= 412/7, x1 < 74}, {x1 = 72, 412/7 < x2, x2 < 60}, {72 < x1, 412/7 < x2, x1 < -(7/4)*x2+177, x2 < 60}, {x1 = -(7/4)*x2+177, 412/7 < x2, x2 < 60}
Region11
{56 <= x2, x2 <= 402/7, 70 < x1, x1 < 72}, {x1 = 70, 56 < x2, x2 < 58}, {70 < x1, 402/7 < x2, x1 < 345/2-(7/4)*x2, x2 < 58}, {x1 = 345/2-(7/4)*x2, 402/7 < x2, x2 < 58}
Region12
{x2 = 54, 70 <= x1, x1 < 72}, {x1 = 70, 54 < x2, x2 < 55}, {54 < x2, 70 < x1, x1 < 333/2-(7/4)*x2, x2 < 55}, {x1 = 333/2-(7/4)*x2, 54 < x2, x2 < 55}
Region13
{x2 = 56, 68 <= x1, x1 <= 137/2}, {x1 = 68, x2 <= 394/7, 56 < x2}, {56 < x2, 68 < x1, x1 < 333/2-(7/4)*x2, x2 < 394/7}, {x1 = 333/2-(7/4)*x2, 56 < x2, x2 < 394/7}

The correct 3d plot:

@Preben Alsholm 

To obtain the 1D version I have pasted the Int in a 2D execution group and then lprint(%).
BTW, the exact and approx value of the integral are computed below.
 

@Carl Love 

OP's function written in a natural and simplified form is

f:=piecewise(0<-w*x*y+u*x and -w*x*y+u*x<1,1)*piecewise(0<-w*x*y+v*y and-w*x*y+v*y<1,1) ;

Now,

int(f,  [x=0..1,y=0..1,u=0..1,v=0..1,w=0..1]);
       7/12

Converting to Heaviside
g:=convert(f,Heaviside);
        g := Heaviside(-w*x*y+u*x)*Heaviside(-w*x*y+v*y)*Heaviside(w*x*y-u*x+1)*Heaviside(w*x*y-v*y+1)
int(g,  [x=0..1,y=0..1,u=0..1,v=0..1,w=0..1]);
    (interrupted)
but ...

int(g,  [u=0..1,v=0..1,w=0..1,x=0..1,y=0..1]);   # ok!
                               7/12

I wanted to check manually the result.
I was surprised to find that Maple is not able to solve the system

[0 < -w*x*y+u*x, -w*x*y+u*x < 1, 0 < -w*x*y+v*y, -w*x*y+v*y < 1, 0 < x, x < 1, 0 < y, y < 1, 0 < u, u < 1, 0 < v, v < 1, 0 < w, w < 1]

using solve, or SolveTools:-SemiAlgebraic

even if it's easy to see that the solution is:

0<x<1, 0<y<1, 0<w<1, w*y<u<1, w*x<v<1

so the integral is indeed
int((1-w*x)*(1-w*y), x=0..1, y=0..1, w=0..1) = 7/12.

@mmcdara 

It's about the user, not the laptop ... :-)

@mmcdara 

Have you tried the code?

discont was added just for a correct plot.

It seems that you actually do not want an answer. In this case you can simply ignore mine.

kernelopts(version);
   Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769

bin.mw

@Klausklabauter 

So, LDLt is actually implemented in

Student:-NumericalAnalysis:-MatrixDecomposition(A, 'method'='LDLt')

Note however that it is too slow (compared with Cholesly for floats) for large matrices.

@Ahmed111 

For p:=1; q:=1; lambda:=2; i:=1  the system has only the trivial solution.
Maple says u1=u2=0  for i = I  too. If you want to check this, try to use the method from the provided worksheet.

@ThU 

Yes, this is what MakeInert does (automatically).

@mmcdara 

I just see that the two parameters are switched compared with some (all?)  standard references.
Have you noticed a bug related to this? As long as Maple is consistent, it seems to be OK.

Even in mathematics the notations are not standardized.
E.g. for some authors N = {0,1,2,...}  but for others N = {1,2,...}.

A much more annoying situation in Maple is that the conventions are not the same everywhere:

spherical coordinates in plot3d:
x = u cos(v) sin(w)
y = u sin(v) sin(w)
z = u cos(w)

spherical coordinates in  VectorCalculus:
x = u cos(w) sin(v)
y = u sin(v) sin(w)
z = u cos(v)

i.e. phi and theta are switched!

@Carl Love 

I agree for frem but not for trunc  which satisfies the identity  frac(x)+trunc(x) = x.
Being built in is not necessarily relevant; probably frac(x) should be more often replaced (internally)  by x - trunc(x).

@Markiyan Hirnyk 
Don't you also think that the last sentence could have been omitted?