Matematické Fórum

Zdravím, snažím se naprogramovat metodu konečných prvků ve Wolframu. Řeším diferencialni rovnici s okrajovými podmínkami a na intervalu .  Postupuju podle tohoto textu. A můj kód vypadá takto:


Clear[x, f, p, q];
n = 20; xt = Table[1/n i, {i, 1, n - 1}];
h = 1/n;
f[y_] := Pi^2 Cos[Pi y];
p[x_] := -1;
q[x_] := 0;
a0 = 1; al = -1; b0 = 0; bl = 0;
K = Table[0, {i, 1, n + 1}, {j, 1, n + 1}];
K[[1, 1]] = 1 + a0;
K[[2, 2]] = 1;
K[[1, 2]] = -1; K[[2, 1]] = -1;
K[[n + 1, n + 1]] = 1 + al;
K[[n, n]] = K[[n, n]] + 1; K[[n + 1, n]] = -1;
K[[n, n + 1]] = -1;
Do[K[[i, i]] = K[[i, i]] + 1 + q[xt[[i - 1]]];
  K[[i + 1, i + 1]] = K[[i + 1, i + 1]] + 1 + q[xt[[i]]];
  K[[i + 1, i]] = K[[i + 1, i]] - 1;
  K[[i, i + 1]] = K[[i, i + 1]] - 1, {i, 2, n - 1}];

F = Table[0, {i, 1, n + 1}];
F[[1]] = h/2 f[0] + b0;
F[[2]] = h/2 f[xt[[1]]];
F[[n + 1]] = h/2 f[1] + bl;
F[[n]] = F[[n]] + h/2 f[xt[[n - 1]]];
Do[F[[i]] = F[[i]] + h/2 f[xt[[i - 1]]];
  F[[i + 1]] = F[[i + 1]] + h/2 f[xt[[i]]], {i, 2, n - 1}];
F = F[[2 ;; n + 1]];
K = K[[2 ;; n + 1, All]];
F = F - 1*K[[All, 1]];
K = K[[All, 2 ;; n + 1]];
F = F[[1 ;; n - 1]];
K = K[[1 ;; n - 1, All]];
F = (F + 1*K[[All, n]]);
K = K[[All, 1 ;; n - 1]];
y = N[LinearSolve[K, F]];
y = Prepend[y, 1]; y = Append[y, -1];
xt = Prepend[xt, 0]; xt = Append[xt, 1];
ListLinePlot[Partition[Riffle[xt, y], 2]]

Výsledek nevypadá ani trochu jako kosinus. Sedím nad tím už dost dlouho a jsem přesvědčen, že tam bude špatně jenom nějaká blbost, za každou radu budu vděčný, děkuji.