clear all
close all
% u''(x)=k*sqrt(1+u'(x)^2)
% u(-1) = 1;
% u(1) = 1;
k = 1;
mrange = [20:10:60];
counter = 0;
for m = mrange
  counter = counter+1;
  x = linspace(-1,1,m)';
  h = 2/(m-1);
  A = toeplitz(sparse([-2/h^2,1/h^2,zeros(1,m-2)]));
  A(1,:) = zeros(1,m);
  A(1,1) = 1/h^2;
  A(m,:) = zeros(1,m);
  A(m,m) = 1/h^2;
  B = toeplitz(sparse([0,-1/(2*h),zeros(1,m-2)]),...
	       sparse([0,1/(2*h),zeros(1,m-2)]));
  B(1,:) = zeros(1,m);
  B(m,:) = zeros(1,m);
  F = @(u) A*u-k*sqrt(1+(B*u).^2)-[1/h^2-k;zeros(m-2,1);1/h^2-k];
  J = @(u) A-k*diag(k*B*u./sqrt(1+(B*u).^2))*B;
  tol = 1e-6;
  u0 = x.^2;
  u = u0;
  errest = -J(u)\F(u);
  while (norm(errest) > tol)
    u = u+errest;
    errest = -J(u)\F(u);
  end
  errorinf(counter) = norm(u-(cosh(k*x)/k-cosh(k)/k+1),Inf);
end
plot(x,u,'*',x,cosh(k*x)/k-cosh(k)/k+1)
figure
loglog(mrange,errorinf,'*',mrange,errorinf(1)*(mrange/mrange(1)).^(-2))
legend('error','order 2')
xlabel('m')
ylabel('error')
