function Kmat=StiffnessMat(NodeNum,ElemNum,ElemProp,NodeCoord,ElemConec)

Kmat=zeros(2*NodeNum,2*NodeNum);

for i=1:ElemNum
    DeltaX=NodeCoord(ElemConec(i,2),1)-NodeCoord(ElemConec(i,1),1);
    DeltaY=NodeCoord(ElemConec(i,2),2)-NodeCoord(ElemConec(i,1),2);
    Leng=sqrt(DeltaX^2+DeltaY^2);
    C=cos(DeltaY/Leng);
    S=sin(DeltaX/Leng);
    kstiff=ElemProp(i,1)*ElemProp(i,2)/Leng;
    ii=ElemConec(i,1);
    jj=ElemConec(i,2);
    Kmat(2*ii-1,2*ii-1)=Kmat(2*ii-1,2*ii-1)+kstiff*C*C;
    Kmat(2*ii-1,2*ii)=Kmat(2*ii-1,2*ii)+kstiff*C*S;
    Kmat(2*ii-1,2*jj-1)=Kmat(2*ii-1,2*jj-1)-kstiff*C*C;
    Kmat(2*ii-1,2*jj)=Kmat(2*ii-1,2*jj)-kstiff*C*S;
    Kmat(2*ii,2*ii-1)=Kmat(2*ii,2*ii-1)+kstiff*C*S;
    Kmat(2*ii,2*ii)=Kmat(2*ii,2*ii)+kstiff*S*S;
    Kmat(2*ii,2*jj-1)=Kmat(2*ii,2*jj-1)-kstiff*C*S;
    Kmat(2*ii,2*jj)=Kmat(2*ii,2*jj)-kstiff*S*S;
    Kmat(2*jj-1,2*ii-1)=Kmat(2*jj-1,2*ii-1)-kstiff*C*C;
    Kmat(2*jj-1,2*ii)=Kmat(2*jj-1,2*ii)-kstiff*C*S;
    Kmat(2*jj-1,2*jj-1)=Kmat(2*jj-1,2*jj-1)+kstiff*C*C;
    Kmat(2*jj-1,2*jj)=Kmat(2*jj-1,2*jj)+kstiff*C*S;
    Kmat(2*jj,2*ii-1)=Kmat(2*jj,2*ii-1)-kstiff*C*S;
    Kmat(2*jj,2*ii)=Kmat(2*jj,2*ii)-kstiff*S*S;
    Kmat(2*jj,2*jj-1)=Kmat(2*jj,2*jj-1)+kstiff*C*S;
    Kmat(2*jj,2*jj)=Kmat(2*jj,2*jj)+kstiff*S*S;
end

end