%
% Example 3.19 and Figure 3.14
%
t=0:0.01:4*pi;
E=5; w0=1;
% first approximation
N=1;
yN=0;
for n=1:N
    Hm=1/sqrt((3-(n*w0)^2)^2+4*(n*w0)^2);
    Hp=-atan(2*n*w0/(3-(n*w0)^2));
    Xn=j*(-1)^n*E/(n*pi);
    Xnm=abs(Xn);
    Xnp=angle(Xn);
    Ynm=Xnm*Hm;
    Ynp=Xnp+Hp;
    yN=yN+2*Ynm*cos(n*w0*t+Ynp);
end
plot(t,yN,':')
hold
% second approximation
N=2;
yN=0;
for n=1:N
    Hm=1/sqrt((3-(n*w0)^2)^2+4*(n*w0)^2);
    Hp=-atan(2*n*w0/(3-(n*w0)^2));
    Xn=j*(-1)^n*E/(n*pi);
    Xnm=abs(Xn);
    Xnp=angle(Xn);
    Ynm=Xnm*Hm;
    Ynp=Xnp+Hp;
    yN=yN+2*Ynm*cos(n*w0*t+Ynp);
end
plot(t,yN,'--')
% third approximation
N=3;
yN=0;
for n=1:N
    Hm=1/sqrt((3-(n*w0)^2)^2+4*(n*w0)^2);
    Hp=-atan(2*n*w0/(3-(n*w0)^2));
    Xn=j*(-1)^n*E/(n*pi);
    Xnm=abs(Xn);
    Xnp=angle(Xn);
    Ynm=Xnm*Hm;
    Ynp=Xnp+Hp;
    yN=yN+2*Ynm*cos(n*w0*t+Ynp);
end
plot(t,yN,'-.')
% fourth approximation
N=4;
yN=0;
for n=1:N
    Hm=1/sqrt((3-(n*w0)^2)^2+4*(n*w0)^2);
    Hp=-atan(2*n*w0/(3-(n*w0)^2));
    Xn=j*(-1)^n*E/(n*pi);
    Xnm=abs(Xn);
    Xnp=angle(Xn);
    Ynm=Xnm*Hm;
    Ynp=Xnp+Hp;
    yN=yN+2*Ynm*cos(n*w0*t+Ynp);
end
plot(t,yN)
xlabel('time')
ylabel('approximations, y1, y2, y3, and y4') 
print -deps fig3_14.eps