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

okrsamoadj[n_,xleft_,xrigth_,Leftbv_,Rigthbv_,p_,q_,f_,pres_]:=
(
h=N[(xrigth-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*x]*h^2+p[xrigth-h/2]*Rigthbv}];
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}];
Print[" matica sustavy:"];
Print[" "];
Print[MatrixForm[A]];
riesenie=LinearSolve[A,N[ps];
)