function Keller()  % patch to make it work in Matlab
x = 0:0.1:3;

function [y,dy] = f(x)
   y  = x.*x-1-cos(x);
   dy = 2*x+sin(x);
end%function

y = f(x);
figure(1);
plot(x,y)
grid on
xlabel('x'); ylabel('f(x)= x^2-1-cos(x)');

if 0
  % use options to get enough correct digits
  clear options
  options.Jacobian = 'on';
  options.TolX   = 1e-12;
  options.TolFun = 1e-12;
  %[z,fvec,info] = fsolve('f',1.0,options)  % Octave
  z = fzero(@f,1.0)             % Matlab
end%if

N = 2001;
%N = 200001;
h = 2/(N+1); % stepsize 

x = (-1:h:1)';

c = 1.176501940;
%c = z;
uexact = log(c^2./(1+cos(c*x)));


%%%% use Newton's method %%%%%%%%
%% build the triangular matrix
di = 2*ones(N,1)/h^2 ;  % main diagonal
up = -ones(N-1,1)/h^2;  % upper diagonal

Niterations = 5;
errorNewton = zeros(Niterations,1);
errorSubst  = errorNewton;

u = zeros(N,1);
for k = 1:Niterations
  g =  diff(diff([0;u;0]))/h^2-exp(u);
  u = u + trisolve(di+exp(u),up,g);
%  u += trisolve(di+exp(u),up,g);
  errorNewton(k) = max(abs(uexact-[0;u;0]));
end
errorNewton

figure(2);
plot(x,[0;u;0],x,uexact)
xlabel('x'); ylabel('approximate and exact solution');
grid on
%%%%%%% use the method of successive substitutions

u = zeros(N,1);
for k = 1:Niterations
  g = -exp(u);
  u = trisolve(di,up,g);
  errorSubst(k) = max(abs(uexact-[0;u;0]));
end

errorSubst
figure(3);
plot(1:Niterations,log10(errorNewton),1:Niterations,log10(errorSubst))
xlabel('iterations'); ylabel('log of error'); grid on
legend('Newton','Substitution')
end%function
