% Legendre and cheb collocation for -u''+x u = f, u(pm 1)=0
n=15; n1=n+1; D=zeros(n1);
method=2; % 1=Legendre  2=Cheb
if method==1 % Legendre
    j=[1:n-2]; d=.5*sqrt(j.*(j+2)./(j+.5)./(j+1.5));
    node=eig(diag(d,1)+diag(d,-1)); % Legendre Lobatto pts
    node=[1; flipud(node); -1];
    legn=legendre(n,node); legn=legn(1,:); % L_n(nodes)
    for i=1:n1 % calculate deriv matrix
        legni=legn(i); nodei=node(i);
        for j=1:n1 if i ~= j D(i,j)=legni/legn(j)/(nodei-node(j)); end; end;
    end
    D(1,1)=n*n1/4; D(n1,n1)=-n*n1/4;
else % Cheb
    j=0:n; d=ones(n1,1); d(1)=2; d(n1)=2;
    node=cos(pi*j/n)'; D(1,1)=(2*n*n+1)/6; D(n1,n1)=-D(1,1);
    for i=1:n1 % calculate deriv matrix
        nodei=node(i); di=d(i);  
        if i>1 & i<n1 D(i,i)=-nodei/2/(1-nodei*nodei); end;
        for j=1:n1 if i ~= j D(i,j)=(-1)^(i+j)*di/d(j)/(nodei-node(j)); end; end;
    end
end

A=-D*D+diag(node); 
A=A(2:n,2:n);
f=(pi^2/4+node).*cos(pi*node/2); f=f(2:n); % analytic f
x=A\f; error=norm(x-cos(pi*node(2:n)/2),inf),
xex=1./(1+node.^2); f=xex.*(2*(1-3*node.^2).*xex.^2+node)-.5*node; f=f(2:n); xex=xex-.5; % C^infty but not analytic f
x=A\f; error=norm(x-xex(2:n),inf),
s=1; xex=sqrt(s-node(2:n).^2); f=3./xex.*(s-2*node(2:n).^2)+node(2:n).*xex.^3; % bad f
x=A\f; error=norm(x-xex.^3,inf),
plot(flipud(node),[0;flipud(x);0],'-'),xlabel('x'); ylabel('u');