##setenv PYTHON python3
syms f1 f2 x y1 y2 s k1 k2 c m1 m2
## m1*s^2*y1 == s*c*y2 -s*c*y1 + k1*y2 -k1*y1
f1 = m1*s^2 * y1 == s*(c*(y2-y1)) + k1*(y2-y1)
f2 = m2*s^2 * y2 == s*(c*(y1-y2)) + k1*(y1-y2) + k2*(x-y2)
[y1a, y2a]=solve(f1,f2,y1,y2)
[n1 d1]=numden(y1a) # this is the bounce of the car body
mm1=sym('500')
mm2=sym('50')
kk2=sym('300000')
kk1=sym('10000')
cc=sym('1000')
d2=vpa(subs(d1,{ m1,m2,k2,k1,c},{ mm1, mm2, kk2,kk1,cc}),20)
d3=simplify(d2)
n2=vpa(subs(n1,{k1,k2,c},{kk1,kk2,cc}),20)
n3=simplify(n2)
f3=n3/d3
f4=factor(f3)
f5=f4/x
[n, d] = numden(f5)
num = double(sym2poly(n,s))
den = double(sym2poly(d,s))
sys=tf(num,den);
step(sys,3)
[n1 d1]=numden(y2a) # this is the bounce of the car wheel
d2=vpa(subs(d1,{ m1,m2,k2,k1,c},{ mm1, mm2, kk2,kk1,cc}),20)
d3=simplify(d2)
n2=vpa(subs(n1,{k1,k2,c,m1},{kk1,kk2,cc,mm1}),20)
n3=simplify(n2)
f3=n3/d3
f4=factor(f3)
f5=f4/x #Remove the x because the step function applies a step to the transfer function
[n, d] = numden(f5)
num = double(sym2poly(n,s))
den = double(sym2poly(d,s))
sys=tf(num,den)
hold on
step(sys,3)
hold off