clear all
close all
p = 3; % smooth case
mrange = 50:100:450;
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(0,1,m)';
  h = x(2)-x(1);
  A = toeplitz(sparse([1,1],[1,2],[2,-1]/h^2,1,m));
  A(1,1) = 1e16;
  A(m,m) = 1e16;
  g = -(2*p+3)*(2*p+1)/4*sqrt(abs(x-1/2)).^(2*p-1);
  g(1) = 0;
  g(m) = 0;
  u = A\g;
  uexact = sqrt(abs(x-1/2)).^(2*p+3)-sqrt(1/2)^(2*p+3);
  err(counter) = norm(u-uexact)/sqrt(m-1);
end
loglog(mrange,mrange.^(-2),mrange,err,'*g')
legend('order 2','error in L2 approx. norm')
axis([10,1e3,1e-7,1e-1])
xlabel('m')
ylabel('error')
print('-depsc','fd1dsmooth.eps')
p = 1; % less smooth case
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(0,1,m)';
  h = x(2)-x(1);
  A = toeplitz(sparse([1,1],[1,2],[2,-1]/h^2,1,m));
  A(1,1) = 1e16;
  A(m,m) = 1e16;
  g = -(2*p+3)*(2*p+1)/4*sqrt(abs(x-1/2)).^(2*p-1);
  g(1) = 0;
  g(m) = 0;
  u = A\g;
  uexact = sqrt(abs(x-1/2)).^(2*p+3)-sqrt(1/2)^(2*p+3);
  err(counter) = norm(u-uexact)/sqrt(m-1);
end
figure
loglog(mrange,mrange.^(-3/2),mrange,err,'*g')
legend('order 3/2','error in L2 approx. norm')
axis([10,1e3,1e-7,1e-1])
xlabel('m')
ylabel('error')
print('-depsc','fd1dnonsmooth.eps')