Solving ODE with Spectral Method using Chebyshev Polynomials

Question

I would like to solve following the basic equation of linear elasticity (for simplicity in 1D) $$ \frac{d}{dx} \left( E \frac{du}{dx} \right) = 0 \quad \textrm{with} \quad u(1)=0, \; u(-1)=b $$

Following the example on page 138 (program 33) in Nick Trefethen's book "Spectral Methods in Matlab" (https://people.maths.ox.ac.uk/trefethen/spectral.html), shouldn't be too complicated to adapt this to the above specified equation of elasticity. According to the boundary conditions, I removed the first column and row of the differentiation matrix. However, something went wrong with the implementation, since I the results is oscillating quite heavy (expecting a linear behavior). I would be appreciated for any suggestions!

clear all; clc;
N = 16; % Number of collocation points
[D,x] = cheb(N); % Set up differentiation matrix

% Insert boundary condition: u(-1) = 0, u(1) = 2
D2 = D(2:N+1,2:N+1);
D2(N,:) = zeros(1,N);
D2(N,N)=ones(1,1);

% right hand side
f = zeros(N,1);
f(N) = 2.0;
u = (D2.*(D2))\f;
u = [0;u];

% Plot solution
clf
%subplot ('position', [.1 .4 .8 .5])
plot(x,u,'.', 'markersize',16)
xx = -1:0.01:1;
uu = polyval(polyfit(x,u,N+1),xx);
line(xx,uu), grid on

with:

% CHEB compute D = differentitation matrix, x = Chebyshev grid

function [D, x] = cheb(N)
  if N == 0, D = 0; x = 1; return, end
  x = cos(pi*(0:N)/N)';
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';
  D = (c*(1./c)')./(dX+(eye(N+1)));
  D = D - diag(sum(D'));

Answer 1

Your code does not solve the BVP you posted. Here is the revised version that works well.

function bvp(N)

  [D,x] = cheb(N); % Set up differentiation matrix
  D2 = D^2;
  
  % Insert boundary condition: u(-1) = 0, u(1) = 2
  D2(1,:) = zeros(1,N+1);
  D2(1,1) = 1;
  D2(N,:) = zeros(1,N+1);
  D2(N,N) = 1;

  % right hand side
  f = zeros(N+1,1);
  f(1) = 2;
  u = D2\f;

  % Plot solution
  clf
  plot(x,u,'.', 'markersize',16)
  xx = -1:0.01:1;
  uu = polyval(polyfit(x,u,N+1),xx);
  line(xx,uu), grid on

end

Call it as bvp(N), where N is the number of Chebyshev points.