A = (({{-52, 24, 0, 0}, {27, -52, 23, 0}, {0, 28, -52, 22}, {0, 0, 29, -52}})); b = {-26.8,-1.6,-2.4,-81.2839};
LinearSolve[A,b]

yt[x_]:=x+Exp[(x^2)]
a=0;b=1;h=0.2;n = ((b-a)/h);
xt = Table[a+i*h,{i,0,n}]

yt[xt]

n=10
A=Table[0,{i,0,n},{j,0,n}]

A//MatrixForm

Съставяме програмен код за решаване на гранична задача за линейно ДУ от II-ри ред

Clear[a,b,h,n,xt,yt,p,q,f,pt,qt,ft]
(* съставяне на мрежата *)
yt[x_]:=x+Exp[(x^2)] (* точното решение *)
a=0;b=1;h=0.2;n = ((b-a)/h);
xt = Table[a+i*h,{i,0,n}];
(* попълваме мрежовите функции p((Subscript[x, i])), q((Subscript[x, i])),f((Subscript[x, i])), i=0,1,...,n *)
p[x_]:=-2x
q[x_]:=-2
f[x_]:=-4x
pt = p[xt];
qt = Table[q[xt],{i,0,n}];
ft = f[xt];
(* задаваме коефициентите в граничните условия *)
alpha1=1;alpha2=1;
beta1=0;beta2=0;
gamma1=1;gamma2=1+E;
(* дефинираме структурата на матрицата А и вектора d на СЛАУ Ay=d *)
A=Table[0,{i,0,n},{j,0,n}];
d = ft;
d[[1]] = gamma1;d[[n+1]] = gamma2;
(* запълване ненулевите елементи на матрицата А *)
(* първи ред: *)
d
A[[1,1]] =alpha1-((3 beta1)/(2 h));
A[[1,2]] =((2 beta1)/( h));
A[[1,3]] =(( beta1)/(2 h));
(* последен ред: *)
A[[n+1,n-1]] =(( beta2)/(2 h));
A[[n+1,n]] =-((2 beta2)/( h));
A[[n+1,n+1]] =alpha2+(( 3beta2)/(2 h));
(* долен диагонал: *)
For[i=1, i≤ n-1,i++,
A[[i+1,i]]=(1/(h^2))-((pt[[i+1]])/(2h))
]
(* главен диагонал: *)
For[i=2, i≤ n,i++,
A[[i,i]]=-(2/(h^2))+qt[[i]]
]
(* горен диагонал: *)
For[i=3, i≤ n+1,i++,
A[[i-1,i]]=(1/(h^2))+((pt[[i-1]])/(2h))
]
A//MatrixForm
y=LinearSolve[A,d]
yt[xt]
points = Table[{xt[[i]],y[[i]]},{i,1,n+1}];
grp = ListPlot[points,PlotStyle->{Black,PointSize[0.02]}];
grfyt = Plot[yt[x],{x,a,b}, PlotStyle->{Red,Thick}];
Show[grfyt, grp]

Li