function [Ht]=Example1
% Fluid line terminating into a fluid capacitance
warning off
Vis=2.0e-5;r=.005;Be=1.896e9;Den=875;
c=sqrt(Be/Den);Zo=Den*c/(pi*r^2);
V=.00142;L=5;Ct=V/Be;R=1.776e8;
syms s C1 C2
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);Dn=Vis*L/(c*r^2);G=Dn*(r^2*s/Vis)/sqrt(1-B);
H=C2/((1+C1*R)*Ct*s+C1);%Trans. fun. to be curve fitted.
H=subs(H,'C1',cosh(G)/(Z*sinh(G)));
H=subs(H,'C2',1/(Z*sinh(G)));
% Generate frequency data for the curve fit. 
% Note, the last value is the range of interest.
w=[(10:1:195),(200:10:980),(1000:20:2500)];
N=length(w);
% Generate a greater frequency range for an accuracy comparison
wc=[(10:1:195),(200:10:980),(1000:20:6000)];
NC=length(wc);
for k=1:N
sw=j*wc(k); % Generate values for s.
TF(k)=subs(H,s,sw); % Generate tf data points for the curve fit.
TFc(k)=TF(k);
end
N1=N+1;
for k=N1:NC %Additional frequencies for accuracy comparisons.
 sw=j*wc(k);
TFc(k)=subs(H,s,sw);% Additional data for the accuracy comparison.
end
[Num,Den]=invfreqs(TF,w,9,10,[],100);%Create a tf by curve fitting.
Ht=tf(Num,Den)
step(Ht,0.03)
% Generate Freq. Resp. plots to determine the accuracy of curve fit.
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');