function [x,sol]=finel(i,a,b,n);
%
% [x,sol]=finel(i,a,b,n)
% risolve -u''=f, u(a)=esatta(a), u(b)=esatta(b)
% con gli elementi finiti
% n=numero di intervalli di suddivisione
% i seleziona il dato f
% prerequisito: avere "f" e "esatta"
%
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)';
%
% uso una formula di Gauss a 8 punti
%
%xx=[-.9061798459 -.5384693101 .5384693101 .9061798459];
xx=[-.9602898565 -.7966664774 -.5255324099 -.1834346425...
     .1834346425  .5255324099  .7966664774  .9602898565];
bb1=xx/2+1/2;
bb2=-xx/2+1/2;
%pp=[.2369268851 .4786286705 .4786286705 .2369268851];
pp=[.1012285363 .2223810345 .3137066459 .3626837834...
    .3626837834 .3137066459 .2223810345 .1012285363];
pp=pp/sum(pp);
g=zeros(n-1,1);
for ii=1:n-1
  xx1=bb1*h+x(ii);
  xx2=xx1+h;
  g(ii)=h*(sum(pp.*f(i,xx1).*bb1)+sum(pp.*f(i,xx2).*bb2));
end
%
sol=zeros(n+1,1);
sol(1)=esatta(i,x(1)); sol(n+1)=esatta(i,x(n+1));
g(1)=g(1)+sol(1)/h; g(n-1)=g(n-1)+sol(n+1)/h;
sol(2:n)=A\g;
