Samoadjungovane okrajove ulohy

Riesme okrajovu ulohu ::  -(1/x^2y')'+4x^4y=6/x^3 ptresne riesenie::    y(x)=1/x-x

f[x_]:=-6/x^3
n=9;
a=1;
b=2;
ya=0;
yb=-1.5;
a0=-1;
p[x_]:=1/x^2;
q[x_]:=4/x^4;
presneriesenie[x_]:=1/x-x;

(*Algoritmus vypoctu*)

Okrsamoadj[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-1}],{f[xleft+n*h]-b*Rightbv}];
diag=Table[p[xleft+i*h-h/2]+p[xleft+i*h+h/2]+h^2*q[xleft+i*h],{i,1,n}];
hdiag=Table[-p[xleft+i*h+h/2],{i,1,n-1}];
A=SparseArray[{Band[{1,2}]->hdiag,Band[{1,1}]->diag,Band[{2,1}]->hdiag},{n,n}];
riesenie=LinearSolve[A,N[ps]];
xs=Table[xleft+(i-1)*h,{i,1,n+2}];
ys=Join[{Leftbv},Table[riesenie[[i]],{i,1,n}],{Rightbv}];
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}];
}

Okrajul[n,a,b,ya,yb,a0,a1,a2,f,presneriesenie]