Dear Users!
First, I define the following polynomial as:

restart; with(LinearAlgebra); nu := 1/2; M1 := 3; M2 := 3; #(any value of M1 and M2)
for k1 from 0 while k1 <= M1-1 do for k2 from 0 while k2 <= M2-1 do
SGP[M2*k1+k2+1] := simplify(sum((-1)^(k1-i1)*GAMMA(k1+i1+2*nu)*GAMMA(nu+1/2)*x^i1*(sum((-1)^(k2-i2)*GAMMA(k2+i2+2*nu)*GAMMA(nu+1/2)*t^i2/(GAMMA(i2+nu+1/2)*factorial(k2-i2)*factorial(i2)*GAMMA(2*nu)), i2 = 0 .. k2))/(GAMMA(i1+nu+1/2)*factorial(k1-i1)*factorial(i1)*GAMMA(2*nu)), i1 = 0 .. k1))
end do end do;
Lambda := `<,>`(seq(SGP[i], i = 1 .. M1*M2));

then, I want to define a square matrix A of order M1M2 by M1M2 after collocating SGP[k](x,y) at x=(i-1)/(M1-1) and y=(j-1)/(M2-1) for i=1,2,3,...M1, j=1,2,3,...M2.

For example, for M1=2 and M3=3 this matrix A is given as:

A:=Matrix(6, 6, {(1, 1) = SGP[1](0, 0), (1, 2) = SGP[2](0, 0), (1, 3) = SGP[3](0, 0), (1, 4) = SGP[4](0, 0), (1, 5) = SGP[5](0, 0), (1, 6) = SGP[6](0, 0), (2, 1) = SGP[1](0, 1), (2, 2) = SGP[2](0, 1), (2, 3) = SGP[3](0, 1), (2, 4) = SGP[4](0, 1), (2, 5) = SGP[5](0, 1), (2, 6) = SGP[6](0, 1), (3, 1) = SGP[1](1/2, 0), (3, 2) = SGP[2](1/2, 0), (3, 3) = SGP[3](1/2, 0), (3, 4) = SGP[4](1/2, 0), (3, 5) = SGP[5](1/2, 0), (3, 6) = SGP[6](1/2, 0), (4, 1) = SGP[1](1/2, 1), (4, 2) = SGP[2](1/2, 1), (4, 3) = SGP[3](1/2, 1), (4, 4) = SGP[4](1/2, 1), (4, 5) = SGP[5](1/2, 1), (4, 6) = SGP[6](1/2, 1), (5, 1) = SGP[1](1, 0), (5, 2) = SGP[2](1, 0), (5, 3) = SGP[3](1, 0), (5, 4) = SGP[4](1, 0), (5, 5) = SGP[5](1, 0), (5, 6) = SGP[6](1, 0), (6, 1) = SGP[1](1, 1), (6, 2) = SGP[2](1, 1), (6, 3) = SGP[3](1, 1), (6, 4) = SGP[4](1, 1), (6, 5) = SGP[5](1, 1), (6, 6) = SGP[6](1, 1)});

for M1=3 and M=2 the matrix A is given as:

A:=Matrix(6, 6, {(1, 1) = SGP[1](0, 0), (1, 2) = SGP[2](0, 0), (1, 3) = SGP[3](0, 0), (1, 4) = SGP[4](0, 0), (1, 5) = SGP[5](0, 0), (1, 6) = SGP[6](0, 0), (2, 1) = SGP[1](0, 1/2), (2, 2) = SGP[2](0, 1/2), (2, 3) = SGP[3](0, 1/2), (2, 4) = SGP[4](0, 1/2), (2, 5) = SGP[5](0, 1/2), (2, 6) = SGP[6](0, 1/2), (3, 1) = SGP[1](0, 1), (3, 2) = SGP[2](0, 1), (3, 3) = SGP[3](0, 1), (3, 4) = SGP[4](0, 1), (3, 5) = SGP[5](0, 1), (3, 6) = SGP[6](0, 1), (4, 1) = SGP[1](1, 0), (4, 2) = SGP[2](1, 0), (4, 3) = SGP[3](1, 0), (4, 4) = SGP[4](1, 0), (4, 5) = SGP[5](1, 0), (4, 6) = SGP[6](1, 0), (5, 1) = SGP[1](1, 1/2), (5, 2) = SGP[2](1, 1/2), (5, 3) = SGP[3](1, 1/2), (5, 4) = SGP[4](1, 1/2), (5, 5) = SGP[5](1, 1/2), (5, 6) = SGP[6](1, 1/2), (6, 1) = SGP[1](1, 1), (6, 2) = SGP[2](1, 1), (6, 3) = SGP[3](1, 1), (6, 4) = SGP[4](1, 1), (6, 5) = SGP[5](1, 1), (6, 6) = SGP[6](1, 1)});

