Dr. Boris J. Lurie
 HOME | CONTACT | RESUME | PUBLICATIONS | CONSULTING | LECTURES | SITE MAP

LIST OF PUBLICATIONS

 Books Patents Journal Publications NASA Tech Briefs Conference Proceedings Publications
 Discussions Readers' Contributions Bode Step Toolbox Function Listings

MATLAB AND SPICE SCRIPTS

These are the listings of all MATLAB and SPICE scripts found in Classical Feedback Control by B. Lurie and P. Enright, exclusive of the functions in Appendix 14. The scripts can be copied freely for nonprofit use. The author bears no responsibility for the results of the scripts' applications.

 Chapter 1 Chapter 2 Chapter 3 Chapter 4 Chapter 5 Chapter 6 Chapter 7 Chapter 10 Chapter 11 Appendix 2 Appendix 5 Appendix 13

## Chapter 1

#### 1.4.1  Gain and phase responses

Example 3.  The frequency responses can be plotted with MATLAB ® script:

```w = logspace(-1, 2);
% log scale of angular
% frequency w
num = 5000;
den = [1 55 250 0];
bode(num, den, w)
a = [1 0]; b = [1 5]; c = [1 50];
ab = conv(a,b); den = conv(ab,c)
```

#### 1.4.2  Nyquist diagram

Example 1.  The L-plane Nyquist diagram is charted with MATLAB script

```num = [20 10]; den = [1 10 20 1 0];
[mag, phase] = bode(num, den);
plot(phase, 20*log10(mag), 'r', -180, 0, 'wo')
title('L-plane Nyquist diagram')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
```

### 1.6  Example of system analysis

The frequency response can be plotted using MATLAB with:

```w = logspace(-2,1);
% frequency range
% 0.01 to 10 rad/sec
num = 1;
den = [1 2.1 0.2 0 0];
bode(num, den, w)
```

The plots can be made in MATLAB with:

```w = logspace(-1,1);             % freq range 0.1 to 10 rad/sec
den = [1 15 50 0 0 0];
num = [0 0 0 50 27.5 0.25];     % equal length of the vectors
g = num + den;                  % makes the addition allowable
bode(num, den, w)               % for T
hold on
bode(g, den, w)                 % for F
bode(num, g, w)                 % for M
hold off
```

Example 1.  The transient response is found with:

```num = [0 0 0 50 27.5 0.25];
den = [1 15 50 0 0 0];
g = num + den; step(num, g)
grid
```

Back

## Chapter 2

### 2.10  Problems

22  (a) Exact solution found with MATLAB script:

```A = [2 0.2 0.3; 0.1 2.1 0.1; 0.04 0.1 1.9]; inv(A)
ans =    0.5038   -0.0443   -0.0772
-0.0235    0.4795   -0.0215
-0.0094   -0.0243    0.5291
```

Back

## Chapter 3

### 3.3  Root locus

Example 1.  The open-loop and closed-loop poles can be calculated with the MATLAB commands

```n = [0 0 10 20]; d = [1 1 1 0];
roots(d)       %open-loop poles
ans = 0       -0.5000 ± 0.8660i
roots(n + d)   %closed-loop poles
ans = -1.6551  0.3275 ± 3.4608i
```

The root loci in Fig. 3.11 are plotted with:

```n = [10 20]; d = [1 1 1 0];
rlocus(n,d)
hold on
k = [0.05 0.1 0.2 0.5 1];
rlocus(n, d, k)
title('k = [0.05 0.1 0.2 0.5 1]');
hold off
```

#### 3.9.6  Phase-gain relation

The weight function is charted with the script

```u = linspace(-3,3,200);
b = log(coth(abs(u/2))); plot(u,b,'w');
grid
```

### 3.14  Problems

14  The following is a SPICE input file:

