%
% DATA ULOHY (uzivame vektorove operace, zapis operace "s teckou")
%
p = 0.25; 

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

%
% OBLAST:
%

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

%
% SIT, volba kroku, stabilita
%

n = 10;
maxk = 100; 
h = (b - a) / n;
tau = T / maxk;

sigma = p * tau / (h * h);

printf('Kriterium pro stabilitu  sigma <= 0.5: je %g <= 0.5 ? \n', sigma);
if (sigma > 0.5)
    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 * xx; % << nastavime rozmer pole pro reseni

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


for k = 1:maxk,
  for i = 2:(n-1)
    U(k + 1, i) = (1 - 2 * sigma) * U(k, i) + sigma * U(k, i + 1) + sigma * U(k, i - 1) + tau * f(k, i);
  end
end

