(*TROJUHOLNIK*)
f[x_,y_]:=-8;
p[x_]:=x^2;
r[y_]:=3*y^2;
q[x_]:=x^2+3*(1-x)^2;
pres[x_,y_]:=If[0≤y≤(1-x),x^2+3*y^2,0];
a=1;
(*vypocet*)

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

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

Ries=LinearSolve[A,PS]

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

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

Plot3D[pres[x,y],{x,0,1},{y,0,1},PlotRange->All]

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

aprox=Table[Ries[[i]],{i,1,n}].....

graprox=BarChart3D[aprox]

Do[apr[j]=Table[{h*j,h*i,Ries[[n*(j-1)+i]]},{i,1,n}],{j,1,n}]

Do[Print[apr[j]],{j,1,n}]

a=Table[apr[i],{i,n}]

aha=Flatten[a,1]

ListPlot3D[aha,Mesh->All]

apr[1]=Table[{0,h*j,r[h*j]},{j,0,n+2}]
Do[
apr[i+2]=
Join[{{h*(i+1),0,p[h*(i+1)]}},
Table[{h*(i+1),h*j,
Ries[[i*n-i*(i-1)/2+j]]},
{j,1,n-i}],
{{h*(i+1),h*(n-i+1),q[h*(i+1)]}}],
{i,0,n-1}
]
apr[n+2]={{h*(n+1),0,p[h*(n+1)]},
{h*(n+1),h,q[h*(n+1)]}}
apr[n+3]={{h*(n+2),0,p[h*(n+2)]}}

Do[Print[apr[i]],{i,1,n+3}]

a=Table[apr[i],{i,1,n+3}]

aha=Flatten[a,1]

ListPlot3D[aha,Mesh->All]