problem:
-(p(x)y')' + q(x)y = f(x)

priklad:
-((1/x^2)*y')'  +  (4/x^4)y - 6/x^3
y(1) = 0
y(2) = -1.5

presne riesenie:
y(x) = 1/x-x

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;
pres[x_] := (1/x) - x;

Needs["LinearAlgebra`Tridiagonal`"]

Okrajsamoadj[n_, xleft_, xright_, LeftBv_, RightBv_, f_, p_, q_, 
        pres_] := (h = 
          N[(xright - xleft)/(n + 1)]; 
pStrana = 
          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]*h^2 + (p[xright - h/2])*
                  RightBv}]; 
diag = 
          Table[p[xleft + i*h - h/2] + p[xleft + i*h + h/2] + 
              q[xleft + i*h]*h^2, {i, 1, n}]; 
hornaDiag = 
          Table[-p[xleft + i*h + h/2], {i, 1, 
              n - 1}]; 

A = 
          IdentityMatrix[
            n]; 
Do[{A[[i, i]] = 
              diag[[i]], 
A[[i + 1, i]] = 
              hornaDiag[[i]], 
A[[i, i + 1]] = 
              hornaDiag[[i]]}, {i, 1, 
            n - 1}]; 
A[[n, n]] = 
          diag[[n]]; 
Print[
          MatrixForm[A]]; 
riesenie = 
          TridiagonalSolve[N[hornaDiag], N[diag], N[hornaDiag], 
            N[pStrana]]; 
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}];)

Okrajsamoadj[n, a, b, ya, yb, f, p, q, pr]

g1 = 
      Plot[pres[x], {x, a, b}, PlotRange -> All, 
        DisplayFunction -> Identity]
g2 = 
      ListPlot[graf, 
        PlotStyle -> {RGBColor[0, 0.7, 0], PointSize[0.02]}, 
        DisplayFunction -> Identity]
Show[g1, g2, DisplayFunction -> $DisplayFunction]