The Method of Graphical Analysis

(* Dr. Koch's Graphical Analysis procedure *)
(* Modified by John to act more like his Maple Code *)
(* Modified by Marcin Bownik to include Manipulate procedure *)
(* Ga[F, {a,b}, n | {n, m}, x0 | {x0, x1, x2, ...}] *) 

(* F should be some function of one variable *)
(* {a,b} is the domain of the function to be plotted *)
(* n is the number of steps to plot. If {n, m} is given instead, m steps will be thrown away, then n steps will be plotted *)
(* x0 is the point to start at. If a list of points is given (ie {x0, x1, x2, ...}), then each will be plotted in a different color. *)

Ga[F_, bound_, steps_, start_] := Module[{g1, g2, g3, g4, l, temp, index, slist, j, n, m},
	(* Generate graph of the function and F[x] == x line *)
	
			  g1 = Plot[{F[x], x},{x, bound[[1]], bound[[2]]},PlotRange->{bound[[1]], bound[[2]]},AspectRatio->1,  DisplayFunction->Identity];
(* Parse steps *)
  If[Length[steps] == 0, n= steps; m=0;, n = steps[[1]]; m = steps[[2]]; ]; 
(* Parse the starting points*)
  	l = Length[start];   
		If[l == 0, slist = {start};l=1;, slist = start];
		g2 ={};
		(* Loop through starting points *)
		For[j = 1, j <= l, j = j + 1, 
		temp = slist[[j]];
			   g2= Append[g2,Hue[(j+1)/(l+1)]];
	(* Loop through the first m iterations, and throw them away *)
For[index =0, index < m, index = index + 1,temp = F[temp];];
			
	(* Draw the lines (x,x) to (x, F[x]) to (F[x], F[x]) *)	
		For[index = 0, index < n, index = index + 1,
				g3 = Append[g2,Line[{{temp, temp}, {temp, F[temp]}}]];
				g2 = Append[g3,Line[{{temp, F[temp]},{F[temp], F[temp]}}]];
				temp = F[temp];
				];
g2= Append[g2,{PointSize[Medium],Point[{temp,temp}]}];
];
	(*	Print[g2]; *)
			  
				g4 = Graphics[g2];
		
		Show[g1, g4,  PlotRange->{{bound[[1]], bound[[2]]}, {bound[[1]], bound[[2]]}},DisplayFunction->$DisplayFunction]
				]; 

GraphOrbit[F_,bound_, endstep_, start_]:=
Manipulate[ Ga[F , bound , n , start], {n,0,endstep,1}]

GraphOrbit[Sqrt,{0,1.2}, 20, 0.1]

GraphOrbit[(#1^3+#1)/6-1&, {-3, 5}, 10, 2.5]

GraphOrbit[(#1^3+#1)/6-1&, {-3, 10}, 4, 2.8]

GraphOrbit[#1^2-1&, {-2, 2}, 10, 0]

GraphOrbit[#1^2-1&, {-2, 2}, 20, 0.2]

GraphOrbit[#1^2-1&, {-2, 2}, 20, 0.3]

(* Here is the quadratic map at c = -2 *)
f[x_] := x^2 - 2
GraphOrbit[f, {-2,2}, 100, {0, 1.55, -1, -.5}]

(* Here it is again using fancier code *)
(* Note that #1^2 - 2&, is just x^2 - 2 *)
(* The #1 is the variable of the function, and the & at the end tells mathematica to treat this as a function, instead of as data. *)

GraphOrbit[#1^2  -2 &,  {-2, 2}, 40, {0, 1.55, -1, -.5}]

(* This is using Nest[] to plot the 3rd iterate of f *)
(* Nest[f, #1, n]& would be the nth iterate of f *)

GraphOrbit[Nest[f, #1, 3]&, {-2, 2}, 20, {0, .5}]

(* This is the third iterate again, using fancy code for x^2 -2.  I like this way better, since it is self contained *)

GraphOrbit[Nest[#1^2  -2  &, #1, 3]&,  {-2, 2}, 20, { 0,.5}]

GraphOrbit[f, {-2,2}, 100, 0.2]

Manipulate[ Ga[#1^2  +c  & , {-3, 3} , n , 0.5 ], {n,0,30,1}, {c,1,-2, -0.1}]

GraphOrbit[#1^2  -1.1 &,  {-2, 2}, 10,Table[-1.0 + n/5, {n, 0, 10}]]