BeginPackage["CSM`Iterate`"]

FunctionIteration::usage =
	"FunctionIteration[f, x0, pre, iter, {xmin, xmax}, opts..]
	shows iter iterations of f with start value x0 (which can be a list
	of values). First pre iterations are computed withough plotting."

LogisticsMap::usage = "LogisticsMap[r] is the map x -> r x(1-x)."

NewtonIteration::usage = "NewtonIteration[f, x0, iter, {xmin, xmax}, opts...]
	shows the Newton iteration of a function f."

Begin["`Private`"]

FunctionIteration[f_, x0_, initial_:0, length_, {xmin_, xmax_}, opts___] :=
    Module[{start, orbits, plot, lines, x},
        start = Nest[f, N[x0], initial];
        If[ Head[start] =!= List, start = {start} ];
        orbits = Transpose[ NestList[f, start, length] ];
        plot = Plot[f[x], {x, xmin, xmax}, DisplayFunction->Identity];
        lines = Line[Partition[
                  Flatten[Transpose[{#, #}], 1], 2, 1]]& /@ orbits;
        Show[ plot,
              Graphics[{PointSize[.02],
                        lines, Point[{#, #}]& /@ start,
                        Line[{{xmin, xmin}, {xmax, xmax}}]}],
              opts,
              AxesOrigin->{xmin, xmin},
              DisplayFunction->$DisplayFunction,
              PlotRange->{{xmin, xmax}, {xmin, xmax}},
              AspectRatio->1
        ]
    ]

LogisticsMap[r_] := r#(1-#)&

NewtonIteration[f_, x0_, length_, {xmin_, xmax_}, opts___] :=
    Module[{start = x0, xvals, plot, lines, x, fp = f'},
        If[ Head[start] =!= List, start = {start} ];
        xvals = Transpose[ NestList[# - f[#]/fp[#]&, start, length] ];
        plot = Plot[f[x], {x, xmin, xmax}, DisplayFunction->Identity];
        lines = Line[Partition[Flatten[{#, 0, #, f[#]}& /@ #, 1], 2]]& /@ xvals;
        Show[ plot,
              Graphics[{Thickness[.0001],  PointSize[.02],
                        lines, Point[{#, 0}]& /@ start}],
              opts,
              DisplayFunction->$DisplayFunction
        ]
    ]

End[]

EndPackage[]
