% FEM solver for d2T/dx2 + q = 0 where q = 50 exp(x)
%
% BC's: T(-1) = 100 and T(1) = 100.
%
% Note: the finite element degrees of freedom are
%       stored in the vector, v.  

% Number of elements
Ne = 10;
x = linspace(-1,1,Ne+1);

% Zero stiffness matrix
K = zeros(Ne+1, Ne+1);
b = zeros(Ne+1, 1);

% Loop over all elements and calculate stiffness and residuals
for ii = 1:Ne,

  kn1 = ii;
  kn2 = ii+1;
  
  x1 = x(kn1);
  x2 = x(kn2);
  
  dx = x2 - x1;
  
  % Add contribution to kn1 weighted residual due to kn1 function
  K(kn1, kn1) = K(kn1, kn1) - (1/dx);

  % Add contribution to kn1 weighted residual due to kn2 function
  K(kn1, kn2) = K(kn1, kn2) + (1/dx);
  
  % Add contribution to kn2 weighted residual due to kn1 function
  K(kn2, kn1) = K(kn2, kn1) + (1/dx);

  % Add contribution to kn2 weighted residual due to kn2 function
  K(kn2, kn2) = K(kn2, kn2) - (1/dx);

  % Add forcing term to kn1 weighted residual 
  b(kn1) = b(kn1) + COMPLETE THIS LINE;
  
  % Add forcing term to kn2 weighted residual
  b(kn2) = b(kn2) + COMPLETE THIS LINE;
  
end


% Set Dirichlet conditions at x=0
kn1 = 1;
K(kn1,:)    = zeros(size(1,Ne+1));
K(kn1, kn1) = 1.0;
b(kn1)      = 100.0;


% Set Dirichlet conditions at x=1
kn1 = Ne+1;
K(kn1,:)    = zeros(size(1,Ne+1));
K(kn1, kn1) = 1.0;
b(kn1)      = 100.0;

  
% Solve for solution
v = K\b;


% Plot solution
plot(x,v,'*-');
xlabel('x');
ylabel('T');

% Plot exact solution
xe = linspace(-1,1,2001);
Te = -50*exp(xe) + 50*xe*sinh(1) + 100 + 50*cosh(1);
hold on; plot(xe,Te); hold off;