EI = 1.0; L = 3;  N = 200;
FList = 0.25:0.25:2;

%%%%%%%%  no modifications necessary beyond this line
h = L/N; % stepsize 
s = (h:h:L)';
alpha1 = 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

errTol = 1e-10;
figure(1);  clf; hold on
yPlot = zeros(length(s)+1,1); xPlot = [0;s];

for F2 = FList
  F = [0 F2]
  errAbs = 2*errTol;
  while errAbs>errTol;
    b = F(2)/EI*sin(alpha1); b(N) = b(N)/2;
    g = ([0;alpha1(1:N-1)]-[2*alpha1(1:N-1);alpha1(N)]+[alpha1(2:N);0])/h^2;
    g = g + F(2)/EI*cos(alpha1);
    g(N) = g(N) - F(2)/EI*cos(alpha1(N))/2;
    phi = trisolve(lo,di+b,up,g);
    alpha1 = alpha1 + phi;
    errAbs = max(abs(phi))
  end%while
  x = cumtrapz([0;s],[0; cos(alpha1)]);
  y = cumtrapz([0;s],[0; sin(alpha1)]);
  yPlot = [yPlot y]; xPlot = [xPlot x];
  plot(xPlot,yPlot)
  grid on; xlabel('x'); ylabel('y');
  axis equal
end%for

hold off 
figure(2);
x = cumtrapz([0;s],[0; cos(alpha1)]);  y = cumtrapz([0;s],[0; sin(alpha1)]);
plot(x,y); grid on; xlabel('x'); ylabel('y');
