function [u,U] = forwardsolv(M,N,el,T,gmma,bta,pnodes,qnodes)
% usage: [u,U] = forwardsolv(M,N,el,T,gmma,bta,pnodes,qnodes)
% description: constructs and returns the solution to 
% our problem at time T.

% local variables: dx,xnodes,n,Ainv,unew

dx = el/M;
dt = T/N;
xnodes = linspace(dx,el,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, next lines simply generates the zero solution for u.
u = zeros(M,1); % since u(x,0) = 0
U = zeros(M,1); % since U(x,0) = 0
u = sm(el,M,dt,gmma,bta,pnodes,qnodes)\rhs(el,M,dt,0,U,u,0,bta,pnodes,qnodes);
% 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));
% now march to the requested solution
for n = 1:N
  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;
end;