%---------------------------------------------------------------------
function de=pde
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 = 40;
J = 400;
T = 40;
C = 0.8;
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

t=0;
init;
du = u;
while t < T
%  lfstep;
%  lwstep;
%  central;
  forward;
  plot(x,un, 'LineWidth', 3);
  axis([-B B 0 1]);
  drawnow;
  u = un;
  if (t==0), pause, end
  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)/4) < pi
      u(j) = (1 + cos(x(j)/4))/2;
%      u(j) = 1;
  end
end
%---------------------------------------------------------------------
function lf=lfstep
% jeden casovy krok Lax-Friendrichsovo schema
global T dt dx J x u un B C a

for j=2:J-1
  un(j) = (u(j+1) + u(j-1))/2 - a*dt*(u(j+1) - u(j-1))/(2*dx);
end
% okrajove podminky
un(1) = u(2);
un(J) = u(J-1);
%---------------------------------------------------------------------
function lw=lwstep
% jeden casovy krok Lax-Wendroff schema
global T dt dx J x u un B C a

for j=2:J-1
  un(j) = u(j) - a*dt*(u(j+1) - u(j-1))/(2*dx) ...
               + a^2*dt^2/(2*dx^2)*(u(j+1)-2*u(j)+u(j-1));
end
% okrajove podminky
un(1) = u(2);
un(J) = u(J-1);
%---------------------------------------------------------------------
function ce=central
% jeden casovy krok Lax-Friendrichsovo schema
global T dt dx J x u un B C a

for j=2:J-1
  un(j) = u(j) - a*dt*(u(j+1) - u(j-1))/(2*dx);
end
% okrajove podminky
un(1) = u(2);
un(J) = u(J-1);
%---------------------------------------------------------------------
function fo=forward
% jeden casovy krok Lax-Friendrichsovo schema
global T dt dx J x u un B C a

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