clear all
close all
mrange = 9; % m has to be odd
a = -1;
b = 1;
xx = linspace(a,b,6001);
f = @(x) abs(x).^(1/2);
f1 = @(x)(1/2)*abs(x).^(-1/2).*sign(x);
counter = 0;
for m = mrange;
  counter = counter+1;
  t = linspace(a,0,(m-1)/2-1)';
  t = [t;0];
  t = [t;linspace(0,b,(m-1)/2-1)'];
  t = [t(1)*ones(3,1);t;t(m-2)*ones(3,1)];
  x=[((t(2:end-2)+t(3:end-1)))/2];
  x((m+3)/2)=[];
  h(counter) = max(diff(x));
  y = f(x); 
  y(2)= f1(x(2)); 
  y(end) = f1(x(end));
  s = spapi(t,x,y); 
  yy = spval(s,xx); 
  yyvin = spline(x(2:end-1),[y(2);y(1);y(3:end)],xx);
  errbspline(counter) = norm(yy-f(xx),inf);
  errvin(counter) = norm(yyvin-f(xx),inf);
end
loglog(h,errbspline,'o',h,errvin,'*',h,h.^(0.5))
grid
axis([1e-3,1,1e-12,1e4])
legend('bspline','vincolated','location','northwest')
xlabel('h')
ylabel('error in infinity norm')
figure
plot(xx,f(xx),xx,yy,xx,yyvin)
legend('function','bspline','vincolated')