function [Q,xi] = wave(BC,ax,bx,n,cfl,q0,out,lim,T,A)
% wave(BC,ax,bx,n,cfl,q0,out,lim,T,u) für q_t + A q_x = 0
%INPUT:
%       BC       Randbedingungen: 0=Reflektierend,1=Periodisch
%       ax, bx   Ränder
%       n        Anz. Zellen in x exklusive ghost cells
%       T        Endzeitpunkt
%       cfl      CFL-Zahl
%       out      Erzeuge Plot nach je "out" schritten. 
%       lim      'superbee'/'minmod'/'MC'/'van_leer'
%       q0       Funkion für AD.
%       A        Matrix in Gleichung

dx = (bx-ax)/n;
xi = linspace(ax+dx/2,bx-dx/2,n);
x = linspace(ax-2*dx+dx/2,bx+2*dx-dx/2,n+4);  % je 2 ghost cells links und rechts

t = 0;

J = 3:n+2;

Q = q0(x);

if BC==0    % refl. Randbed.
  M = [1 0; 0 -1];
  Q(:,1) = M*Q(:,4); Q(:,end)   = M*Q(:,end-3);
  Q(:,2) = M*Q(:,3); Q(:,end-1) = M*Q(:,end-2);
else        % periodische Randbed.
  Q(:,1) = Q(:,end-3); Q(:,end)   = Q(:,4);
  Q(:,2) = Q(:,end-2); Q(:,end-1) = Q(:,3);
end
k=1;

% Anfangsdaten plotten
if out~=0
  subplot(2,1,1);
  g = plot(xi,Q(1,J));
  axis([ax bx 0 1]);
  title('p');

  subplot(2,1,2);
  h = plot(xi,Q(2,J));
  axis([ax bx -1 1]);
  title('u');
end

F = zeros(length(A),n+2);

while t<T
    % Löse Riemann-Problem
    [W,lambda,alpha,AmDQ,ApDQ] = rp_ac_1d(Q(:,1:end-1),Q(:,2:end),A) ;
    
    dt = cfl*dx/max(abs(diag(lambda)));  % Zeitschrittbeschränkung
    
    k=k+1;
    if t+dt>T       % letzter Schritt ist evtl. zu groß
      dt = T-t;
    end
    
    % Upwind anwenden
    Q(:,J) = Q(:,J)-dt/dx*(ApDQ(:,2:end-2) + AmDQ(:,3:end-1));
    
    % Bestimme Anstiegslimiter
    Te  = tet(alpha,diag(lambda));
    phi = limiter(Te,lim);
    
    % Berechne Korrektur
    for j = 2:n+2
      F(:,j) = W(:,:,j)*(0.5*abs(diag(lambda)).*(1-dt/dx*abs(diag(lambda))).*phi(:,j-1));
    end
    
    % Korrektur anwenden
    Q(:,J) = Q(:,J) - dt/dx*(F(:,3:end) - F(:,2:end-1));
    
    t=t+dt;
    
    if BC==0    % refl. Randbed.
      M = [1 0; 0 -1];
      Q(:,1) = M*Q(:,4); Q(:,end)   = M*Q(:,end-3);
      Q(:,2) = M*Q(:,3); Q(:,end-1) = M*Q(:,end-2);
    else        % periodische Randbed.
      Q(:,1) = Q(:,end-3); Q(:,end)   = Q(:,4);
      Q(:,2) = Q(:,end-2); Q(:,end-1) = Q(:,3);
    end
    
    % Ausgabe
    if mod(k,out)==0
        set(g,'ydata',Q(1,J));
        set(h,'ydata',Q(2,J));
        pause(0.3*dt);
        %input('Hit <return> to continue  ');
    end
end