(*zadanie, oblast, prava strana, okrajove podm.*)
f[x_,y_]:=-6x-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;
(*vypocet*)

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],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,h*j]]},{j,1,n}],{i,0,n-1}];
Print[pre];

vysledky=Table[{Ries[[k]],pre[[k]],Abs[(Ries[[k]]-pre[[k]])]},{k,1,n*n}];
Print["   num. ries.    pres. ries.    odchylka"]
PaddedForm[TableForm[vysledky],{5,4}]

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

<<"BarCharts`";
<<"Histograms`";

aprox=Table[i,{i,1,n*n},{j,1,n}]

graprox=BarChart3D[aprox]