@peace951 

Also Tried this way: 

f := cos(y);
                         

a := 0;
b := Pi;
c := 0;
d := x;
n := 4;
m := 8;
dx := (b - a)/n;
dy := (d - c)/m;
u := i;
v := j -> c + (j - 1/2)*dy;
APV := sum(sum(f(u(i), v(i)), j = 1 .. m), i = 1 .. n)*dx*dy;
evalf(APV);
                             a := 0

                            b := Pi

                             c := 0

                             d := x

                             n := 4

                             m := 8

                                 1   
                           dx := - Pi
                                 4   

                                 1  
                           dy := - x
                                 8  

                             u := i

  v := proc (j) options operator, arrow; c+(j-1/2)*dy end proc

               1  /        /   1   \           /   3   \
        APV := -- |8 cos(y)|1, -- x| + 8 cos(y)|2, -- x|
               32 \        \   16  /           \   16  /

                     /   5   \           /   7   \\     
           + 8 cos(y)|3, -- x| + 8 cos(y)|4, -- x|| Pi x
                     \   16  /           \   16  //     

                   /         /   1   \            /   3   \
     0.09817477044 |8. cos(y)|1, -- x| + 8. cos(y)|2, -- x|
                   \         \   16  /            \   16  /

                   /   5   \            /   7   \\  
        + 8. cos(y)|3, -- x| + 8. cos(y)|4, -- x|| x
                   \   16  /            \   16  //  

EV := Int*(Int*(f(x, y), x = a .. b), y = c .. d);
    EV := Int (Int (cos(y)(x, y), x = 0 .. Pi), y = 0 .. x)

evalf(EV);
  Int (Int (cos(y)(x, y), x = 0. .. 3.141592654), y = 0. .. x)

@Rouben Rostamian  

I didn't have any intetion to use someelses work as mine: I tried for couple of hours to figure it out but my code hade problem to run and then I wrote it here. here is what I have done: 

f := (x, y) -> cosy;
APV := proc(f, a, b, c, d, n, m)

global A_end, A_even, A_odd, i, ya, yb, B_even, B_odd, j;

h = (b - a)/n;

A_end := 0;

A_even := 0;

A_odd := 0;

for i from 0 to n do

x = a + i*h;

ya := c(x);

yb := d(x);

k = (d(x) - c(x))/m;

B_end = f(x, ya) + f(x, yb);

B_even := 0; B_odd := 0;

for j to m - 1 do

y = j*k + ya; z = f(x, y);

if rem(j, 2) = 0 then B_even = B_even + z;

else B_odd = B_odd + z;

end if;

end do;

L = 1/3*(B_end + 4*B_odd + 2*B_even)*k;

if i = 0 or i = n then A_end = A_end + L;

else

if rem(i, 2) = 0 then A_even = A_even + L;

else A_odd = A_odd + L;

end if;

end if;

end do;

Dint = 1/3*(A_end + 4*A_odd + 2*A_even)*h;

print(Dint);

end proc;

evalf(APV(f, 0, pi, 0, x, 4, 8));
Error, (in APV) invalid input: rem uses a 3rd argument, x, which is missing

Exact*Value;
EV := int(int(f, y = c .. d), x = a .. b);
                              1        2
                        EV := - cosy Pi 
                              2         

evalf(EV);
                        4.934802202 cosy

@Kitonum 

Can't find any answer!!!

@mmcdara 

Thank you so much.

@Kitonum 

 there is no 6 in the matrix in the question, [1, 2, 3] in the first row and [3, 4, 5] in the 2nd row. should I use randperm($1..5)? 

3 repeats? 

Thank you so much.