function q1 = hw4_fv(q,q0,x,t,method)
  
% function q1 = hw4_fv(q,q0,x,t,method) 
% Solve the traffic flow problem on gridpoints x with
% discretization times t given left BC q0(t). Output q1 = q(xmax,t)

  nx = length(x);
  nt = length(t);
  dx = x(2)-x(1);
  dt = t(2)-t(1);
  
  q1 = zeros(1,nt); % Initialize q(1,t)
  
  
  
  % Implement FV method with mean value q(j) and general slope sig(j)
  % in each cell j = 1:nx. Here cell 1 corresponds to x = 0 (the BC)
  % and cell N corresponds to x = 1 (the outflow
  % boundary). The fluxes are f(j+0.5), j = 1, nx
  
  for it = 2:nt
    qh = 0.5*(q(1:nx)+q([2:nx nx]));
    ch = 1 - 2*qh;
    muh = ch*dt/dx;
    if (method == 'L')
      % Downwind slope like Lax-Wendroff.  At x=1,
      % assume downwind value is the same as at x=1.
      
      sig = [q(2:nx)-q(1:(nx-1)) q(nx)-q(nx)]/dx;
    else
      % method = 'M' - MC method
      
      sig1 = [q(2)-q(1) q(3:nx)-q(1:(nx-2)) q(nx)-q(nx-1)]/(2*dx);
      sig2 = 2*[0 q(2:nx)-q(1:(nx-1))]/dx;
      sig3 = 2*[q(2:nx)-q(1:(nx-1)) 0]/dx;
      sigsum = sign(sig1) + sign(sig2) + sign(sig3);
      allpos = sigsum>2;
      allneg = sigsum<-2;
      posneg = abs(sigsum)<2;
      sig = zeros(1,nx);
      if(length(sig1(allpos)) > 0) 
	sig(allpos) = min(sig1(allpos),sig2(allpos));
	sig(allpos) = min(sig(allpos),sig3(allpos));
      end
      if(length(sig1(allneg)) > 0) 
	sig(allneg) = max(sig1(allneg),sig2(allneg));
	sig(allneg) = max(sig(allneg),sig3(allneg));
      end
      if(length(sig1(posneg)) > 0) 
        sig(posneg) = 0;
      end
    end
    fh= q.*(1-q) + 0.5*dx*ch.*sig.*(1 - muh); 
    q(1) = q0(it); % BC at x = 0, t = (it-1)*dt
    q(2:nx) = q(2:nx) - (dt/dx)*(fh(2:nx) - fh(1:(nx-1)));

    q1(it) = q(nx);
  end