%code2_1 for lecture 1; Simpson's rule
N=input(' N= ');
eps=input( 'eps= ');
T=input(' T= ');
hh=2/N;
guo=1.0/(eps*T+1);
for k=1:2:21
ss=0;
xx0=-1;
tt0=(-1)^k;
va0=guo*eps*(xx0+tan(xx0*0.5*guo));
while xx0 < 1
xx1=xx0+hh;
tt1=chebyshev(k,xx1);
va1=guo*eps*(xx1+tan(xx1*0.5*guo));
xx2=xx1+hh;
tt2=chebyshev(k,xx2);
va2=guo*eps*(xx2+tan(xx2*0.5*guo));
ss=ss+hh*(va0*tt0+4*va1*tt1+va2*tt2)/3;
tt0=tt2;
xx0=xx2;
va0=va2;
end
guo=xx0
xx(k)=k;
yy(k)=abs(ss);
end
fprintf(1, '%9.0f %10.4e ', [xx; yy]);
