%---------------------------------------------------------------------
function de=adv_cn
global T t dt dx J x u un B C a 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 + a u_x = 0
B = 20;
J = 800;
T = 5;
C = 2;
a = 1;
solve;
%---------------------------------------------------------------------
function so=solve
% resi PDR u_t = a u_x
global T t dt dx J x u un B C a M

t=0;
init;
M = zeros(J,J);
M(1,1) = 1;
M(1,2) = -1;
%M(1,J-1) = -1;
for j=2:J-1
  M(j,j-1) = - a*dt/(4*dx);
  M(j,j)   = 1;
  M(j,j+1) =   a*dt/(4*dx);
end
M(J,J-1) = -1;
M(J,J) = 1;
%M(J,2) = -1;
M = inv(M);
while t < T
  imstep;
  plot(x,un);
%  axis([-B B 0 1]);
  drawnow;
  %pause
  u = un;
  t = t + dt;
end
%---------------------------------------------------------------------
function ii=init
% inicializace, pocatecni podminky
global T dt dx J x u un B C a

dx = 2*B/J;
x = zeros(J,1);
for j=1:J
  x(j) = -B + j*dx;
end
dt = C*dx/abs(a);
u = zeros(J,1);
un = zeros(J,1);
for j=1:J
  if abs(x(j)) < 10 %pi
      %u(j) = (1 + cos(x(j)))/2;
      u(j) = 1;
  end
end
%---------------------------------------------------------------------
function im=imstep
% jeden casovy krok implicitni schema
global T dt dx J x u un B C a M

uu = u;
uu(1) = 0;
uu(J) = 0;
for j=2:J-1
  uu(j) = u(j) - a*dt/(4*dx)*(u(j+1) - u(j-1));
end
un = M*uu;
%---------------------------------------------------------------------


