NR=10;

Coordinates & Differentiation - precomputed

calcCheb[NR_]:=(

NN=2*NR-1;

ChebR=Table[Cos[((π ii)/NN)],{ii,NN,0,-1}];

x[ii_]:=ChebR[[ii+1]];

Dmat=Table[If[ii!=jj,((If[ii==0||ii==NN,2,1])/(If[jj==0||jj==NN,2,1]))(((-1)^(ii+jj))/((x[ii]-x[jj]))),0],{ii,0,NN},{jj,0,NN}];

Do[
Dmat[[ii+1,ii+1]]=((-x[ii])/(2(1-x[ii]^2)));
,{ii,1,NN-1}];

Dmat[[1,1]]=-((2(NN^2)+1)/6);

Dmat[[NN+1,NN+1]]=((2(NN^2)+1)/6);

D2mat=Dmat.Dmat;

)

Check

f[y_]:=Exp[y^2](1+(y^2-(1/3))^2Log[(y^2-(1/3))^2]);

f[y_]:=Exp[y^2](1+(y^2-(4/3))^2Log[(y^2-(4/3))^2]);

f[y_]:=Exp[y^2](1+(y^2-(4/3))^2);

f[y_]:=Exp[y^2](1+Log[(y^2-(4/3))^2]);

dat={};

Table[

calcCheb[num];

fdat=f[ChebR];

val=Max[N[Dmat.fdat-f'[ChebR],30]];

tmp={num,Log[10,val]};

Print[tmp];

AppendTo[dat,tmp];

,{num,4,40,2}];

dat

ListPlot[dat]