%求單跨梁振型及頻率（兩端任意約束）
b=input('please input b.' );
h=input('please input h.' );
Rho=input('please input Rho.' );
E=input('please input E.' );
ds=input('How many elments do you want to devide?');
l=input('Pleaes input the length');
Lres=input('Pleaes input the left restraint,Fix->3,Hinged2->2,Hinged1->1,Free->0 ');
Rres=input('Pleaes input the left restraint,Fix->3,Hinged2->2,Hinged1->1,Free->0 ');
I=b*(h^3)/12;
EI=E*I;

x1=0;x2=sym('L');
x=sym('x');
j=0:3;v=x.^j;
a=[1,x1,x1^2,x1^3;
0,1,2*x1,3*x1^2;
1,x2,x2^2,x2^3;
0,1,2*x2,3*x2^2];
d=v/a;
dt=d';
m=dt*d;
M=Rho*b*h*int(m,x,0,'L');
M=subs(M,'L',l/ds);
M=double(M);
Ni=diff(d,x,2);
Nt=Ni';
k=Nt*Ni;
KL=EI*int(k,x,0,'L');
KL=subs(KL,'L',l/ds);
KL=double(KL);

KGC=cell(1,ds);
for n=1:ds
KGC{n}=zeros((ds+1)*2);
end
for n=1:ds;
for i=1:4
for j=1:4
KGC{n}(i+2*n-2,j+2*n-2)=KL(i,j);
end
end
end 
KG=zeros((ds+1)*2);
for n=1:ds
KG=KGC{n}+KG;
end

MGC=cell(1,ds);
for n=1:ds
MGC{n}=zeros((ds+1)*2);
end
for n=1:ds;
for i=1:4
for j=1:4
MGC{n}(i+2*n-2,j+2*n-2)=M(i,j);
end
end
end
MG=zeros((ds+1)*2);
for n=1:ds
MG=MGC{n}+MG;
end

if Lres==3;
MG=MG(3:end,3:end);KG=KG(3:end,3:end);
for n=1:ds;
KGC{n}=KGC{n}(3:end,3:end);
end
elseif Lres==2|Lres==1;
MG=MG(2:end,2:end);KG=KG(2:end,2:end);
for n=1:ds;
KGC{n}=KGC{n}(2:end,2:end);
end
else MG=MG;KG=KG; 
end

if Rres==3;
MG=MG(1:end-2,1:end-2);KG=KG(1:end-2,1:end-2);
for n=1:ds;
KGC{n}=KGC{n}(1:end-2,1:end-2);
end
elseif Rres==2|Rres==1;
MG=MG([1:end-2,end],[1:end-2,end]);KG=KG([1:end-2,end],[1:end-2,end]);
for n=1:ds;
KGC{n}=KGC{n}([1:end-2,end],[1:end-2,end]);
end
else MG=MG;KG=KG;
end