NMVaje3

%{
-2*y1 + y2 = h^2
  y1 - 2y2 + y3 = h2
......................

      y(n-1) - 2yn = h^2
%}

%x = linspace(0,1);
n = 100;
x = linspace(0,1,n);
y = (x.*x)/2 - x/2;

plot(x,y, 'b');
hold on

A = diag(-2*ones(n,1)) + diag(ones(n-1,1),1) + diag(ones(n-1,1),-1);

h = x(2) - x(1);

b = ones(n) * h^2;

y2 = A\b;

plot(x, y, 'r');

enke = ones(n-1,1);
dvojke = -2*ones(n,1);

M=[[0;enke],dvojke,[enke;0]];

y3 = tridiag(M, b);

%%%%%%%%%%%%%%%


funkcija


function y = tridiag (M, d);
%{
Resi tridiagonalen sistem A*y=b
kjer n x n matriko A predstavimo y n x 3 matriko M
%}

n=length(d);

for i=1:n-1
  k = -M(i+1,1)/M(i,2);
  M(i+1,1)=0;
  M(i+1,2) = M(i+1,2)+ k*M(i,3);
  d(i+1) = d(i+1)+k*d(i);
endfor

y=zeros(n,1);
y(n)=d(n)/M(n,2);

for i=n:-1:2
  y(i-1)=(d(i-1)-M(i-1,3)*y(i))/M(i-1,2);
endfor


endfunction