function lam = triop(m,flag)
% for a discretized differential operator: build, apply, solve system, and find eigenvalues
% examples:  >> tic, triop(4000); toc      % 0.004 secs
%            >> tic, triop(1e5); toc       % 0.04 secs
%            >> tic, triop(1e7); toc       % 5.2 secs ... A*b and A\b scale very well
%            >> tic, triop(40,'i'); toc    % 0.005 secs
%            >> tic, triop(1000,'i'); toc  % 3.5 secs
%            >> tic, triop(4000,'i'); toc  % 240 secs ... inv() scales *really* badly
%            >> tic, triop(40,'e'); toc    % 0.009 secs
%            >> tic, triop(1000,'e'); toc  % 0.5 secs
%            >> tic, triop(4000,'e'); toc  % 38 secs ... eig() scales badly

x = (1:m) / (m+1);              % points in the unit interval, except 0 and 1, ...
dx = 1.0 / (m+1);               % with this spacing

d = 2 + dx * dx * 9 * cos(x);   % on the diagonal
ab = - ones(1,m-1);             % super and sub diagonal
A = spdiag(ab,-1) + spdiag(d,0) + spdiag(ab,1); % use sparse storage

b = dx * dx * x';
c = A * b;
v = A \ b;
[norm(b,'inf') norm(c,'inf')  norm(v,'inf')]

if (nargin > 1) && (flag == 'i')
  B = inv(A);                   % "inv" is not suited for large m (certainly not m=10^5)
end

if (nargin > 1) && (flag == 'e')
  lam = sort(eig(A));           % "eig" is not suited for large m (certainly not m=10^5)
                                % Octave does not have "eigs", but Matlab does
  [lam(1:4)' lam(end-3:end)']
end