``` ** unst_pl.cir for feedback loop with unstable plant
* compensator, inverting, zero s=10, pole s=40,
* crossover s=20
GC 2 0 1 0 1
RC1 1 0 1MEG    ; for SPICE, to make node (1) non-floating
RC2 2 0 1
RC3 2 3 0.33333
LC 3 0 0.033333 ; Z2=(s+10)/(s+4)
* plant loop, non inverting, closed loop 0.1s/(s^2+0.2s+4)
EP 5 0 2 4 10
GP1 4 0 5 0 1
RP 4 0 0.5
CP1 4 0 10
LP 4 0 0.025    ; Z5=0.1s/(s^2+0.05s+4)
* plant integrators, (4) to (6) to (7)
GP2 6 0 0 5 80
GP3 7 0 0 6 1
RP2 5 0 1MEG    ; for SPICE
RP3 6 0 1MEG    ; for SPICE
RP4 7 0 1MEG    ; for SPICE
CP2 6 0 1
CP3 7 0 1
* loop closing resistor
RL 7 1 1MEG     ; to close the loop, place semicolon
; after M to reduce RL
*
VIN 8 0 AC 1
RIN 8 1 1
.AC DEC 100 0.1 20
.PROBE         ; (if using PSPICE)
.END
```

Back

## Chapter 4

#### 4.2.3  Example of a system having a loop response with a Bode step

Example 1.  The loop transfer function is plotted with

```n = conv([11 55 110 36], [-1 10]);
d = conv([1 7.7 34 97 83 0 0], [1 10]);
w = logspace(-1,1,200);
bode(n,d,w)
```

The L-plane Nyquist diagram is plotted with

```[mag, phase] = bode(n,d,w);
plot(phase,20*log10(mag),'w', -180,0,'wo'); grid
```

The closed-loop transfer function M = n/g where g = n + d is plotted with

```n = conv([0 0 0 1],n);
g = n + d;
w = logspace(-1,1,200);
bode(n,g,w)
```

The prefilter response is plotted with

```nr = [1 1 0.81]; dr = [1 2 0.81]; bode(nr, dr, w)
```

and the closed-loop response with the prefilter is plotted with

```nc = conv(nr, n); gc = conv(dr, g);
bode(nc, gc, w)
```

By compensating for the hump with a small, 1.6 dB additional notch with

```nr2 = [1 0.5 1]; dr2 = [1 0.6 1];
nc2 = conv(nr2, nc); gc2 = conv(dr2, gc);
bode(nc2, gc2, w)
step(nc2, gc2)
```

The effects of variations in k can be calculated as follows:

```n = conv([11 55 110 36], [-1,10]);
d = conv([1 7.7 34 97 83 0], [1 10]);
k = 1; n = conv(n, k);                 % specify k
w = logspace(-1,1,200);
bode(n, d, w)                          % open-loop response
[mag, phase] = bode(n, d, w);
plot(phase, 20*log10(mag))             % Nyquist diagram
grid
n = conv([0 0 1], n);
g = n + d;
bode(n, g, w)                          % closed-loop response
nr = [1 1 0.81]; dr = [1 2 0.81];
bode(nr, dr, w)                   % prefilter notch response
nc = conv(nr, n); gc = conv(dr, g);
bode(nc, gc, w)          % closed-loop response with notch
nr2 = [1 0.5 1]; dr2 = [1 0.6 1];
bode(nr2 ,dr2, w)                 % second notch response
nc2 = conv(nr2, nc); gc2 = conv(dr2, gc);
bode(nc2, gc2, w)        % closed-loop response with notches
step(nc2, gc2)                    % step response with notches
```

Example 2.  The MATLAB code for this transfer function's numerator and denominator is

```k = 0.7; res = 2; ja = 0.027; a = k*k/(res*ja);
nc1 = conv(res*ja/k,[1 a]); dc1 = [1 0];
```

Example 3.  The numerator `ns` and the denominator `ds` can be found with

```n = conv([11 55 110 36], [-1,10]); wb = 75.4;
d = conv([1 7.7 34 97 83 0 0], [1 10]);
format short e; [ns, ds] = lp2lp(n, d, wb)
```

and the numerators and denominators of the prefilter with

```nr = [1 1 0.81]; dr = [1 2 0.81];
[nrs, drs] = lp2lp(nr, dr, wb)
nr2 = [1 0.5 1]; dr2 = [1 0.6 1];
[nr2s, dr2s] = lp2lp(nr2, dr2, wb)
```

#### 4.3.6  Lightly damped flexible plants; collocated and non-collocated control

Example 2.  The gain and phase responses are plotted in Fig. 4.36 with

```n = conv([1 0 10],[1 0 50]);
d = conv([1 0],[1 0 11]); d = conv(d,[1 0 54]);
w = logspace(0, 1, 1000); bode(n,d,w)
```

### 4.6  Problems

#### Answers to selected problems

