function q = finalp1_FEMsolve(x,dt,T)

%  function q = finalp1_FEMsolve(x,dt,T)
%  Use FEM with vector of N+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 + dq/dx = 0, 0<x<1, 0<t
%        q(x(1:nx),0) = 0, q(0,t) = sin(50t).
%  
%  and return the solution vector q(x(1:N),T) at time T.
  
  nt = round(T/dt);
  nx = length(x) - 1;
  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) = sin(50*t(n)).

  %  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),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]; % Automatically handles right boundary element.
  Ijj   = 2*(Ijjm1 + Ijjp1);
  I = (spdiags([Ijjp1 Ijj Ijjm1],-1:1,nx,nx))'; % Specify diagonals of I to
                                                % start all on row 1.
                                                
  Jjjm1 = -0.5;
  Jjjp1 = 0.5;
  e = ones(nx,1);
  J = (spdiags([Jjjp1*e 0*e Jjjm1*e],-1:1,nx,nx))'; % ...as for I
  J(nx,nx) = 0.5; % Nonzero for right boundary element 
                  % since no cancelling x > x(nx) contribution there.

  M = I/dt +0.5*J;

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

  % Timestepping loop

  q = zeros(nx,1); % Initial condition

  for np1 = 1:nt
    tnp1 = np1*dt;
    q0np1 = sin(50*tnp1);         % Left BC at new time
    q0n   = sin(50*(tnp1-dt));
    rhs =   (I/dt -0.5*J)*q;
    r10  = -I10*(q0np1-q0n)/dt - 0.5*J10*(q0n+q0np1);
    rhs(1) = rhs(1)+ r10;
    
    q = M\rhs; % Update q to new time level
  end