and for M1=M2=3 the matrix A should be the following form:

A:=Matrix(9, 9, {(1, 1) = SGP[1](0, 0), (1, 2) = SGP[2](0, 0), (1, 3) = SGP[3](0, 0), (1, 4) = SGP[4](0, 0), (1, 5) = SGP[5](0, 0), (1, 6) = SGP[6](0, 0), (1, 7) = SGP[7](0, 0), (1, 8) = SGP[8](0, 0), (1, 9) = SGP[9](0, 0), (2, 1) = SGP[1](0, 1/2), (2, 2) = SGP[2](0, 1/2), (2, 3) = SGP[3](0, 1/2), (2, 4) = SGP[4](0, 1/2), (2, 5) = SGP[5](0, 1/2), (2, 6) = SGP[6](0, 1/2), (2, 7) = SGP[7](0, 1/2), (2, 8) = SGP[8](0, 1/2), (2, 9) = SGP[9](0, 1/2), (3, 1) = SGP[1](0, 1), (3, 2) = SGP[2](0, 1), (3, 3) = SGP[3](0, 1), (3, 4) = SGP[4](0, 1), (3, 5) = SGP[5](0, 1), (3, 6) = SGP[6](0, 1), (3, 7) = SGP[7](0, 1), (3, 8) = SGP[8](0, 1), (3, 9) = SGP[9](0, 1), (4, 1) = SGP[1](1/2, 0), (4, 2) = SGP[2](1/2, 0), (4, 3) = SGP[3](1/2, 0), (4, 4) = SGP[4](1/2, 0), (4, 5) = SGP[5](1/2, 0), (4, 6) = SGP[6](1/2, 0), (4, 7) = SGP[7](1/2, 0), (4, 8) = SGP[8](1/2, 0), (4, 9) = SGP[9](1/2, 0), (5, 1) = SGP[1](1/2, 1/2), (5, 2) = SGP[2](1/2, 1/2), (5, 3) = SGP[3](1/2, 1/2), (5, 4) = SGP[4](1/2, 1/2), (5, 5) = SGP[5](1/2, 1/2), (5, 6) = SGP[6](1/2, 1/2), (5, 7) = SGP[7](1/2, 1/2), (5, 8) = SGP[8](1/2, 1/2), (5, 9) = SGP[9](1/2, 1/2), (6, 1) = SGP[1](1/2, 1), (6, 2) = SGP[2](1/2, 1), (6, 3) = SGP[3](1/2, 1), (6, 4) = SGP[4](1/2, 1), (6, 5) = SGP[5](1/2, 1), (6, 6) = SGP[6](1/2, 1), (6, 7) = SGP[7](1/2, 1), (6, 8) = SGP[8](1/2, 1), (6, 9) = SGP[9](1/2, 1), (7, 1) = SGP[1](1, 0), (7, 2) = SGP[2](1, 0), (7, 3) = SGP[3](1, 0), (7, 4) = SGP[4](1, 0), (7, 5) = SGP[5](1, 0), (7, 6) = SGP[6](1, 0), (7, 7) = SGP[7](1, 0), (7, 8) = SGP[8](1, 0), (7, 9) = SGP[9](1, 0), (8, 1) = SGP[1](1, 1/2), (8, 2) = SGP[2](1, 1/2), (8, 3) = SGP[3](1, 1/2), (8, 4) = SGP[4](1, 1/2), (8, 5) = SGP[5](1, 1/2), (8, 6) = SGP[6](1, 1/2), (8, 7) = SGP[7](1, 1/2), (8, 8) = SGP[8](1, 1/2), (8, 9) = SGP[9](1, 1/2), (9, 1) = SGP[1](1, 1), (9, 2) = SGP[2](1, 1), (9, 3) = SGP[3](1, 1), (9, 4) = SGP[4](1, 1), (9, 5) = SGP[5](1, 1), (9, 6) = SGP[6](1, 1), (9, 7) = SGP[7](1, 1), (9, 8) = SGP[8](1, 1), (9, 9) = SGP[9](1, 1)});

