function [x,sol]=diffreaz(i,e,a,b,n,f)
% [x,sol]=diffreaz(i,e,a,b,n,f)
% solves -e u''+u=f, u(a)=0, u(b)=0
% with finite elements
% if i=1, standard scheme
% if i=0, mass lumping
h=(b-a)/n;
A=(diag(2*ones(1,n-1))-diag(ones(1,n-2),-1)-diag(ones(1,n-2),1))/h;
x=(a:h:b)';
g=h*feval(f,(x(1:n-1)+x(2:n))/2);
sol=zeros(n+1,1);
if i==1
  A=e*A+(1/6*diag(ones(1,n-2),1)+2/3*diag(ones(1,n-1))+ ...
         1/6*diag(ones(1,n-2),-1))*h;
elseif i==0
  A=e*A+(diag(ones(1,n-1)))*h;
else
  error('i does not belong to {0,1}')
end 
sol(2:n)=A\g;
