function biharmonic_cheby1d ( N )

%*****************************************************************************80
%
%% biharmonic_cheby1d() solves the 1D biharmonic equation by Chebyshev methods.
%
%  Discussion:
%
%    In the original distribution, this is p38.m.
%
%    Solve u_xxxx = exp(x), 
%    u(-1) = u(1) = u'(-1) = u'(1) = 0
%
%  Modified:
%
%    24 May 2021
%
%  Author:
%
%    Lloyd Trefethen
%
%  Reference:
%
%    Lloyd Trefethen,
%    Spectral methods in MATLAB,
%    SIAM, 2000,
%    ISBN: 0898714656,
%    LC: QA377.T65.
%
%  Input:
%
%    integer N: the number of grid points.
%
  if ( nargin < 1 )
    N = 15;
  end
%
%  Construct the spectral approximation to the discrete biharmonic operator:
%
  [ D, x ] = cheb ( N );
  S = diag ( [0; 1./(1-x(2:N).^2); 0] );
  D4 = ( diag(1-x.^2)*D^4 - 8*diag(x)*D^3 - 12*D^2 ) * S;
  D4 = D4(2:N,2:N);
%
%  Solve the boundary-value problem in the interior, 
%  and then pad solution with zero boundary values.
%
  f = exp ( x(2:N) ); 
  u = D4 \ f; 
  u = [ 0.0; u; 0.0 ];
%
%  Fit a polynomial to the U data.
%
  xx = (-1.0:.01:1.0)';
  uu = ( 1.0 - xx.^2 ) .* polyval ( polyfit(x,S*u,N), xx );  
%
%  Plot the computed result as dots:
%
  clf ( ); 
  hold ( 'on' );
    plot ( xx, uu, 'b-', 'linewidth', 2 );
    plot ( x, u, 'r.', 'markersize', 20 );
  hold ( 'off' );
  axis ( [ -1.0, 1.0, -0.01, 0.06 ] );
  grid ( 'on' );
  title ( 'biharmonic cheby1d' );
  filename = 'biharmonic_cheby1d.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Determine the exact solution, which has the form
%    u(x) = a + bx + cx^2 + dx^3 + exp(x)
%  Because u(x) must satisfy 4 boundary conditions, we
%  can determine the polynomial coefficients.
%
  A = [ 1, -1,  1, -1; ... 
        0,  1, -2,  3; ...
        1,  1,  1,  1; ...
        0,  1,  2,  3 ];
  V = vander ( xx ); 
  V = V(:,end:-1:end-3);
  c = A \ exp ( [ -1; -1; 1; 1 ] ); 
  exact = exp ( xx ) - V * c;
%
%  Compute the maximum approximation error:
%
  err = norm ( uu - exact, inf );
  fprintf ( 1, '  Maximum error = %g\n', err );

  return
end

