function BeamNL()
nu = 0.3;
L = 3;
F = 0.2;

function res = EA(x)
  res = (2-sin(x/L*pi))/2;
end%function

function y = fRHS(x)
  y = zeros(size(x));
end%function
%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 100; h = L/(N+1); % stepsize 

x = (h:h:L)';
f = [fRHS(x)];
u = zeros(size(f));
g = h^2*f;  g(N+1) = g(N+1)/2 + h*F;

cc = 1; % color counter
cstring = 'rgbcmykrgbcmyk'; % color sequence

figure(1)
clf
hold on
grid on
axis([0 3 0 1.4])
xlabel('position x'); ylabel('displacement u')
relDiffTol = 1e-5;  % choose you relative difference tolerance
Differences = []; relDiff = 2*relDiffTol;
while relDiff > relDiffTol
  a = EA(x-h/2).*(1-nu*diff([0;u])/h).^2; % compute coefficients
  di = [a(1:N)+a(2:N+1); a(N+1)];         % diagonal entries
  up = -a(2:N+1);                         % upper diagonal entries
  uNew = trisolve(di,up,g);               % solve the linear symmetric system 
  plot([0;x],[0;u],'linewidth',2,'color',cstring(cc)); cc = cc+1;
  pause(0.5)
  relDiff = max(abs(u-uNew))/max(abs(uNew)) % determine relative difference
  Differences = [Differences;relDiff];      % store the relative differences
  u = uNew;                                 % prepare for restart
end%while
axis([0 3 0 1.4])
xlabel('position x'); ylabel('displacement u')

figure(2)
semilogy(Differences)
xlabel('iterations'); ylabel('relative difference')
end%function
