```\\ adapted from an original idea by Ilya Zakharevich

\\ generate an RGB color triple from a "magnitude" between 0 and 255
\\ (low = close to a cold blue, high = close to a hot red).
\\ To generate simple colormaps.
rgb(mag) =
{ my(x = mag/255., B, G, R);
B = min(max(4*(0.75-x),     0), 1);
R = min(max(4*(x-0.25),     0), 1);
G = min(max(4*abs(x-0.5)-1, 0), 1);
return (floor([R, G, B]*255));
}
default(graphcolormap, concat(["white","black","blue"], vector(25,i,rgb(10*i))));
default(graphcolors, vector(25,i,i+2));

\\ plot Taylor polynomials of f (sin by default),
\\ of index  first + i*step <= ordlim, for x in [xmin,xmax].
plot_taylor(f = x->sin(x), xmin=-5, xmax=5, ordlim=16, first=1, step=1) =
{
my(T,s,t,w,h,dw,dh,cw,ch,gh, extrasize = 0.6);
my(Taylor_array);

default(seriesprecision,ordlim+1);
T = f('q);
ordlim = (ordlim-first)\step + first;
Taylor_array = vector(ordlim+1);
forstep(i=ordlim+1, 1, -1,
T += O('q^(1 + first + (i-1)*step));
Taylor_array[i] = truncate(T)
);

t = plothsizes();
w=floor(t[1]*0.9)-2; dw=floor(t[1]*0.05)+1; cw=t[5];
h=floor(t[2]*0.9)-2; dh=floor(t[2]*0.05)+1; ch=t[6];

plotinit(2, w+2*dw, h+2*dh);
plotinit(3, w, floor(h/1.2));
\\ few points (but Recursive!), to determine bounding box
s = plotrecth(3, x=xmin,xmax, f(x),
"Recursive|no_X_axis|no_Y_axis|no_Frame", 16);
gh=s[4]-s[3];

plotinit(3, w, h);
plotscale(3, s[1], s[2], s[3]-gh*extrasize/2, s[4]+gh*extrasize/2);
plotrecth(3, x=xmin,xmax, subst(Taylor_array, 'q, x), "no_Rescale");
plotclip(3);
plotcopy(3, 2, dw, dh);

plotmove(2, floor(dw+w/2-15*cw), floor(dh/2));
plotstring(2, "Multiple Taylor Approximations");
plotdraw([2, 0, 0]);
}

\p9
plot_taylor()
plot_taylor(x->exp(x),-3,3)
plot_taylor(x->besselk(2,x), 1,5)
plot_taylor(x->1/(1+x^2), -1.2,1.2)
```