function [G]=example2
% Hydraulic brake system for a front-end loader.
warning off
Vis=10e-6;r=.0032;Be=1.7235e8;Den=820;c=sqrt(Be/Den);
Zo=Den*c/(pi*r^2);V=2.46e-6;A=0.0111;K=3.78e6;
La=.559;Dna=Vis*La/(c*r^2);Ld=2.134;Dnd=Vis*Ld/(c*r^2);
Lb=.762;Dnb=Vis*Lb/(c*r^2);Lc=.7112;Dnc=Vis*Lc/(c*r^2);
syms s
B=2*besselj(1,j*sqrt(r^2*s/Vis))/(j*sqrt(r^2*s/Vis)*besselj(0,j*...
sqrt(r^2*s/Vis))); Z=Zo/sqrt(1-B);
Ga=Dna* r^2*s/Vis/sqrt(1-B);Gd=Dnd* r^2*s/Vis/sqrt(1-B);
Gc=Dnc* r^2*s/Vis/sqrt(1-B);Gb=Dnb* r^2*s/Vis/sqrt(1-B);
Sol=solve('Q2=C2a*1-C1a*P2','Q3=C1b*P2-C2b*P3',...
'Q4=C2b*P2-C1b*P3','Q5=C1c*P3-C2c*P4',...
'Q6=C2c*P3-C1c*P4','Q7=C1d*P2-C2d*P5',...
'Q8=C2d*P2-C1d*P5','Q9=C1c*P5-C2c*P6',...
'Q10=C2c*P5-C1c*P6','Q2=Q3+Q4','Q4=2*Q5',...
'Q8=2*Q9','Q6=Cb*s*P4','Q10=Cb*s*P6',...
'P2,P3,P4,P5,P6,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10');
H=Sol.P6;collect(H,s);
H=subs(H,'C1a',cosh(Ga)/(Z*sinh(Ga)));
H=subs(H,'C2a',1/(Z*sinh(Ga)));
H=subs(H,'C1d',cosh(Gd)/(Z*sinh(Gd)));
H=subs(H,'C2d',1/(Z*sinh(Gd)));
H=subs(H,'C1c',cosh(Gc)/(Z*sinh(Gc)));
H=subs(H,'C2c',1/(Z*sinh(Gc)));
H=subs(H,'C1b',cosh(Gb)/(Z*sinh(Gb)));
H=subs(H,'C2b',1/(Z*sinh(Gb)));
H=subs(H,'Cb',(V/Be+A^2/K));
w=[(1:.1:9.9),(10:.2:19.8),(20:1:99),(100:5:195),(200:20:300)];
N=length(w);
wc=[(1:.1:9.9),(10:.2:19.8),(20:1:99),(100:5:195),(200:20:1000)];
NC=length(wc);
for k=1:N
sw=j*wc(k); %Generate values for s.
TF(k)=subs(H,s,sw);
TFc(k)=TF(k);
end
N1=N+1;
%Additional frequencies generated for accuracy comparisons.
for k=N1:NC 
sw=j*wc(k);
 TFc(k)=subs(H,s,sw);
end
%Create a transfer function by curve fitting.
[Num,Den]=invfreqs(TF,w,9,10,[],100);
G=tf(Num,Den);
step(G)
figure
Hc=freqs(Num,Den,wc);
MH=20*log10(abs(Hc));
AH=angle(Hc)*180/pi;
MTF=20*log10(abs(TFc));
ATF=angle(TFc)*180/pi;
semilogx(wc,MH,'r',wc,MTF,'b');
title('Magnitude Comparison Plots');
xlabel('Frequency rad/sec');ylabel('Decibels');
figure
semilogx(wc,AH,'r',wc,ATF,'b')
title('Phase Comparison Plots');
xlabel('Frequency rad/sec');ylabel('Degrees');