% p38_mymod.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);
  D4 = D^4;
  D4 = [D(1,2:N);D4(3:(N-1),2:N);D(N+1,2:N)];
  
% Solve boundary-value problem and plot result:
  f = exp(x(2:N));
  f = [0; exp(x(3:N-1)); 0];
  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); 
  uu = polyval(polyfit(x,u,N),xx);
  line(xx,uu)

% 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)