(*-(troj^-1).(xy(troj^-1)u)=9xy^2-12xy*)(*samoadujn *)

u(0,y)=-y^3
u(x,0)=3x^2
u(1,y)=3-y^3
u(x,1)=3x^2-1
na omega=1x1
presne=3x^2-y^3

f[x_,y_]:=9*x*y^2-12*x*y;
p[x_]:=3*x^2;
q[x_]:=3*x^2-1;
r[y_]:=-y^3;
s[y_]:=3-y^3;
t[x_,y_]:=y*x;
pres[x_,y_]:=3*x^2-y^3;
a=1;(*rozmer oblasti*)

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

g1=Plot3D[pres[x,y],{x,0,1},{y,0,1}]

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

aprox=Table[Ries[[n(j-1)+i]],{i,1,n},{j,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]

g2=ListPlot3D[aha,Mesh->All]

Show[g1,g2]