%---------------------------------------------------------------------
function de=heat
global T t dt dx J x u un B C b dc sip sim theta
% T je konecny cas
% dt,dx jsou casovy a prostorovy krok
% J je pocet bunek
% x(j) je prostorova souradnice j-teho bodu
% u(j) je reseni v bode x(j) na casove vrstve n
% un(j) je reseni v bode x(j) na casove vrstve n+1
% resime na intervalu (-B,B)
% C je CFL cislo
% resime rovnici u_t = b u_xx
B = 20;
J = 200;
T = 10;
C = 0.45;
b = 2;
%figure;
solve;
%---------------------------------------------------------------------
function so=solve
% resi PDR u_t = a u_x
global T t dt dx J x u un B C a ermax

t=0;
init;
n=0;
while t < T
  central;
%  LW;
%  LF;
  plot(x,un);
%  axis([-B B 0 1]);
  drawnow;
 %pause;
  u = un;
  t = t + dt;
  n = n + 1;
end

n = n

%---------------------------------------------------------------------
function ii=init
% inicializace, pocatecni podminky
global T dt dx J x u un B C b

dx = 2*B/J;
x = zeros(J,1);
for j=1:J
  x(j) = -B + j*dx;
end
dt = C*dx^2/abs(b);
u = zeros(J,1);
un = zeros(J,1);
for j=1:J
%  u(j) = exp(-x(j)^2);
%  if abs(x(j)) < pi
  if x(j) < 0
%      u(j) = (1 + cos(x(j)))/2;
      u(j) = 1;
  end
end
%---------------------------------------------------------------------
function do=central
% jeden casovy krok  schema
global T dt dx J x u un B C b

for j=2:J-1
  un(j) = u(j) + b*dt*(u(j+1) - 2*u(j) + u(j-1))/dx^2;
end
% okrajove podminky
un(1) = un(2);
un(J) = un(J-1);
