function [x u uexact e] = bvp1d(level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solves 1-d linear boundary value problem
%
%       u''(x) = f(x)  on  x = [0,1]; u(0) = u_L, u(1) = u_R
%
% using second-order finite difference techniques and sparse matrix
% operations.
%
% Input arguments
%
%    level (integer scalar): Disretization level: nx = 2^level + 1
%
% Output arguments
%
%    x (real vector): Row vector of length nx containing discrete x coords.
%    u (real vector): Row vector of length nx containing approximate soln.
%    e (real vector): Row vector of length nx containing soln error.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   levelmin = 3;
   levelmax = 20;

   xmin = 0.0;
   xmax = 1.0;

   omega = 12*pi;

   trace = 0;
   if trace 
      fprintf('bvp1d: In function level=%d\n', level);
   end

   % Default values for output arguments.
   x = zeros(1,0);
   u = zeros(1,0);
   e = zeros(1,0);

   % Input argument checking 
   err = 0;
   if level ~= fix(level) 
      fprintf('bvp1d: level=%g but should be an integer\n');
      err = 1;
   end 
   if level < levelmin | level > levelmax
      fprintf('bvp1d: level=%d must be in range %d .. %d\n', level, ...
         levelmin, levelmax);
      err = 1
   end
   if err
      return
   end

   % Define mesh and associated parameters
   nx = 2^level + 1;
   x = linspace(xmin, xmax, nx);
   delx = x(2) - x(1);

   % Compute exact solution and right hand side of equations
   uexact = sin(omega * x);
   b = (-omega^2) * uexact';

   % Define 3 diagonals for interior equations ...
   lower = (1 / delx^2) * ones(nx,1);
   main  = -2 * lower;
   upper = lower;

   % ... and adjust end values to implement boundary conditions.
   upper(2)    = 0;
   main(1)     = 1;
   b(1)        = uexact(1);

   main(nx)    = 1;
   lower(nx-1) = 0;
   b(nx)       = uexact(nx);

   if trace
      fprintf('bvp1d: lower diagonal\n');
      lower
      fprintf('bvp1d: main diagonal\n');
      main
      fprintf('bvp1d: upper diagonal\n');
      upper
      fprintf('bvp1d: rhs\n');
      b
   end

   % Set up triagonal system ...
   A = spdiags([lower main upper], -1:1, nx, nx);

   if trace
      fprintf('bvp1d: compressed tridiagonal matrix\n');
      A
      fprintf('bvp1d: full tridiagonal matrix\n');
      full(A)
   end

   % ... solve it, and take transpose so u is a row vector
   u = (A\b)';

   % Compute solution error
   e = uexact - u;

   if trace
      fprintf('bvp1d: exact soln\n');
      uexact
      fprintf('bvp1d: computed soln\n');
      u
      fprintf('bvp1d: soln error\n');
      e
   end

end
