clear
EI = 1.0; L = 3; F = [0 0.1];
F = [0 2] % large force
N = 300;

%%%%%%%%  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

g = F(2)/EI*ones(size(s)); g(N) = 0;

alpha = trisolve(lo,di,up,g);

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

[max(y) F(2)/(3*EI)*L^3]
