Neumanova podmienka...
-(xy')' + xy = x^2 - 2xsin x - cos x -1
y(0)=0
y'(π)=0
n - pocet vnutornych uzlov

f[x_] := x^2 + 2*x*Sin[x] - Cos[x] - 1;
n = 9;
a = 0;
b = π;
ya = 0;
yb = 0;
p[x_] := x
q[x_] := x
pres[x_] := Sin[x] + x

Needs["LinearAlgebra`Tridiagonal`"]

Okrajsamoadjnp[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}], {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}]; 
hornaDiag = 
          Table[-p[xleft + i*h + h/2], {i, 1, 
              n}]; 
dolnaDiag = 
          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])}]; 

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

Okrajsamoadjnp[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]