function[xvect,Tmego]=vdm1(Nx)
%
%Véges differenciák módszere egy L hosszúságú rúd stacionárius
% hőeloszlásának kiszámítására. Ez a 
% 
%     d^2 T   
%     -----  + c (Tinf - T) = 0
%      dx^2
%típusú kétpontos peremérték feladatok megoldását igényli, amelyet a
%szokásos  differenciasémával approximálunk. 
%A peremfeltételeket közvetlenül beírhatjuk a sémába.
%Nx: a térbeli osztásrészek száma
Ta   = 300;      % bal oldali végponteli peremfeltétel
Tb   = 400;     % jobb oldali végponteli peremfeltétel
Tinf = 200;      % külső hőmérséklet
c    = 0.05;    % állandó
L      = 10;    % a rúd hossza
ndivs  = Nx;
nunknowns = ndivs - 1;
deltax = L/ndivs;
A = -(2 + deltax^2*c);
B = -deltax^2*c*Tinf;
T0=Ta; T1=Tb;
matrix = zeros(nunknowns);  % initialize

matrix(1,1) = A;                          % first row
matrix(1,2) = 1;
rhs(1)      = B - Ta;


for i=1:Nx+1
xvect(i)=(i-1)*deltax;
end



for i = 2:nunknowns - 1                   % intermediate rows
  matrix(i,i-1) = 1;
  matrix(i,i)   = A;
  matrix(i,i+1) = 1;
  rhs(i)        = B;
end;

matrix(nunknowns, nunknowns-1) = 1;       % last row
matrix(nunknowns, nunknowns)   = A;
rhs(nunknowns)                 = B - Tb;

T = matrix\rhs';                          % solve for unknowns

Tmego(1)               = Ta;              % reconstruct full vector
Tmego(2:1 + nunknowns) = T(:);
Tmego(nunknowns + 2)   = Tb;

lambda=sqrt(c); s1=exp(lambda*L); s2=1/s1; szam1=(T0-Tinf)*s2-(T1-Tinf);
nev=s2-s1;
szam2=-(T0-Tinf)*s1+(T1-Tinf);
A1=szam1/nev; B1=szam2/nev; 
Tpontos=pontosvdm(Nx,xvect,A1,B1);
hiba=norm(Tmego-Tpontos,inf);
disp('lépésköz és  hiba max. normában:'),deltax, hiba
i=1:Nx+1;
plot(xvect(i),Tmego,'r', xvect(i), Tpontos,'b')
xlabel('vezeték'), ylabel('véges differenciás módszer (piros), pontos megoldás (kék)')