function x = LinearFEM(n)
%       n = number of finite elements

% A function solving the problem 
%           -(k(x) u')' + q(x)u = f(x)   in (0,l)
%           u(0) = lbc, u(l) = rbc
% In particular, this function can be used to solve your
%    first programming assigment. (This is actually it's
%    only purpose - see the notes below.)
  

%NOTE: I have written a "C-style" code ("for"-cycles and such,
%      instead of MatLab's built in vector operations), so 
%      that those of you who use C,C++ or Java can (almost)
%      copy it.

%NOTE: This is by no means a "good" FEM code. As discussed in
%      class, a good FEM code should *NOT* assemble the full 
%      matrix (as it is VERY sparse (tridiagonal in 1D)), and
%      should *NOT* use precalculated element matrices (but use
%      numerical integration instead, so it is easily portable 
%      to elements of arbitrary order).

% Coefficient definitions (the functions k,q and f are defined
%  at the end of the file.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  l=50; 
  lbc=0; rbc=0;

  % Evaluate the element size.
  h=l/n; 
  
  % In all of the following three cycles we go ofer each
  %   finite element and add its contribution to the respective
  %   object (stiffness matrix, mass matrix or right hand side).
  % Naturally those 3 cycles can be combined in one.
  
  % Since we're using linear finite elements, on the j-th
  %   element, only the j-th and (j+1)-th basis funcitons are
  %    non-zero.
  
  % Populate the stiffness matrix.
  % The quantities +/- 1/h are the integrals of phi_i'*phi_j'
  %     over an element (phi_i - basis functions)
  A1=zeros(n+1,n+1);
  for i=1:n
      x = i*h-h/2;
      A1(i,i) = A1(i,i) + 1/h*k(x);
      A1(i,i+1) = A1(i,i+1) - 1/h*k(x);
      A1(i+1,i) = A1(i+1,i) - 1/h*k(x);
      A1(i+1,i+1) = A1(i+1,i+1) + 1/h*k(x);
  end      

  % Populate the mass matrix.
  % The quantities h/3,h/6 are the integrals of phi_i*phi_j
  %     over an element (phi_i - basis functions)
  A0=zeros(n+1,n+1);
  for i=1:n
      x = i*h-h/2;
      A0(i,i) = A0(i,i) + h/3*q(x);
      A0(i,i+1) = A0(i,i+1) + h/6*q(x);
      A0(i+1,i) = A0(i+1,i) + h/6*q(x);
      A0(i+1,i+1) = A0(i+1,i+1) + h/3*q(x);
  end

  % ... and the right hand side.
  % Again, h/2 is an integral of a basis function over
  %     an element.
  B = zeros(n+1,1);
  for i=1:n
      x = i*h-h/2;
      B(i) = B(i) + h/2*f(x);
      B(i+1) = B(i+1) + h/2*f(x);
  end

  % The global matrix is now:
  A=A0+A1;

  % Imposing the boundary conditions.
  %   Left,
  B(1)=lbc; B(2)=B(2)-A(2,1)*B(1);   
  A(2,1)=0; A(1,1)=1;  A(1,2)=0;   
  %   Right.
  B(n+1)=rbc;  B(n)=B(n)-A(n,n+1)*B(n+1);
  A(n,n+1)=0; A(n+1,n)=0; A(n+1,n+1)=1;
 
  % Solving, 
  x=A\B

  % and plotting the solution against a uniform mesh.
  plot(0:h:l,x)
  
  discrete_max = max(abs(x'-w(0:h:l)));
  
  % 02/06 update - code that computes the H1-norm of the error.
  %  We do element by element integration using the Simpson's rule.
  l2_norm = 0;
  for i=1:n
      mid = i*h-h/2;
      l2_norm = l2_norm + h/6*( (x(i)-w(mid-h/2))^2 + 4*((x(i)+x(i+1))/2 - w(mid))^2 + (x(i+1) - w(mid+h/2))^2 );
  end
  l2_norm = sqrt(l2_norm);
  
  h1_seminorm = 0;
  for i=1:n
      mid = i*h-h/2;
      der = (x(i+1)-x(i))/h;
      h1_seminorm = h1_seminorm + h/6*( (der-w_prime(mid-h/2))^2 + 4*(der - w_prime(mid))^2 + (der - w_prime(mid+h/2))^2 );
  end
  h1_seminorm = sqrt(h1_seminorm);
  
  h1_norm = sqrt(l2_norm^2 + h1_seminorm^2)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % The variable coefficents of the problem we're solving.
  function y = k(x)
      y = 1;
  end
  
  function y = q(x)
      y = 100/(8.8*10^7);
  end
  
  function y = f(x)
    y=200/(2*8.8*10^7)*x*(l-x);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % The solution (added for computing the errors)
  function y = w(x)
      a = 100/(8.8*10^7)*l^2;
      b = 200/(2*8.8*10^7)*l^4;
      t = x/l;
      
      y = b/a * (-t.^2+t-2/a+2/(a*sinh(sqrt(a)))*(sinh(sqrt(a)*t)+sinh(sqrt(a)*(1-t))));
  end

  % (Approximate) Derivative of the solution (added for computing the errors)
  function y = w_prime(x)
    delta = 0.001;
    y = (w(x+delta)-w(x-delta))/(2*delta);
  end

end