Please help to define general matrix A for any values of M1 and M2. I will be very grateful for you. 

Dear Users!
Hope everyone fine here. I tried (as given bellow) to find the solution of nonlinear system of PDEs via FDM. To solve system of nonlinear equations I used newton raphson method and for higher value of like Mx > 8 the matrix G and G1 (mentioned as red) take alot of time. Can anyone help me to reduce the computational time? Becuase I have to evealuate the solution for Mx = 50.

restart; Digits := 30; with(LinearAlgebra);
T := 1; L := 3; N := 30; Mx := 5; `&Delta;x` := L/(1.*Mx); `&Delta;t` := T/(1.*N);
for i from 0 while i <= Mx do
u[i, 0] := 0.; u[i, -1] := u[i, 1]; tau[i, 0] := 0.; theta[i, 0] := 0.; theta[i, -1] := theta[i, 1]
end do;
for n from 0 while n <= N do u[0, n] := 0.; u[Mx, n] := 0.; theta[0, n] := 1.; theta[Mx, n] := 0.
end do;
for n from 0 while n <= N-1 do
print("Simulation in proccess at time-level n", n+1);
for i while i <= Mx-1 do
Ru[i, n] := simplify((u[i+1, n+1]-u[i+1, n])/`&Delta;t`+(u[i+1, n+1]-2*u[i+1, n]+u[i+1, n-1])/`&Delta;t`^2-(u[i+1, n+1]-2*u[i, n+1]+u[i-1, n+1])/`&Delta;x`^2+25.*(u[i+1, n+1]+(u[i+1, n+1]-u[i+1, n])/`&Delta;t`)-1.5*(theta[i, n]+(theta[i, n+1]-theta[i, n])/`&Delta;t`));
`R&theta;`[i, n] := simplify((theta[i+1, n+1]-theta[i+1, n])/`&Delta;t`+(theta[i+1, n+1]-2*theta[i+1, n]+theta[i+1, n-1])/`&Delta;t`^2-(theta[i+1, n+1]-2*theta[i, n+1]+theta[i-1, n+1])/((15.)*`&Delta;x`^2)-((u[i, n+1]-u[i-1, n+1])/`&Delta;x`)^2/(3.)) end do;
for i while i <= Mx-1 do
`R&tau;`[i, n] := simplify(tau[i+1, n+1]+(tau[i+1, n+1]-tau[i+1, n])/`&Delta;t`-1.5^(-1/4)*(u[i+1, n+1]-u[i, n+1])/`&Delta;x`)
end do;
Sys := `<,>`(seq(Ru[i, n], i = 1 .. Mx-1), seq(`R&tau;`[i, n], i = 1 .. Mx-1), seq(`R&theta;`[i, n], i = 1 .. Mx-1));
V := `<,>`(seq(u[i, n+1], i = 1 .. Mx-1), seq(theta[i, n+1], i = 1 .. Mx-1), seq(tau[i, n+1], i = 2 .. Mx));
G := Matrix(3*(Mx-1), proc (i, j) options operator, arrow; diff(Sys[i], V[j]) end proc); G1 := MatrixInverse(G);
X[n, 0] := Vector(1 .. 3*(Mx-1), 1);
for k1 from 0 to r do
X[n, k1+1] := eval(V-G1 . Sys, Equate(V, X[n, k1]))
end do;
Sol[n] := Equate(V, X[n, r+1]); assign(op(Sol[n]));
if n > 0 then
U := eval(`<,>`(seq(u[i1, n+1], i1 = 1 .. Mx)-seq(u[i1, n], i1 = 1 .. Mx))); Noru[n+1] := Norm(%, 2); print("L[&infin;] norm of &tau;(x,y,t) at time level = ", %);
Theta := eval(`<,>`(seq(theta[i1, n+1], i1 = 0 .. Mx)-seq(theta[i1, n], i1 = 0 .. Mx))); `Nor&theta;`[n+1] := Norm(%, 2); print("L[&infin;] norm of &theta;(x,y,t) at time level = ", %) else print("n < 0")
end if end do

