% beam.m calculates the position of a loaded beam

n = 80;                 % use n+1 data points
dx = 1/n;               % width of spatial intervals
x = 0:dx:1;             % x coordinates

D2 = zeros(n-1,n+1);    % define the second derivative matrix
for i=1:n-1,
  D2(i,i)   = 1;
  D2(i,i+1) = -2;
  D2(i,i+2) = 1;
end
D2 = D2/(dx^2);         % scale the derivative correctly

S = [ [1 zeros(1,n)]; D2; [zeros(1,n) 1]];  % full matrix for the system

% EXAMPLE 1 ----------------------------------------------------------------

f = zeros(n-1,1);       % the force
b = [2; f; 3];
y = inv(S) * b;

subplot(4,2,1)
plot(x,[0; f; 0]);
title('Force on the beam (zero in this case)')
axis([0 1 0 50])

subplot(4,2,2)
plot(x,y);
title('Unloaded beam')
axis([0 1 0 4])

% EXAMPLE 2 ----------------------------------------------------------------

f = 30*rand(n-1, 1);
b = [2; f; 3];
y = inv(S) * b;

subplot(4,2,3)
plot(x,[0; f; 0]);
title('Force on the beam')
axis([0 1 0 50])

subplot(4,2,4)
plot(x,y);
title('Loaded beam')
axis([0 1 0 4])

% EXAMPLE 3 ----------------------------------------------------------------

f = 80*binornd(2,5/(2*n),n-1, 1);
b = [2; f; 3];
y = inv(S) * b;

subplot(4,2,5)
plot(x,[0; f; 0]);
title('Force on the beam')
axis([0 1 0 160])

subplot(4,2,6)
plot(x,y);
title('Loaded beam')
axis([0 1 0 4])

% EXAMPLE 4 ----------------------------------------------------------------

f = 80*binornd(2,5/(2*n),n-1, 1);
b = [2; f; 3];
y = inv(S) * b;

subplot(4,2,7)
plot(x,[0; f; 0]);
title('Force on the beam')
axis([0 1 0 160])

subplot(4,2,8)
plot(x,y);
title('Loaded beam')
axis([0 1 0 4])