edmap:=proc(f,L,a,max)
 #
 #This is a procedure that provides
 #a geometric interpretation of delta
 #and epsilon values for a function
 #near a point.
 #
 # f is the function (expression in x).
 #
 # L is the expected limit of the function
 #     at the value a (real).
 #
 # a denotes the point (real).
 #
 # max is the maximum acceptable epsilon (real).
 #
 local fn,acc,A,B,XL,XR,E,con,j,k,n,xstep;
 with(plots):with(plottools):
 #
 #Set number of frames
 fn:=16:
 #
 #Set number of steps for approximation of delta
 acc:=100:
 xstep:=max/acc:
 #
 #Sequence of epsilon values
 E:=seq(n*max/fn, n=1..fn):
 #Function to convert numbers to strings for textplots
 con:=n->convert(evalf(n,4),string):
 #
 #
 #Sequence of left endpoints
 XL[0]:=a:
 j:=1:
 for n from 1 to fn do
      XL[n]:=XL[n-1]:
      for k from j to acc
           while eval(abs(subs(x=XL[n]-xstep,f)-L))<=E[n] do
           XL[n]:=XL[n]-xstep:
           j:=j+1:
      od:
 od:
 #
 #Sequence of right endpoints
 XR[0]:=a:
 j:=1:
 for n from 1 to fn do
      XR[n]:=XR[n-1]:
      for k from j to acc
           while eval(abs(subs(x=XR[n]+xstep,f)-L))<=E[n] do
           XR[n]:=XR[n]+xstep:
           j:=j+1:
      od:
 od:
 #
 #Declare array of plots
 A:=array(1..2,1..2):
 #
 #Create animation for upper left panel
 B:=seq(display([
      plot([[a-max,L-E[n]],[a+max,L-E[n]]],color=green),
      plot([[a-max,L+E[n]],[a+max,L+E[n]]],color=green),
      plot([[XL[n],L-max],[XL[n],L+max]],color=blue),
      plot([[XR[n],L-max],[XR[n],L+max]],color=red)]),n=1..fn):
 #
 #Complete upper left panel (graph of f with animated lines)
 A[1,1]:=display({
      plot(f,x=XL[fn]..XR[fn],color=black),
      plot([[a,L-max],[a,L+max]],color=black),
      plot([[a-max,L],[a+max,L]],color=yellow),
      display(B,insequence=true)},
      axes=framed):
 #
 #Create upper right panel (text description of upper graph)
 A[1,2]:=display([seq(
        textplot([[0,1,`If |x-`.(con(a)).`|<`.(con(min(a-XL[n],XR[n]-a))).`, then`],
                 [0,.9,`x is in the interval (`.(con(XL[n])).`,`.(con(XR[n])).`),`],
                 [0,.8,`and |f(x)-`.(con(L)).`|<`.(con(E[n]))],
                 [0,.6,`Center of the y-axis:  `.(con(L))],
                 [0,.5,`Green to yellow:  `.(con(E[n]))],
                 [0,.4,`Center of the x-axis: `.(con(a))],
                 [0,.3,`Blue vertical line:  `.(con(XL[n]))],
                 [0,.2,`Red vertical line:  `.(con(XR[n]))],
                 [0,.1,`Minimum of |`.(con(a)).`-`.(con(XL[n])).`| and`],
                 [0,0,`|`.(con(a)).`-`.(con(XR[n])).`| is  `
                                .(con(min(a-XL[n],XR[n]-a)))]],
                 align=LEFT, font=[TIMES,ROMAN,12]),
                 n=1..fn)],insequence=true,axes=none):
 #Create lower left panel (delta as a function of epsilon)
 A[2,1]:=display({
           plot([seq([E[k],a-XL[k]],k=1..fn)],color=blue),
           plot([seq([E[k],XR[k]-a],k=1..fn)],color=red),
           display([seq(plot([[E[n],0],[E[n],max]],color=black),n=1..fn)]
                  ,insequence=true,labelfont=[SYMBOL,18],labels=[e,d])},
           labelfont=[SYMBOL,18],labels=[e,d],axes=framed,xtickmarks=2):
 #Create lower right panel (text description of lower graph)
 A[2,2]:=display([seq(
       textplot([[0,.8,`Delta as a function of epsilon...`],
                [0,.7,`Black vertical line:  `.(con(E[n]))],
                [0,.6,`Blue line height:  |`.(con(a)).`-`.(con(XL[n])).`|=`.(con(a-XL[n]))],
                [0,.5,`Red line height:  |`.(con(XR[n])).`-`.(con(a)).`|=`.(con(XR[n]-a))]],
                align=LEFT,font=[TIMES,ROMAN,12]),n=1..fn)],insequence=true,axes=none):
 display(A,insequence=true);
 end: