clear all
close all
a = 1;
m = 100;
x = linspace (-1, 1, m)';
h = 2 / (m - 1);
D2 = toeplitz (sparse ([1, 1], [1, 2], ...
    [-2, 1] / h ^ 2, 1, m));
D2(1, 1:2) = [1 / h ^ 2, 0];
D2(m, m - 1:m) = [0, 1 / h ^ 2];
D1 = toeplitz (sparse (1, 2, ...
    -1 / (2 * h), 1, m), sparse (1, 2, ...
    1 / (2 * h), 1, m));
D1(1, 2) = 0;
D1(m, m - 1) = 0;
b = [1 / h ^ 2 - a; zeros(m - 2, 1); ...
    1 / h ^ 2 - a];
F = @(u) D2 * u - ...
    a * sqrt (1 + (D1 * u) .^ 2) - b;
JF = @(u) D2 - a * diag ((D1 * u) ./ ...
    sqrt (1 + (D1 * u) .^ 2)) * D1;
u0 = ones (m , 1);
tol = h ^ 2 / 100;
u = u0;
delta = -JF(u) \ F(u);
while (norm (delta, inf) > tol)
    u = u + delta;
    delta = -JF(u) \ F(u);
end
u = u + delta;
plot (x, u, '*', ...
    x, cosh (a * x) / a - cosh (a) / a + 1)
norm (u - ...
    (cosh (a * x) / a - cosh (a) / a + 1), inf)
