%  Finite Element Code for BVP 

  clear all
  syms xsym hsym

  N = 10;
  h = 3/N;
  hsym = h;
  x = 0:h:3;
  u = zeros(N+1,1);
  Nq = 4;

  v0 = h*(2/3)*(pi^2)*ones(N-1,1) + (2/h)*ones(N-1,1);
  v1 = h*(1/6)*(pi^2)*ones(N-2,1) + (-1/h)*ones(N-2,1);

  A = diag(v0,0) + diag(v1,1) + diag(v1,-1); 
  f = bvp_rhs(N,Nq,h);

  us = A\f;
  u(2:N,:) = us;

  ftest = double(int((pi^2/hsym)*xsym*sin(pi*xsym/3),0,hsym) + int((pi^2/hsym)*(2*hsym-xsym)*sin(pi*xsym/3),hsym,2*hsym));
  abs(ftest - f(1))

  figure(1)
  plot(x,u,'linewidth',2) 
  %hold on
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')

  P = (-.9/2.25)*(x - 1.5).^2 + 0.9; 
  figure(2)
  plot(x,u,x,P,'--r','linewidth',2)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  legend('Approximate Solution','Polynomial Fit')