Knowing how other CAS do this might help you.
To my knowledge, Mathematica uses a variation on the following basic algorithm for plotting a one-variable function $f(x)$ or a parametric curve $(x(t), y(t))$ (I'm going to assume $f(x)$ for this description).
Start with a regularly spaced grid of points on the plotting domain. (In Mathematica, there's a parameter to control how many to take, called PlotPoints
.)
Look at each pair of successive line segments (defined by three successive points, $(x_1, f(x_1)),\; (x_2, f(x_2)),\; (x_3, f(x_3))$) and insert a new sampling point in the middle of both segments ($\frac{x_1+x_2}{2}$ and $\frac{x_2+x_3}{2}$) if their angle is larger than a threshold.
If we haven't reached the iteration limit yet (set by MaxRecursion
in Mathematica), repeat from step 2.
Some of this is discussed in the book Mathematica in Action by Stan Wagon, which you can see here on Google Books.
I implemented this algorithm before to have better control over how many times my expensive to calculate function was evaluated. Here's the Mathematica code for step 2:
nd[{points_, values_}] :=
Transpose@{(Drop[points, 1] + Drop[points, -1])/2,
Differences[values]/Differences[points]}
subdivide1d[result_, resolution_, maxAngle_: 10] :=
Module[
{deriv, angle, dangle, pos, nf},
deriv = nd[result\[Transpose]];
angle = ArcTan[#2] & @@@ deriv;
dangle = Differences[angle];
pos = Flatten@Position[dangle, d_ /; Abs[d] > maxAngle/180 Pi];
pos = Union[pos, pos + 1];
nf = Nearest[result[[All, 1]]];
Select[deriv[[pos, 1]], Abs[# - First@nf[#]] > resolution &]
]