%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  \n', [xx;  yy]);