оригиналната задача Пример 5 от стр. 96 - сменяме с общи смесени гранични условия, т.е. използване (Subscript[α, 1]), (Subscript[α, 2]),β1,β2,γ1,γ2

(* Mrejov metod *)
p[x_]:=x
q[x_]:=-2.
f[x_]:=30.
utochno[x_]:=16.x^2+1
a=0;b=1;n=5.;h=((b-a)/n);

(* предефинираме коефициентите в граничните условия: *)
α1 = 1.;
α2 = 1.;
β1 = 0.;
β2 = 0.;
γ1 = 1.;
γ2 = 17.;

(* добавяме и дясното гранично  условие като още един ред в матрицата *)
A=Table[0,{i,0,n},{j,0,n}];B=Table[0,{i,0,n}];

For[i=2,i≤ n,i++,x=a+(i-1) h;(Subscript[A, [[i,i-1]]])=(1/(h^2))-((p[x])/(2 h));(Subscript[A, [[i,i]]])=-(2/(h^2))+q[x];(Subscript[A, [[i,i+1]]])=(1/(h^2))+((p[x])/(2 h));(Subscript[B, [[i]]])=f[x]];

(* предефинираме коефициентите от първия ред на матрицата *)
(Subscript[A, [[1,1]]])=(α1-((3. β1)/(2 h)));(Subscript[A, [[1,2]]])=((2. β1)/h);(Subscript[A, [[1,3]]])=-((β1)/(2 h));(Subscript[B, [[1]]])=γ1;
(* предефинираме коефициентите от  последния ред на матрицата *)
(Subscript[A, [[n+1,n-1]]])=((β2)/(2 .h));(Subscript[A, [[n+1,n]]])=-((2. β2)/h);(Subscript[A, [[n+1,n+1]]])=(α2+((3. β2)/(2 h)));(Subscript[B, [[n+1]]])=γ2;


Print[MatrixForm[A]," * u = ",MatrixForm[B]]
Print["Приближеното решение е: ",u=LinearSolve[A,B]]
Print["Точното решение е: ",tochno=Table[utochno[x],{x,a,b,h}]//N]
Print["Истинската грешка между приближеното решение и точното решение е: ",Abs[u-tochno]]
Print[" Теоретичната грешка между приближеното решение и точното решение е: ",(h^2)//N]
x=Table[a+i h,{i,0,n}]
xu=Table[{(Subscript[x, [[i]]]),(Subscript[u, [[i]]])},{i,1,n+1}]
gr1=ListPlot[xu,PlotStyle->{PointSize[0.015],Hue[.7]}];
xtochno=Table[{(Subscript[x, [[i]]]),(Subscript[tochno, [[i]]])},{i,1,n+1}];
gr2=ListPlot[xtochno,PlotStyle->{PointSize[0.015],Hue[.9]}];
Show[gr1,gr2]

Решаваме задачата с n=20

(* Mrejov metod *)
p[x_]:=x
q[x_]:=-2.
f[x_]:=30.
utochno[x_]:=16.x^2+1
a=0;b=1;n=20.;h=((b-a)/n);

(* предефинираме коефициентите в граничните условия: *)
α1 = 1.;
α2 = 1.;
β1 = 0.;
β2 = 0.;
γ1 = 1.;
γ2 = 17.;

(* добавяме и дясното гранично  условие като още един ред в матрицата *)
A=Table[0,{i,0,n},{j,0,n}];B=Table[0,{i,0,n}];

For[i=2,i≤ n,i++,x=a+(i-1) h;(Subscript[A, [[i,i-1]]])=(1/(h^2))-((p[x])/(2 h));(Subscript[A, [[i,i]]])=-(2/(h^2))+q[x];(Subscript[A, [[i,i+1]]])=(1/(h^2))+((p[x])/(2 h));(Subscript[B, [[i]]])=f[x]];

(* предефинираме коефициентите от първия ред на матрицата *)
(Subscript[A, [[1,1]]])=(α1-((3. β1)/(2 h)));(Subscript[A, [[1,2]]])=((2. β1)/h);(Subscript[A, [[1,3]]])=-((β1)/(2 h));(Subscript[B, [[1]]])=γ1;
(* предефинираме коефициентите от  последния ред на матрицата *)
(Subscript[A, [[n+1,n-1]]])=((β2)/(2 .h));(Subscript[A, [[n+1,n]]])=-((2. β2)/h);(Subscript[A, [[n+1,n+1]]])=(α2+((3. β2)/(2 h)));(Subscript[B, [[n+1]]])=γ2;


Print[MatrixForm[A]," * u = ",MatrixForm[B]]
Print["Приближеното решение е: ",u=LinearSolve[A,B]]
Print["Точното решение е: ",tochno=Table[utochno[x],{x,a,b,h}]//N]
Print["Истинската грешка между приближеното решение и точното решение е: ",Abs[u-tochno]]
Print[" Теоретичната грешка между приближеното решение и точното решение е: ",(h^2)//N]
x=Table[a+i h,{i,0,n}]
xu=Table[{(Subscript[x, [[i]]]),(Subscript[u, [[i]]])},{i,1,n+1}]
gr1=ListPlot[xu,PlotStyle->{PointSize[0.015],Hue[.7]}];
xtochno=Table[{(Subscript[x, [[i]]]),(Subscript[tochno, [[i]]])},{i,1,n+1}];
gr2=ListPlot[xtochno,PlotStyle->{PointSize[0.015],Hue[.9]}];
Show[gr1,gr2]