function mat = matdif (xx,eps)
global nb nv barras restric massa x0

mat=zeros(2*nv);
for j=1:2*nv
  ej=zeros(2*nv,1);ej(j)=1;
  dforce = (res(xx+eps*ej)-res(xx-eps*ej))/(2*eps);
  mat(:,j)=dforce;
end