ULOHA :
-Δu= - 6x -2
u(0,y) = y^2
u(x,0) = x^3
u(4,y) = 64 + 8y + y^2
u(x,4) = 16 + 8x + x^3
u(x,y) = x^3 + 2xy + y^2

f[x_, y_] := -6*x - 2;
p[x_] := x^3;
r[y_] := y^2;
q[x_] := 16 + 8*x + x^3;
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}];

Do[
Do[{
k = 
            i*n + j; 
PS[[k]] = 
            N[h^2*f[(i + 1)*h, j*h]]; 
A[[k, k]] = 
            4; 
If[(i == 0), 
            PS[[k]] = PS[[k]] + r[j*h], 
            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], 
            A[[k, k - 1]] = -1]; 
If[(j == 
                n), PS[[k]] = PS[[k]] + q[(i + 1)*h], 
            A[[k, k + 
                    1]] = -1];
}, 
{j, 1, n}], {i, 0, 
        n - 1}]

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

Ries = LinearSolve[A, PS]

pre = Table[0, {i, 1, n*n}];
Do[
        Do[{k = i*n + j, pre[[k]] = N[pres[(i + 1)*h, j*h]]}, {j, 1, 
            n}, {i, 0, n - 1}]];

Print[pre];

vysledky = 
        Table[{i, pre[[i]], Ries[[i]], 
            Abs[(pre[[i]] - Ries[[i]])]}, 
{i, 
            1, n*n}];
Print["cislo uzlu    Pres ries    Num ries    chyba"]
PaddedForm[TableForm[vysledky], {5, 3}]

<< Graphics`Graphics3D`

plopr = Plot3D[pres[x, y], {x, 0, 4}, {y, 0, 4}]

aprox = {Table[Ries[[i]], {i, 1, n}], 
        Table[Ries[[n + i]], {i, 1, n}], 
        Table[Ries[[2*n + 1]], {i, 1, n}]}

graprox = BarChart3D[aprox]

aprpl = ListPlot3D[aprox]