`help/text/poly` := TEXT(
`poly	-- return 3rd degree polynomial satisfying the specified `,
`	end conditions`):

# ==================================================================
# Rational approximation on [a,infinity); give arguments in the form 
# (ffrom, which), where which = 1,2,3 or 4, and 
# 1 = 1,0,0,0
# 2 = 0,1,0,0  
# 3 = 0,0,1,0
# 4 = 0,0,x,1

a2inf := proc(ffrom, which)
   local c1;
   if which =1 then;
      # 	case 1:
      c1 := subs(x=x-ffrom, (x+1)/(x^2+x+1));
   elif which =2 then;
      # 	case 2:
      c1 := subs(x=x-ffrom, x/(x^2+x+1));
   elif which =3 then;
      # 	case 3:
      c1 := subs(x=x-ffrom, x^2/(x^2+x+1))
   elif which =4 then;
      # 	case 4:
      c1 := subs(x=x-ffrom, x^2/(x+1));
   else;
      ERROR(`invalid type argument`);
   fi;
   c1;
end:

# 	Tested with the following:
# > t1 := a2inf(a, 4);
# > subs(x=a, t1);
# > subs(x=a, diff(t1, x));
# > limit(t1/x, x=infinity);
# > limit(diff(t1, x), x=infinity);
# ==================================================================

poly := proc(lval,ldval, rval, rdval, ffrom, tto)
   local p,k, c, sy1,sy2,sy3,sy4,pp, t1, t2;
   p := sum('x**k*c.k', k=0..3);
   pp := diff(p, x);
   sy1 := eval(subs(x=ffrom, p)) = lval;
   sy2 := eval(subs(x=ffrom, pp)) = ldval;
   sy3 := eval(subs(x=tto, p)) = rval;
   sy4 := eval(subs(x=tto, pp)) = rdval;
   t1 := solve({sy1,sy2,sy3,sy4}, {c0,c1,c2,c3});
   eval(subs(t1, p));
end:		# ok [Nov-06-1994]


FF := proc(gpar)	# THE ARG TO FF MUST BE IN TERMS OF X !
   local i, j, a, m, indi,g;
   m := M+N+1;
   a := matrix(m,1,0);
   g := unapply(gpar, x);
   for i from -M to N do;
      indi := i + M +1;
      a[indi,1] := evalf(g(psi(i*h)));	# evaluate g at the node points
   od;
   op(a);
end:		# ok [Nov-06-1994]

DD := proc(gpar)	# THE ARG TO DD MUST BE IN TERMS OF X !
   local i, j, a, m, indi,g;
   m := M+N+1;
   a := matrix(m,m,0);
   g := unapply(gpar, x);
   for i from -M to N do;
      indi := i + M +1;
      a[indi,indi] := evalf(g(psi(i*h)));	# evaluate g at the node points
   od;
   op(a);
end:		# ok [Nov-06-1994]


delta := proc(num, j,k);
   if (num=0) then;
      if (j=k) then;
	 1;
      else;
	 0;
      fi;
   elif (num=1) then;
      if (j=k) then;
	 0;
      else;
	 (-1)**(k-j)/(k-j);
      fi;
   elif (num=2) then;
      if (j=k) then;
	 -Pi**2/3;
      else;
	 -2*(-1)**(k-j)/(k-j)**2;
      fi;
   fi;
end:   		# ok [Nov-06-1994]
   

sinc := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(1-1/6*Pi^2*x^2+1/120*Pi^4*x^4);
      else;
	 evalhf(sin(Pi*x)/(Pi*x));
      fi;
   else;
      'sinc'(y);
   fi;
end:

`diff/sinc` := proc(a,x);
   sinc_p(a)*diff(a,x);
end:

sinc_p := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(-1/3*Pi^2*x+1/30*Pi^4*x^3-1/840*Pi^6*x^5);
      else;
	 evalhf(cos(Pi*x)/x-sin(Pi*x)/Pi/x^2);
      fi;
   else;
      'sinc_p'(y);
   fi;
end:

`diff/sinc_p` := proc(a,x);
   sinc_pp(a)*diff(a,x);
end:

sinc_pp := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(-1/3*Pi^2+1/10*Pi^4*x^2-1/168*Pi^6*x^4);
      else;
	 evalhf(-sin(Pi*x)*Pi/x-2*cos(Pi*x)/x^2+2*sin(Pi*x)/Pi/x^3);
      fi;
   else;
      'sinc_pp'(y);
   fi;
end:

plot_func := proc(a,b);		# provide for plotting of function
   plot({u(x),exa(x)}, x=a..b, title = `function plotted at random points`.` M= `.M.` N = `.N);
end:

plot_func_at_gridpoints := proc(a,b)
   local i, xi, t1,t2,plt, plt1;
   plt := [];
   plt1 := [];
   for i from -M-1 to N+1 do;
      xi := psi(i*h);
      if (a<xi and xi<b) then;
	 plt := [op(plt), [xi, exa(xi)]];	# exact answer
	 plt1 := [op(plt1), [xi, u(xi)]];	# approximation
      fi;
   od;
   plot({plt, plt1}, style = point, title = `function plotted at gridpoints`.` M= `.M.` N = `.N);
end:
      
plot_der1_at_gridpoints := proc(a,b)
   local i, xi, t1,t2,plt, plt1;
   plt := [];
   plt1 := [];
   t1 := unapply(diff(exa(x),x),x);
   t2 := unapply(diff(u(x),x),x);
   for i from -M-1 to N+1 do;
      xi := psi(i*h);
      if (a<xi and xi<b) then;
	 plt := [op(plt), [xi, t1(xi)]];	# exact answer
	 plt1 := [op(plt1), [xi, t2(xi)]];	# approximation
      fi;
   od;
   plot({plt, plt1}, style = point, title = `derivative plotted at gridpoints, M =  `.M.` N = `.N);
end:

plot_der_1 := proc(a,b)		# plot derivatives
   local t1, t2;
   t1 := diff(exa(x),x);
   t2 := diff(u(x),x);
   plot({t1, t2}, x=a..b, title = `derivative plotted at random points`.` M= `.M.` N = `.N);
end:
   
plot_der_2 := proc(a,b)		# plot derivatives
   local t1, t2;
   t1 := diff(exa(x),x,x);
   t2 := diff(u(x),x,x);
   plot({t1, t2}, x=a..b);
end:

macro(pl= plot_func, pl1 = plot_der_1, pl2= plot_der_2, plg = plot_func_at_gridpoints , plg1 = plot_der1_at_gridpoints );