clear
EI = 1.0; L = 3; F = [0 0.1];
N = 100;
%F = [0,2];
%%%%%%%%  no modifications necessary beyond this line
h = L/N;  s = (h:h:L)';

%% build the tridiagonal matrix
di = 2*ones(N,1)/h^2;  % diagonal
di(N) = 1/h^2;
up = -ones(N-1,1)/h^2;     % upper diagonal
lo = -ones(N-1,1)/h^2;     % lower diagonal

figure(1)
clf; hold on
alpha = zeros(size(s));
for k = 1:10
  g = F(2)/EI*cos(alpha); g(N) = g(N)/2;
  alpha = trisolve(lo,di,up,g);
  plot(s,alpha)
end%for
hold off
x = cumtrapz([0;s],[1; cos(alpha)]); y = cumtrapz([0;s],[0; sin(alpha)]);
figure(2)
plot(x,y)
xlabel('x'); ylabel('y'); grid on

MaximalAngle = [alpha(N) F(2)/(2*EI)*L^2]
MaximalDeflections = [max(y) F(2)/(3*EI)*L^3]
