p38.m

% p38.m - solve u_xxxx = exp(x), u(-1)=u(1)=u'(-1)=u'(1)=0
% (compare p13.m)
% Construct discrete biharmonic operator:
N = 15; [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 boundary-value problem and plot result:
f = exp(x(2:N)); u = D4\f; u = [0;u;0];
clf, subplot('position',[.1 .4 .8 .5])
plot(x,u,'.','markersize',16)
axis([-1 1 -.01 .06]), grid on
xx = (-1:.01:1)';
uu = (1-xx.^2).*polyval(polyfit(x,S*u,N),xx);
line(xx,uu,'linewidth',.8)
% Determine exact solution and print maximum error:
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;
title(['max err = ' num2str(norm(uu-exact,inf))],'fontsize',12)