EA = 1.0;  L = 3; uM = 0.2;   N = 20;
fRHS = @(x) sin(0.5*x); % define the function for the RHS

h = L/(N+1);      % stepsize 
x = (h:h:L-h)';   f = fRHS(x);
g = h^2/EA*f;  g(N) = g(N)+uM;

%% build the tridiagonal, symmetric matrix
di = 2*ones(N,1);      % diagonal
up = -ones(N-1,1);     % upper and lower diagonal

u = trisolve(di,up,g);  % use the special solver
figure(1);
plot([0;x;L],[0;u;uM])  % plot the displacement
xlabel('distance'); ylabel('displacement')
grid on

du = diff([0;u;uM])/h;
figure(2);
plot([0;x]+h/2,EA*du) % plot the force at the cross sections
xlabel('distance'); ylabel('force')
grid on
