function testBeam()
clear EA
nu = 0.3; L = 3;

function res = EA(x)
  res = (2-sin(x/L*pi))/2;
end%function

function y = G(z,T)
  nu = 0.3;
  y = nu^2*z.^3-2*nu*z.^2 + z -T;
end%function 

function y = dG(z)
  y = 3*nu^2*z.^2-4*nu*z + 1;
end%function 

N = 500;

%%%%%%%% no modifications necessary beyond this line
h = L/(N+1); % stepsize 

x = (0:h:L)';

clf
relErrorTol=1e-10;  % choose you relative error tolerance

w = zeros(size(x));
wNew = w;

FList = 0.01:0.01:0.24;
% FList = 0.01:0.02:0.3;
maxAmp = zeros(size(FList));
k = 0;
for F = FList
  k = k+1;
  T = F./EA(x);
  relError = 2*relErrorTol;
  while relError>relErrorTol;
    wNew = w-G(w,T)./dG(w);
    relError = max(abs(w-wNew))/max(abs(wNew));
    w = wNew;
  end%while
  maxAmp(k) = trapz(x,wNew);
end%for

u = cumtrapz(x,wNew);
figure(1);
plot(x, u);
grid on
xlabel('position x'); ylabel('displacement u');

figure(2);
plot(maxAmp,FList);
grid on
xlabel('u(L)'); ylabel('F');

end%function