4

```fst = 120; Q = 10; fb = fst/2^[20*log10(Q)/18 + 2]
fb = 13.8881
```

Back

## Chapter 5

### 5.2  Asymptotic Bode diagram

Example 1.  The polynomial coefficients can be found from the roots with

```rn = [-0.5 -2 -20];
rd = [-0.2 -0.3 -10];
num = poly(rn);
denum = poly(rd);
```

and the Bode plot can be obtained with

```bode(num, denum)
```

or by using the command `freqs`.

After multiplication of the polynomials in the numerator and denominator:

```n = 4 * 2.8 * conv([1 0.42],[1 2 4])
n = 11.2000   27.1040   54.2080   18.8160
d = conv(conv([1 0.06],[1 1.4]), conv([1 2.8],[1 2.4 9 0]))
d = 1.0000  6.6600  23.3960  48.5880 38.1125  2.1168  0
```

The Nyquist diagram on the L-plane is plotted with

```w = logspace(-1, 1);    [mag, phase] = bode(n, d, w);
plot(phase,20*log10(mag),'r', -180, 0,'wo')
title('L-plane Nyquist diagram')
set(gca,'XTick',[-270 -240 -210 -180 -150 -120 -90])
grid
```

### 5.8  Simulation of a PID controller

The SPICE input file is shown below.

```*** PID example Figs. 5.19, 5.20 ***
ES   3 0 1 2 1            ; input signal summer
RSR1 1 0 1MEG             ; leakage resistor
RSP2 2 0 1MEG             ; leakage resistor
***
GSAT 5 0 0 3 0.001        ; saturation in I-path
RSAT 5 0 1K               ; threshold = (0.7+VT1)*GI1/GSAT
D1   5 6 DIODE
D2   7 5 DIODE
.MODEL DIODE D
VT1  6 0 1V
VT2  0 7 1V
***
GI1  8 0 0 5 1            ; I-path
CI2  9 0 0 8 10           ; integral coefficient
RSP8 8 0 1MEG             ; leakage resistor
***
GP   9 0 0 3 2            ; proportional coefficient
***
GD1  4 0 0 3 3            ; differential coefficient
LD   4 0 1
GD2  9 0 0 4 1
***
RS   9 0 1                ; summing resistor
***
GA1  10 0 9 0 1           ; actuator gain coefficient
***
RSATA 10 0 1K             ; saturation in actuator
D3    10 11 DIODE
D4    12 10 DIODE
VT3   11 0 1V
VT4   0 12 1V
GA2   13 0 0 10 1         ; force source
***
CP1   13 0 5              ; mass of the main body
RSP13 13 0 1MEG           ; leakage resistor
LP2   13 14 0.1           ; spring of flexible mode
CP2   12 0 0.5            ; mass of second body
RP    14 15 0.02          ; losses in the flexible mode
GINT  2 0 0 13 1
CINT  2 0 1               ; integrator to generate position
***   V10 is force, V13 is velocity, V2 is position
***
VTEST 1 0 AC 1            ; use only when frequency responses
; are tested
.AC DEC 20 0.01 10 ; use only when frequency responses tested
** Pulse            (Vmin Vmax delay rise fall width period)
* VPULSE 10 PULSE ( 0V  10V   0S    0S   0S   500    500 )
; when transient responses tested
* .TRAN 0.1 10            ; when transient responses tested
.PROBE                   ; or other graphical postprocessor
.END
```

### 5.10  Digital compensator design

Example 2.  The same problem is solved with MATLAB by

```n = [1 5];
d = [1 10];
fs = 100;                     % sampling frequency in Hz
[nd,dd] = bilinear(n,d,fs)
nd =  0.9762   -0.9286
dd =  1.0000   -0.9048
```

### 5.12  Problems

#### Answers to selected problems

15  MATLAB function `conv` is used to multiply the polynomials in the numerator:

```a = ; b = [1 0.88]; c = [1 3.5 19.4];
ab = conv(a,b); num = conv(ab,c)
```

and in the denominator:

```d = [1 0]; e = [1 0.22]; f = [1 5.5]; g = [1 10.7 157];
de = conv(d,e); def = conv(de,f); deff = conv(def,f);
den = conv(deff,g)
```

16 (1)  This response plotted with MATLAB commands

