(* samoadjungovane OU  - NEUMANNOVE PODMIENKY*)
(* riesi/me -(1/x^2y')'+s/x^4y=6/x^3 , y(1)=0, y(2)=-1.5, presne : y(x)=1/x-x *)
f[x_]:=x^2+2*x*Sin[x]-Cos[x]-1;
n=5;
a=0;
b=Pi;
ya=0;
yb=0;
p[x_]:=x;
q[x_]:=x;
presne[x_]:=Sin[x]+x;

(* vypocet *)
OkrsamoadjNP[n_,xleft_,xright_,Leftbv_,Rightbv_,p_,q_,f_,pres_]:=
(
h=N[(xright-xleft)/(n+1)];
ps=Join[
{f[xleft+h]*h^2+p[xleft+h/2]*Leftbv},
Table[h^2*f[xleft+i*h],{i,2,n}],
{f[xleft+(n+1)*h]*h^2+p[xright-h/2]*Rightbv*2*h}
];
diag=
Table[p[xleft+i*h-h/2]+p[xleft+i*h+h/2]+h^2*q[xleft+i*h],{i,1,n+1}];
Print[" diagonala"]; 
Print[diag];
hdiag=Table[-p[xleft+i*h+h/2],{i,1,n}];
Print[" horna diagonala"]; 
Print[hdiag];
ddiag=Join[Table[-p[xleft+i*h+h/2],{i,1,n-1}],{-(p[xleft+(n+1)*h-h/2]+p[xleft+(n+1)*h+h/2])}];
Print[" dolna diagonala"]; 
Print[ddiag];
(*ddiag=Table[c,{n-1}];*)
A=SparseArray[{Band[{1,2}]->hdiag,Band[{1,1}]->diag,Band[{2,1}]->ddiag},{n+1,n+1}];Print[" matica"];
Print[" "];
Print[MatrixForm[A]];
(*MatrixForm[Normal[A]];*)
riesenie=LinearSolve[A,N[ps]];


Print[" prava strana"]; 
Print[ps];
xs=Table[xleft+(i-1)*h,{i,1,n+2}];
ys=Join[{Leftbv},Table[riesenie[[i]],{i,1,n+1}]];
graf=Table[{xs[[i]],ys[[i]]},{i,1,n+2}];
vysledky=Table[{xs[[i]],ys[[i]],pres[xs[[i]]],Abs[(ys[[i]]-pres[xs[[i]]])]},{i,1,n+2}];
)

OkrsamoadjNP[n,a,b,ya,yb,p,q,f,presne]

g1=Plot[presne[x],{x,a,b},PlotRange->All,PlotStyle->RGBColor[0,1,0]];

g2=ListPlot[graf,PlotRange->All,PlotStyle->RGBColor[0,0,1]];

Show[g1,g2]