function [qT,qr] = advect_FEMsolve(x,c,dt,q0,ql)

%  function q = advect_FEMsolve(x,c,dt,q0,ql)
%  Use FEM with vector of nx+1 monotonic increasing nodes x (first node x(0)=0,
%  last node x(nx)=1) and trapezoidal time-differencing with timestep dt
%  to solve:
%        dq/dt + c(x)*dq/dx = 0, 0<x<1, 0<t
%        q(x(1:nx),0) = q0(x), 
%        q(0,t) = ql(t)
%  Inputs:
%    x  Row vector of nodes (including x=0) of length nx+1
%    c  If c is a scalar then it is the constant advection speed. 
%       If c is a vector of the same length as x, c is the advection
%       speed at the vector of positions x.
%    dt timestep
%    q0 (column vector of length nx, at the nodes, excluding x=0)
%    ql (specified at times n*dt, n=0,1,.,nt.)
%  Returns:
%    qT   Solution column vector q(x(1:N),T) at time T=dt*nt. 
%    qr  Solution q(1,t) at the right boundary for t = 0,dt,...,nt*dt   
  
  nt = length(ql)-1; % Number of timesteps to take
  nx = length(x) - 1;  % Number of unknowns
  dx = diff(x)'; % Turns dx into a column vector, 

  %  Trapezoidally differenced FEM is:
  %  (   (I(j,j-1)*q(j-1,n+1) + I(j,j)*q(j,n+1) + I*(j,j+1)*q(j+1,n+1)) 
  %    - (I(j,j-1)*q(j-1,n  ) + I(j,j)*q(j,n  ) + I*(j,j+1)*q(j+1,n  )))/dt
  %  +  0.5*( J(j,j-1)*q(j-1,n+1) + J(j,j)*q(j,n+1) + J(j,j+1)*q(j+1,n+1)
  %          +J(j,j-1)*q(j-1,n  ) + J(j,j)*q(j,n  ) + J(j,j+1)*q(j+1,n)) 
  %  for j = 1,...,nx-1 with  q(0,n) = ql(n+1)

  %  We set up the FEM in matrix notation for simplicity:
  %
  %    I*(q(:,n+1)-q(:,n))/dt + 0.5*J*(q(:,n+1)+q(:,n))= r10 (j=1)
  %                                                      0   (j>1)
  %
  %  Here I is the matrix of basis fn inner products I(j,n) = <phij(x),phin(x)>
  %  and  J is the matrix of inner products J(j,n) = <phij(x),c*dphin(x)/dx>
  %  where j,n = 1,...,N. The right hand side comes from the j-1=0
  %  terms, which are known from the left BC:
  %  r10  = -I(1,0)*(q0(n+1)-q0(n))/dt - 0.5*J(1,0)*(q0(n+1)+q0(n));
    
  Ijjm1 = dx/6;
  Ijjp1 = [Ijjm1(2:nx); 0]; % Last entry handles right boundary element j=nx,
                            % which has no right neighbor.
  Ijj   = 2*(Ijjm1 + Ijjp1);
  I = (spdiags([Ijjp1 Ijj Ijjm1],-1:1,nx,nx))'; % Specify diagonals of I to
                                                % start all on row 1.

  if (length(c)==1)
    c = c*ones(1,nx+1);
  end
  
  % Note c(j+1) corresponds to index j
  
  cjm1 = c(1:nx);
  cj = c(2:nx+1);
  cjp1 = [c(3:nx+1) NaN];
  
  Jjjm1 = -(cjm1/6 + cj/3)';
  Jjjp1 = (cj/3 + cjp1/6)'; % Note J_nx,nx+1 is not used
  Jjjp1(nx) = 0; % Takes care of right boundary
  Jjj   = -(Jjjm1 + Jjjp1);
  e = ones(nx,1);
  J = (spdiags([Jjjp1 Jjj Jjjm1],-1:1,nx,nx))'; % ...as for I

  M = I/dt +0.5*J;

  I10 = Ijjm1(1); % Note i0, j0 each have only one nonzero element I10, J10
  J10 = Jjjm1(1); % This is handled without vector notation, since easier.

  % Timestepping loop

  q = q0; % Initial condition
  qr = zeros(1,nt+1); % Solution at right boundary to be output each timestep
  qr(1) = q0(nx);  % Initial condition at right boundary
  for np1 = 1:nt
    %  Note ql(np1) = ql at time level n= np1-1, so that ql(1) = ql at t=0.
    tnp1 = np1*dt;
    rhs =   (I/dt -0.5*J)*q;
    r10  = -I10*(ql(np1+1)-ql(np1))/dt - 0.5*J10*(ql(np1)+ql(np1+1)); % left BC
    rhs(1) = rhs(1)+ r10;
    
    q = M\rhs; % Update q to new time level
    qr(np1+1) = q(nx); % Store q at right boundary
  end
  qT = q;
  