```n = [10.8 37.9 52.8]; d = [1 6.4 25.6 64 0 0];
w = logspace(-1,1);
bode(n,d,w)
```

With:

```n1 = conv(3.3,[ 1 0.3]);
% during iterations, adjust the zero
% (try  also zero 0.75, pole 1.5)
d1 = conv([1 4 0], [1 1]);
% during iterations, adjust the pole
d2 = [1 2.4 16];
d1d2 = conv(d1,d2);
d = conv(d1d2,[1 0])
n = conv(n1,d2) + conv(7.5,d1)
% vectors have equal length since
% both polynomials are cubic
w = logspace(-1,1);
bode(n,d,w)
```

(2)  With MATLAB, the function can be decomposed:

```num = [601 2632 13510 10260];
den = [1 21.9 309.7 2117.8 5200.4 1044.8];
[Res,Pol,K] = residue(num, den)
Res =                    Pol =                     K = []
1.0e+002 *
0.1598 - 0.2042i     -5.3460 +11.3430i
0.1598 + 0.2042i     -5.3460 -11.3430i
-0.1684 - 4.4410i     -5.4940 + 0.1396i
-0.1684 + 4.4410i     -5.4940 - 0.1396i
0.0172 - 0.0000i     -0.2200
```

The polynomials can be found as follows:

```num_prod = [2*a  (-2*(a*c + b*d))]
den_prod = [1 (-2*c) (c*c + d*d)]
```

When d is the result of calculation inaccuracy, then, approximately,

```num_prod = [2*a  (-2*a*c)]
den_prod = [1 (-2*c) (c*c)]
```

(3)  With:

```num = conv([1 0.88],[ 1 3.5 19.4])
d1 = conv([1 0.22], [1 5.5]);   den = conv(d1, [1 10.7 157])
```

With MATLAB, the function can be decomposed as follows:

```num = [1 4.38 22.48 17.072];
den = [1 16.42  219.414  910.987 189.97];
[Res,Pol,K] = residue(num, den)

Res =                    Pol =                     K = []
0.3889 + 0.2973i     -5.3500 +11.3304i
0.3889 - 0.2973i     -5.3500 -11.3304i
0.2072               -5.5000
0.0151               -0.2200
```

The sum of the two complex pole fractions can be found as follows:

```a = 0.3889; b = 0.2973; c = -5.35; d = 11.3304;
prod_num = [2*a  (-2*(a*c + b*d))]
prod_den = [1 (-2*c) (c*c + d*d)]
```

Back

## Chapter 6

#### 6.2.1  Cauer and Foster RC two-poles

In MATLAB, the response can be found using

```r = [r1 r2 .. ri ..];  p = [p1 p2 .. pi ..];
[num, den] = residue(r,p); bode(num,den)
```

Example 2.  The MATLAB code calculates the residues:

```num = [0.000001 0.0005]; den = [1 2100 200000];
[r, p] = residue(num, den)
```

Back

## Chapter 7

### 7.7  Examples of system modeling

Example 1.

```[A,B,C,D]=linmod('file_name'); bode(A,B,C,D,1)
```

Back

## Chapter 10

#### 10.7.3  Design examples

Example 1.  The logarithmic Nyquist diagram for TP is plotted in Fig. 10.24(b) with the `nyqlog` function from the Bodestep toolbox described in Appendix 14:

```np = [2 1]; dp = [1 2 0 0]; nyqlog(1,np,dp)
```

(A)  The L-plane Nyquist diagrams are plotted with the following script:

```np = [2 1]; dp = [1 2 0 0]; nyqlog(1,np,dp); hold on
ne = 2; de = [1 2 0]; nyqlog(1,ne,de); hold on
ng = 1; dg = [1 2 2 0]; nyqlog(0.5,ng,dg); zoom off
gtext('TP');  gtext('TE');  gtext('G') %place labels with mouse
```

Example 2.

(C)

```npg = [0 0 2 1.2 0.1]; dpg = [1 2 1.6 0.16 0];
gpg = npg + dpg; sstep(npg,gpg); grid

ne = [0.4 1.04 0.1]; de = [1 2 1.6 0.16 0];
ge = n + d; sstep(n,g); grid
```

Back

## Chapter 11

### 11.7  NDC with a single nonlinear nondynamic link

Example 1.  The SPICE simulation input file is shown below.

