f[x_,y_]:=-6*x-2;
p[x_]:=x^3;
q[x_]:=16+8*x+x^3;
r[y_]:=y^2;
s[y_]:=64+8*y+y^2;
pres[x_,y_]:=x^3+2*x*y+y^2;
a=4;

n=3;
h=a/(n+1);
A=IdentityMatrix[n*n];
PS=Table[0,{i,1,n*n}];
Print[MatrixForm[A]];
Print[PS];

Do[
{Do[
{k=i*n+j,
PS[[k]]=N[h^2*f[(i+1)*h,h*j]],
A[[k,k]]=4,
If[(i==0),
{PS[[k]]=PS[[k]]+r[j*h],Print[j],Print[MatrixForm[A]]},
A[[k,k-n]]=-1],
If[(i==(n-1)),
PS[[k]]=PS[[k]]+s[j*h],
A[[k,k+n]]=-1],
If[(j==1),
{PS[[k]]=PS[[k]]+p[(i+1)*h],Print[MatrixForm[A]]},
{Print[j],Print[MatrixForm[A]],A[[k,k-1]]=-1,Print[j],Print[MatrixForm[A]]}],
If[(j==n),
{PS[[k]]=PS[[k]]+q[(i+1)*h],Print[j],Print[MatrixForm[A]]},
{Print[j],Print[MatrixForm[A]],A[[k,k+1]]=-1}]
},
{j,1,n}],
Print[MatrixForm[A]]},
{i,0,n-1}];

Print[MatrixForm[A]];
Print[PS];

Ries=LinearSolve[A,PS]