% finite element method demo #1
% -u''(x) = 1, u(0)=u(1)=0
% exact solution: u(x) = -(x^2)/2 + x/2;

 N = input('input N: ');
 h = 1/(N+1);

 x = chebfun('x',[0 1]);
 f = x.^0;                   % right-hand size, f(x) = 1
 u = -(x.^2)/2 + x/2;        % exact solution for f(x) = 1
                             % build stiffness matrix
 A = (1/h)*(2*eye(N)-diag(ones(N-1,1),-1)-diag(ones(N-1,1),1));

 ff = zeros(N,1);            % build load vector
 clear Phi
 for j=1:N 
    Phi{j} = hat(j,N);
    ff(j) = f'*Phi{j};       % f_j = (f,phi_j)
 end  

 c = A\ff;

 u_N = 0*x;                  % initialize u_N = 0
 for j=1:N
    u_N = u_N + c(j)*Phi{j}; % build u_N term-by-term
 end

 figure(1),  clf
 plot(u, 'b-','linewidth',2)
 hold on
 plot(u_N, 'r-','linewidth',2)
 leg = legend('$u(x)$', '$u_N(x)$','orientation','horizontal','location','northoutside');
 set(leg,'fontsize',16,'interpreter','latex')

 figure(2),  clf
 plot(u-u_N, 'r-','linewidth',2)
 leg = legend('$u(x)-u_N(x)$','orientation','horizontal','location','northoutside');
 set(leg,'fontsize',16,'interpreter','latex')
