%MatLab fudged rework Dasgupta Columbia Fin Elm E6332 17SEP81 MLab homework
d1=2;d2=5;l=30;fend=10;e=3e4;n=10;uniloc=[1 -1; -1 1];kgl=zeros(n+1,n+1);x=0;
for el=1:n
  area=(pi/4)*(d2)^2+((pi/4)*(2*d2-d1)^2-(pi/4)*(d2)^2)*(el-1)/n;
  xold=x;     x=((((pi/4)*(d2+(d2-d1)*el/l)^2/(pi/4))^.5)-d2)*l/(d2-d1);
              if el==1, delx=x; else delx=(x-xold);end
  kloc=e*area/delx*uniloc;
  for i=1:2,for j=1:2,kgl(el+i-1,el+j-1)=kgl(el+i-1,el+j-1)+kloc(i,j);end;end
end
force=fend*(1:n)/n;force(1)=0;force(n)=fend;disp(1)=0;
force(n-1)=force(n-1)-kgl(n-1,n)*fend;kgl(n-1,n)=0;kgl(n,n-1)=0;kgl(n,n)=1;
disp(2:n)=inv(kgl(2:n,2:n))*force(2:n)';sum(disp)