```**** ch9ex1.cir for iso-w simulation of NDC Figs. 11.21(b), **** 11.22
*** input integrator
G2 2 0 0 1 1
C2 2 0 1
R2 2 0 1MEG
*** feedback summer
G3 3 0 7 2 1
R3 3 0 1
***  kDF path:  G5 =.1, 1, or 10
G5 5 0 0 3 10
R5 5 0 1
*** forward path, inverting
G4 4 0 3 0 1
C4 4 0 1
R4 4 0 1MEG
*** forward summer, the output is VDB(6)
G6 6 0 4 5 1
R6 6 0 1
*** feedback path
G7 7 0 0 5 1
C7 7 0 1
R7 7 0 1MEG
***
VIN 1 0 AC 1
RIN 1 0 1MEG
.AC DEC 20 .001 10
.PROBE
.END
```

#### Answers to selected problems

1(a)  When employing SPICE, the input file is

```*ch9oc1.cir for determining oscillation condition
*T = 200(s+300)(s+600)/[s(s+20)(s+50)(s+100)], open loop
G2 2 0 0 1 200  ; gain, zero at 300
L2 2 8 1
G3 3 0 0 2 1    ; zero at 600
L3 3 9 1
R3 9 0 600
G4 4 0 0 3 1    ; pole at 20
R4 4 0 0.05
C4 4 0 1
G5 5 0 0 4 1    ; pole at 50
R5 5 0 0.02
C5 5 0 1
G6 6 0 0 5 1    ; pole at 100
R6 6 0 0.01
C6 6 0 1
G7 7 0 0 6 1    ; integrator
R7 7 0 1MEG     ; de-floating resistor
C7 7 0 1
VIN 1 0 AC 1
RIN 1 0 1MEG    ; source loading resistor
.AC DEC 20 1 100
.PROBE
.END
```

Back

## Appendix 2

### A2.4  Laplace transfer function

Example 2.  The time-function can be plotted either with

```num = [5 7];
den = [1 3 2 0];
impulse(num, den)
```

or with

```num = [5 7];
den = [1 3 2];
step(num, den)
```

Back

## Appendix 5

The programs' listing follows:

Function `find_phase2`

```function [phase] = find_phase2(magdb, natfreq)
% function [phase] = find_phase2(magdb ,freq)
% This routine uses the Bode Integral to generate phase data.
% from a magnitude response of a m.p. function.
% magdb:  row gain vector given in dB
% freq:   row frequency vector given in rad/sec
% The Magnitude and Frequency vectors must be the same length.
% Before running this function prepare table with table_maker.

[row,col] = size(magdb); if col == 1, magdb = magdb'; end;
[row,col] = size(natfreq); if col == 1, natfreq = natfreq'; end;

%%% table_maker % calls the function creating the table
load table  % load data needed for toe and tail calculations:
% table con hilimit lolimit numentries numintstep
points = length(natfreq);
numsteps = points - 1;

% The following variables are for the lookuptable (u domain)
ilnfreq = log(natfreq(1));
flnfreq = log(natfreq(points));
toeslope = (magdb(2) - magdb(1))/(log(natfreq(2))-ilnfreq);
tailslope = (magdb(points) - magdb(numsteps))/(flnfreq -log(natfreq(numsteps)));

dmagdb = magdb(2:points) - magdb(1:numsteps);
dnatfreq = natfreq(2:points) - natfreq(1:numsteps);

deriv = dmagdb./dnatfreq;

nnfreq = natfreq(1:numsteps);
w1 = nnfreq';w2 = natfreq(2:points)';

for I = 1:points,

% The next lines perform the integration...
wc = natfreq(i)*ones(w1);
weights = (w2+ wc).*log(w2 + wc) - (w1 + wc).*log(w1 + wc) +...
(wc-w2).*log(abs(wc-w2+eps)) - (wc-w1).*log(abs(wc-w1+eps));

looplnfreq = log(natfreq(i));
u = [flnfreq ilnfreq] - [looplnfreq]*[1 1] + [eps eps];
ind = (log10(abs(u))+[-log10(lolimit)-log10(lolimit)])/con+[1 1];
ind = max([1 1;ind]);ind=min([numentries numentries;ind]);
tailtoe = abs(pi^2/4*[1;1]-table(ind',2));

phase(i) = (deriv*weights + [tailslope toeslope]*tailtoe)/pi;
end
```

