%clear
EI = 1.0; L = 3; F = [0 0.1];  % try values of 0.5 1.5 and 2
F = [0 2];
N = 200;

%%%%%%%%   no modifications necessary beyond this line
h = L/N; % stepsize 
s = (h:h:L)';
alpha = zeros(size(s));

%% build the tridiagonal matrix
di = 2*ones(N,1)/h^2; di(N) = 1/h^2;   % diagonal
up = -ones(N-1,1)/h^2; lo = up;        % lower and upper diagonal
DiffTol = 1e-10;  DiffAbs = 2*DiffTol;
while DiffAbs>DiffTol
  b = F(2)/EI*sin(alpha); b(N) = b(N)/2;
%  g = -[0; lo.*alpha(1:N-1)]-di.*alpha-[up.*alpha(2:N);0];
  g = ([0;alpha(1:N-1)]-[2*alpha(1:N-1);alpha(N)]+[alpha(2:N);0])/h^2;
  g = g + F(2)/EI*cos(alpha);
  g(N) = g(N) - F(2)/EI*cos(alpha(N))/2;
  phi = trisolve(lo,di+b,up,g);
  alpha = alpha+phi;
  DiffAbs = max(abs(phi));    % display maximal difference
  MaxAngle = max(abs(alpha)); % display maximal angle
  disp(sprintf('maximal angle = %7.4f, difference = %7.4e',max(abs(alpha)),max(abs(phi))))
end%while

x = cumtrapz([0;s],[0; cos(alpha)]);  y = cumtrapz([0;s],[0; sin(alpha)]);
plot(x,y)
xlabel('x'); ylabel('y'); grid on
