%
% DATA ULOHY (uzivame vektorove operace, zapis operace "s teckou")
%
c2 = 4; 

fxt = '0 * x';
u0  = 'sin(pi * x)';
u1  = '0';
alfa = '0';
beta = '0';

%
% OBLAST:
%

b = 1; a = 0; 
T = 3; 

%
% SIT, volba kroku, stabilita
%

n = 15;
maxk = 80; 
h = (b - a) / n;
tau = T / maxk;

sigma2 = c2 * (tau * tau) / (h * h);

printf('sigma2 = %g \n', sigma2);
%
% SIT
%
xx = a:h:b; tt = 0:tau:T;

[xsit, tsit] = meshgrid(xx, tt);

x = xsit; t = tsit; f = eval(fxt);

%
% pocatecni a okrajove podminky
%

U = 0 * xsit; % << nastavime rozmer pole pro reseni       (poloha)

x = xx; 
U(1,:) = eval(u0); 
U(2,:) = U(1,:) + tau * eval(u1); 

t = tt; U(:,1) = eval(alfa); U(:,end) = eval(beta);



% Matice prechodu z jedne casove vrstvy na dalsi:

%%   d2 u /dt2 = c2 d2 u /dx2 + f(x,t);
%%   (u2 - 2 u1 + u0 ) / tau2 = c2 (u2(i-1) - 2 u21i +  u2(i+1)) / ( 2 h2 ) + c2 (u0(i-1) - 2 u01i +  u0(i+1)) / ( 2 h2 ) + f (x,t)
%%

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

A(1,1:2)      = [ 1, 0];
A(end,end-1:end)  = [0, 1];

[L1,U1] = lu(A);

for k = 2:maxk,
  rhs = (2 * U(k,:)' - A * U(k -1,:)' + tau * tau * f(k,:)');           
  rhs(1)   = U(k+1,1);
  rhs(end) = U(k+1,n + 1);
  newU = U1 \ (L1 \ rhs); 
  U(k+1,:) = newU';
end



