function [u,U] = forwardsolvrex(M,N,el,T,gmma,bta,pnodes,qnodes,pnodesh,qnodesh)
% usage: [u,U] = forwardsolvrex(M,N,el,T,gmma,bta,pnodes,qnodes,pnodesh,qnodesh)
% description: constructs and returns the solution to 
% our problem at time T, using an end Richardson extrapolation
% in both time and space.

% local variables: dx,xnodes,pnodes,qnodes,p1,U,u,Uh,uh,n,Ainv,unew

dx = el/M;
dt = T/N;
xnodes = linspace(dx,el,M)';
xnodesh = linspace(dx/2,el,2*M)';
% the initial value of u may *not* be zero if source term
% F(x,0) is not 0. Handle this separately. Of course, if 
% F(x,0)=0, this simply generates the zero solution for u.
u = zeros(M,1); % needed for U
U = zeros(M,1); % since U(x,0) = 0
u =sm(el,M,dt,0,bta,pnodes,qnodes)\rhs(el,M,dt,0,U,u,0,bta,pnodes,qnodes); 
uh = zeros(2*M,1); % needed for U
Uh = zeros(2*M,1); % since U(x,0) = 0 
uh = sm(el,2*M,dt/2,0,bta,pnodesh,qnodesh)\ ... 
        rhs(el,2*M,dt/2,0,Uh,uh,0,bta,pnodesh,qnodesh); 
% get the coefficient matrix A only once. To speed things
% up in *this* code, we invert A. In production code one 
% would obtain LU factorization only once, then repeatedly
% forward and backward solve the system.
Ainv = inv(sm(el,M,dt,gmma,bta,pnodes,qnodes));
Ainvh = inv(sm(el,2*M,dt/2,gmma,bta,pnodesh,qnodesh));
% now march to the requested solution
for n = 1:N
  % first the full step
  unew = Ainv*rhs(el,M,dt,n,U,u,gmma,bta,pnodes,qnodes);
  U = U/exp(gmma*dt) + dt/2*(u/exp(gmma*dt) + unew);
  u = unew;
  % then two half steps
  unew = Ainvh*rhs(el,2*M,dt/2,2*n-1,Uh,uh,gmma,bta,pnodesh,qnodesh);
  Uh = Uh/exp(gmma*dt/2) + dt/4*(uh/exp(gmma*dt/2) + unew);
  uh = unew;
  unew = Ainvh*rhs(el,2*M,dt/2,2*n,Uh,uh,gmma,bta,pnodesh,qnodesh);
  Uh = Uh/exp(gmma*dt/2) + dt/4*(uh/exp(gmma*dt/2) + unew);
  uh= unew; 
end;
% finally, extrapolation
u = 4/3*uh(2:2:2*M) - 1/3*u;
U = 4/3*Uh(2:2:2*M) - 1/3*U;