%
% equazione della catenaria risolta con il metodo di Newton
% 
l = 3; % >=2
a0 = sqrt(3*(l-2));
a = a0;
delta = -(2*sinh(a)-a*l)/(2*cosh(a)-l);
tol = 1e-8;
while (abs(delta) > tol)
  a = a+delta;
  delta = -(2*sinh(a)-a*l)/(2*cosh(a)-l);
end
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*spdiags(B*u./sqrt(1+(B*u).^2),0,m,m)*B;
  tol = h^2/100;
  u0 = ones(m,1);
  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);
figure(1)
plot(x,u,'*',x,uexact(x),'k')
xlabel('x')
ylabel('u(x)')
h1 = gca();
%f2 = figure(2);
figure(2)
loglog(mrange,errorinf,'*',mrange,errorinf(end)*(mrange/mrange(end)).^(-2),'k')
axis([15,70,1e-5,1e-3])
legend('errore','h^2')
xlabel('m')
ylabel('errore in norma infinito')
h2 = gca();
% check lunghezza corda
l = sum(sqrt(diff(u).^2+diff(x).^2))
