% Solves -d^2u/dx^2 =x   u(0)=u(1)=0 using piecewise linear FEM
% Compares to the true solution
N = input('size of fem subspace  ');
% use uniform grid
h = 1/(N+1);
% build stiffness matrix (the Gram matrix with respect to energy inner product)
A = (1/h)*(2*eye(N)-diag(ones(N-1,1),1)-diag(ones(N-1,1),-1));
% right hand side
f=h^2*(1:N);
% solve for unknown coefficients of solutions
u=A\f';
% coefficients are just the nodal values of the solution, can plot directly
xnodes=0:h:1;
% add on the zero boundary values
u= [ 0 ; u ; 0];
% plot true solution on a fine grid so it looks smooth
finex=0:0.01:1;
plot(xnodes,u,finex,finex./6 -finex.^3./6);
legend('fem solution', 'true solution');