Special request to:
@acer @Carl Love @Kitonum @Preben Alsholm

Dear Users!

I want to find the solution of the solution of PDEs as given bellow:

restart;
PDE1 := diff(u(y, t), t)+diff(u(y, t), t, t) = diff(u(y, t), y, y)-u(y, t);
PDE2 := v(y, t)+diff(tau(y, t), t) = diff(u(y, t), y);
ICandBC := {tau(y, 0) = 0, u(0, t) = 0, u(3, t) = 0, u(y, 0) = 0, (D[2](u))(y, 0) = 0};
pds := pdsolve({PDE1, PDE2}, ICandBC, numeric);

But got the following error

Error, (in pdsolve/numeric/process_PDEs) number of dependent variables and number of PDE must be the same
Kindly help me to fix this error. I shall be very thankful. 

Special request to:
@acer @Carl Love @Kitonum @Preben Alsholm

Dear Users! 

I want to evaluate following integral (f[i]'s) for different values of rho as shown. But find some problem and can't find f[3], f[4], f[5] and f[6]. If the integral from direct way is hard to calculate please guide me how I can evaluate it using numerical integration?

restart; HAM := [1, 2, 3, 4, 5, 6]; U := y^6*sin(x); alpha := 1/3;

for i while i <= nops(HAM) do

rho := op(i, HAM);

f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), q = 0 .. y), p = 0 .. x)

end do

Thanks in advance!

Special request to @acer @Carl Love @Kitonum @Preben Alsholm

Dear Users!

Hope you would be fine. I have some problem in execution the last loops (highlighted as red) where sumation is present. When NN>3 it takes alot of time more than 12 hours. Is there any alternative command to reduce the query. I am waiting for your response. Thanks in advance. 

restart; with(LinearAlgebra); Digits := 30; NN := 2; nu := 1; M1 := NN; M2 := NN; M3 := NN;

for k1 from 0 while k1 <= M1-1 do for k2 from 0 while k2 <= M2-1 do for k3 from 0 while k3 <= M3-1 do

SGP[M3*(M2*k1+k2)+k3+1] := simplify(sum((-1)^(k1-i1)*GAMMA(k1+i1+2*nu)*x^i1*(sum((-1)^(k2-i2)*GAMMA(k2+i2+2*nu)*y^i2*(sum((-1)^(k3-i3)*GAMMA(k3+i3+2*nu)*z^i3/(GAMMA(i3+nu+1/2)*factorial(k3-i3)*factorial(i3)), i3 = 0 .. k3))/(GAMMA(i2+nu+1/2)*factorial(k2-i2)*factorial(i2)), i2 = 0 .. k2))/(GAMMA(i1+nu+1/2)*factorial(k1-i1)*factorial(i1)), i1 = 0 .. k1)) end do end do end do;

SGPxyz := `<,>`(seq(seq(seq(SGP[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

Lambda := `<,>`(seq(seq(seq(chi[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

for i while i <= NN^3 do for j while j <= NN^3 do for k while k <= NN^3 do

q[i, j, k] := int(int(int(SGP[i]*SGP[j]*SGP[k]*(-x^2+x)^(nu-1/2)*(-y^2+y)^(nu-1/2)*(-z^2+z)^(nu-1/2), z = 0 .. 1), y = 0 .. 1), x = 0 .. 1) end do end do end do;

U := Matrix(NN^3, NN^3, 0);

for j while j <= NN^3 do for k while k <= NN^3 do U[j, k] := simplify(sum(chi[i1]*q[i1, j, k], i1 = 1 .. NN^3)) end do end do;

F := simplify(evalm(U));

Special request to @acer @Carl Love @Kitonum @Preben Alsholm