Budeme riesit rovnice:
a2y'' + a1y' + a0y = f[x]
y[a]=ya      y[b]=yb       DIRICHLETOVE PODMIENKY
a0, a1, a2 ∈ R
PRIKLAD:
y'' - y = x^3 + 6x^2 + 1
y[0]=1    y[2]=5

pr[x_] := 
      8.286134462  Exp[x] + 5.713865537  Exp[-x] - x^3 - 6  x^2 - 6  x - 
        13
f[x_] := x^3 + 6  x^2 + 1
a2 = 1;
a1 = 0;
a0 = -1;
a = 0;
b = 2;
ya = 1;
yb = 5;

Needs["LinearAlgebra`Tridiagonal`"]

n - pocet krokov, xleft - zaciatok, xright - koniec, LeftBu - finkcna hodnota v xleft, RightBu - funkcna hodnota v xright
h - krok,   k, l, m - cisla na diagonalach

OkrajUl[n_, xleft_, xright_, LeftBv_, RightBv_, a0_, a1_, a2_, f_, 
        pres_] := Module[{h, a, b, c, pStrana, diag, hornaDiag, dolnaDiag, 
          riesenie}, 
        h = N[(xright - xleft)/(n + 1)]; 
a = 
          N[a0 - ((2*a2)/(h^2))]; 
b = 
          N[(a2/h^2) + (a1/(2*h))]; 
c = 
          N[(a2/h^2) - (a1/(2*h))]; 
pStrana = 
          Join[{f[xleft + h] - c*LeftBv}, 
            Table[f[xleft + (i*h)], {i, 2, n - 1}], {f[xleft + n*h] - 
                b*RightBv}]; 
diag = 
          Table[a, {n}]; 
hornaDiag = 
          Table[b, {n - 1}]; 
dolnaDiag = 
          Table[c, {n - 1}]; 
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}], {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[5, 0, 2, 1, 5, -1, 0, 1, f, pr]

g1 = 
      Plot[pr[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]