PROGRAM d5r4(input,output);
(* driver for routine CHEBFT *)
LABEL 10,99;
CONST
   nval=40;
   pio2=1.5707963;
   eqs=1.0e-6;
TYPE
   glcarray=ARRAY [1..nval] OF real;
VAR
   a,b,dum,f : real;
   t0,t1,term,x,y : real;
   i,j,mval : integer;
   c : glcarray;

FUNCTION func(x: real): real;
BEGIN
   func := sqr(x)*(sqr(x)-2.0)*sin(x)
END;

(*$I MODFILE.PAS *)
(*$I CHEBFT.PAS *)

BEGIN
   a := -pio2;
   b := pio2;
   chebft(a,b,c,nval);
(* test result *)
10:   writeln;
   writeln('How many terms in Chebyshev evaluation?');
   write('Enter n between 6 and ',nval:2,
         '. (n := 0 to end).  ');
   readln(mval);
   IF ((mval <= 0) OR (mval > nval)) THEN GOTO 99;
   writeln;
   writeln('x':9,'actual':14,'chebyshev fit':16);
   FOR i := -8 to 8 DO BEGIN
      x := i*pio2/10.0;
      y := (x-0.5*(b+a))/(0.5*(b-a));
(* evaluate chebyshev polynomial without using routine chebev *);
      t0 := 1.0;
      t1 := y;
      f := c[2]*t1+c[1]*0.5;
      FOR j := 3 to mval DO BEGIN
         dum := t1;
         t1 := 2.0*y*t1-t0;
         t0 := dum;
         term := c[j]*t1;
         f := f+term
      END;
      writeln(x:12:6,func(x):12:6,f:12:6)
   END;
   GOTO 10;
99:
END.