a = 1;
mrange = [20:10:60];
counter = 0;
uexact = @(x) cosh(a*x)/a-cosh(a)/a+1;
for m = mrange+1
  counter = counter+1;
  x = linspace(-1,1,m)';
  h = 2/(m-1);
  A = toeplitz(sparse([1,1],[1,2],[-2/h^2,1/h^2],1,m));
  A(1,1:2) = [1/h^2,0];
  A(m,m-1:m) = [0,1/h^2];
  B = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
  B(1,2) = 0;
  B(m,m-1) = 0;
  F = @(u) A*u-a*sqrt(1+(B*u).^2)-[1/h^2-a;zeros(m-2,1);1/h^2-a];
  J = @(u) A-a*diag(a*B*u./sqrt(1+(B*u).^2))*B;
  tol = h^2/100;
  u0 = x.^2;
  u = u0;
  errest = -J(u)\F(u);
  while (norm(errest,inf) > tol)
    u = u+errest;
    errest = -J(u)\F(u);
  end
  errorinf(counter) = norm(u-uexact(x),Inf);
end
f1 = figure(1);
plot(x,u,'*',x,uexact(x),'b')
xlabel('x')
ylabel('u(x)')
h1 = gca();
f2 = figure(2);
loglog(mrange,errorinf,'*',mrange,errorinf(end)*(mrange/mrange(end)).^(-2),'b')
axis([15,70,1e-5,1e-3])
legend('error','order 2')
xlabel('m')
ylabel('error in infinity norm')
h2 = gca();