conec=[1 4; 1 8; 1 2; 2 3; 3 5; 3 6; 2 7];
coord=[1 1; 2 1; 3 1; 0 0;4 0;4 2;2 2;0 2];
kmola=ones(7,1); kmola(7)=1;
m=[1; 2; 3]; 

nb=length(conec);
nv=length(coord);
Kglo=zeros(2*nv,2*nv);
for ib=1:nb
Kloc=zeros(4,4);
na=conec(ib,1); nb=conec(ib,2);
Xa=coord(na,:)';
Xb=coord(nb,:)';
d=Xb-Xa;
l0=norm(d);
kk=kmola(ib)/(l0^2);
aux=d*d';
Kloc(1:2,1:2)=kk*aux;
Kloc(1:2,3:4)=-kk*aux;
Kloc(3:4,1:2)=-kk*aux;
Kloc(3:4,3:4)=kk*aux;
%Kloc agora contem a matriz de rigidez 
% da mola ib
loc2glo=[2*na-1,2*na,2*nb-1,2*nb];
for j=1:4
for k=1:4
jglo=loc2glo(j);
kglo=loc2glo(k);
Kglo(jglo,kglo)=Kglo(jglo,kglo)+Kloc(j,k);
end
end
end

Kmat=Kglo(1:6,1:6);
Mmat=diag([m(1) m(1) m(2) m(2) m(3) m(3)]);