function retval = operatorA(M,N,el,T,gmma,bta,gnodes,pnodes,qnodes)
% usage: psi = operatorA(M,N,el,T,gmma,bta,gnodes,pnodes,qnodes)
% description: constructs and returns the value of operator A 
% on the function q(x), given as a column vector of discrete 
% values at the abscissa array xnodes = linspace(el/M,el,M).
% The function named "mu" should really be mu'(t) for this routine.
% gnodes is used to define the operator.  The gnodes vector
% gives values of g at nodes linspace(0,el,M+1).

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

dx = el/M;
dt = T/N;
xnodes = linspace(dx,el,M)';
u = zeros(M,1); % since u_t(x,0) = 0
U = zeros(M,1); % since U(x,0) = 0
% 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. Or use a band method.
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;
% we now have U at time T on xnodes
% calculate the numerator of returned value
retval = ((1+pnodes(1:M-1)*dx/2).*gnodes(1:M-1) - ...
         2*gnodes(2:M) + (1-pnodes(1:M-1)*dx/2).*gnodes(3:M+1));
% do the last entry by quadratic interpolation
retval = [retval;0]/(dx*dx);
retval = retval./U;
retval(M) = 3*retval(M-1)-3*retval(M-2)+retval(M-3);