Function `table_maker.m`

```function table_maker
% table_maker.m, called by function find_phase2. This routine
% creates the lookup table for calculation the toe and tail
% 1e-15 is used for zero, 650 is used for infinity
% reasonable limits:  1e-9 to 100
% This contains nearly all of the area (to 6 places).

lolimit = 1e-9;  hilimit = 100;
numentries = 4001;  % number of table entries, should be odd
numintstep = 100;   % number of steps for each integration

vector = logspace(log10(lolimit),log10(hilimit),numentries);
table(:,1) = vector';

[table(1,2)] = integral_u(1e-15,lolimit, numintstep);

for k = 2:length(vector),
[table(k,2)]=table(k-1,2)+integral_u(table(k-1,1),table(k,1), numintstep);
end

clear vector k

con = log10(table(7,1))-log10(table(6,1));

save table con hilimit lolimit numentries numintstep table
clear con hilimit lolimit numentries numintstep table
```

Function `integral_u`

```function [trap]=integral_u(u1,u2,numsteps)
% function [trap]=integral_u(u1,u2,numsteps)
% This routeens integrates in u domain.
% It may not handle some special cases, however.
% vector=linspace calculation(u1,u2,numsteps+1);

if (u1==0),
u1=1e-15;
end
if (u2==0),
u2=1e-15;
end
vector = logspace(log10(u1),log10(u2),numsteps+1);
% This next calculation is the logu_function

valvector = log(abs((exp(vector)+1)./(exp(vector)-1)));
delta = vector(2:numsteps+1)-vector(1:numsteps);

%up =   valvector(1:numsteps)*delta';
%down = valvector(2:numsteps+1)*delta';

trap=(.5*(valvector(1:numsteps)+valvector(2:numsteps+1)))*delta';
if trap==NaN,
trap=0;
disp('Warning:  NaN found in integration, result set to 0')
end
```

Back

## Appendix 13

### A13.12  Conceptual design of an antenna attitude control

The Bode diagram is plotted with

```w = logspace(1,3);
[A,B,C,D] = linmod('mlsplan');
bode(A,B,C,D,1,w)
hold on
% remove the connection to the input of block Z2
bode(A,B,C,D,1,w)
```

### A13.13  Pathlength control of an optical delay line

#### Block diagram with rational transfer functions

The Bode diagrams are obtained with the following script:

```npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28];
dpzt = [1  2.83e5  7e9  9.18e13  6.5e17  1.684e21  0  0];
nvc = [2.5e9]; dvc = [1 40 400 0];

ntot1 = conv(npzt,dvc);
ntot2 = conv(dpzt,nvc);
ntot = ntot1 + [adz ntot2];
dtot = conv(dvc,dpzt);

w = logspace(1,4.5);
bode(npzt,dpzt,w);
hold on
bode(nvc,dvc,w);
bode(ntot,dtot,w);
hold off
zoom on
```

The logarithmic plane Nyquist diagrams are plotted with `nyqlog`:

```nyqlog(3700, nvc,dvc)
hold on
nyqlog(3700, npzt,dpzt)
hold on
nyqlog(3700, ntot,dtot)
hold off
```

#### Limit cycle of the VC loop with the NDC

The Bode diagram is obtained with the script

```npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28];
dpzt = [1  2.83e5  7e9  9.18e13  6.5e17  1.684e21  0  0];
nvc = [2.5e9]; dvc = [1 40 400 0];

n1 = conv(npzt,dvc);
d1 = conv(dvc,dpzt);
d11 = [0 n1] + d1;

d2 = conv(dpzt,nvc);
deq = d11 + [0 0 0 d2];

w = logspace(0,4);
% nyqlog(1260,-n1,deq)
bbode(n1,deq,w)
```

#### Global stability of a system with the NDC

The Bode diagram of this transfer function found with the script

```npzt = [2.33e-10 6.7e-6 -5.14e11 1.26e17 1.5e21 9.5e24 1.2e28];
dpzt = [1  2.83e5  7e9  9.18e13  6.5e17  1.684e21  0  0];
nvc = [2.5e9]; dvc = [1 40 400 0];

n1 = conv(nvc,dpzt);
d1 = conv(dvc,dpzt);

d2 = conv(npzt,dvc);
deq = d1 + [0 d2];

w = logspace(1,4);
bode(n1,deq,w)
```

Back