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

fxt = 'x .* (1 -x) .* exp(-t)';
u0  = 'sin(pi * x)';
u1  = '0';
alfa = '0';
beta = '0';

%
% OBLAST:
%

b = 1; a = 0; 
T = 5; 

%
% SIT, volba kroku, stabilita
%

n = 51;
maxk = 800; 
h = (b - a) / n;
tau = T / maxk;

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

printf('Kriterium pro stabilitu  sigma2 <= 1: je %g <= 1 ? \n', sigma2);
if (sigma2 > 1)
    printf('\nPOZOR, NESPLNENO!\n\nVypocet nestabilni metodou!\n');
else
    printf('Splneno.\n');
end
%
% 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);


%%
%%   d2 u /dt2 = c2 d2 u /dx2 + f(x,t);
%%   (u2 - 2 u1 + u0 ) / tau2 = c2 (u1(i-1) - 2 u1i +  u1(i+1)) / h2 + f (x,t)
%%
for k = 2:maxk,
  for i = 2:(n-1)
    U(k + 1, i) = (2 - 2 * sigma2) * U(k, i) + sigma2 * U(k, i + 1) + sigma2 * U(k, i - 1) - U(k - 1, i) + tau * tau * f(k, i);
  end
end

