function interp_erf
	x = [0,0.5,1]
	y = erf(x)
%	plot(x,y,'b*.');
%   	hold on;

	n = numel(x)
	c = newtdd(x,y,n);
	
	xx = linspace(x(1),x(n),200);
	yy = nest(n-1,c,xx,x);
	err = yy-erf(xx);
        errbnd = (4/6/sqrt(pi)).*(xx).*(xx-0.5).*(xx-1);
%	plot(xx,yy,'r');
%       plot(xx,erf(xx),'k');
	plot(xx,err,'r');	hold on;
	plot(xx, errbnd,'g');
	plot(xx,-errbnd,'g');   hold off;
	title('Error bound and actual error (red) of  quadratic interpolant for erf on [0,1]');
end