function retval = sm(el,M,dt,gmma,bta,pnodes,qnodes)
% usage: A = sm(el,M,dt,gmma,bta,pnodes,qnodes)
% description: returns the system matrix for this problem.

% local variables: dx, ii

retval = zeros(M,M);
dx = el/M;
pnodes = dx/2*pnodes;
qnodes = dx*dx*qnodes;
% do first row separately inserting Dirichlet condition
retval(1,1) = 2 + qnodes(1)*(1-gmma*dt/2);
retval(1,2) = -1 + pnodes(1);
% now the interior rows
for ii = 2:M-1
  retval(ii,ii-1) = -1 - pnodes(ii);
  retval(ii,ii) = 2 + qnodes(ii)*(1-gmma*dt/2);
  retval(ii,ii+1) = -1 + pnodes(ii);
end;
% fix the last row by centered differences, etc.
retval(M,M-1) = -2;
retval(M,M) = 2 + qnodes(M)*(1-gmma*dt/2) + (1-pnodes(M))*2*dx/bta;
% or do last row by a one-sided derivative formula (more accurate)
%retval(M,M-2) = bta;
%retval(M,M-1) = -4*bta;
%retval(M,M) = 2*dx + 3*bta;