(*Adapted by Rafael Luís from Sebastian M.Marotta
updated to Mathematica v.10 by Ross Moore *)

g[z_]:=(Sqrt[z-c])
h[z_]:=-(Sqrt[z-c])
Equation=Function[z,If[Random[]>0.5,g[z],h[z]]];
c=-0(*+0.0I*);
Iterates=50000;
DropIterates=100;
Orbit=Drop[NestList[Equation,Random[],Iterates],DropIterates];
DataPlot=Table[{AbsolutePointSize[.1],Point[{Re[Orbit[[i]]],Im[Orbit[[i]]]}]},{i,1,Length[Orbit]}];
$TextStyle={FontSize->12};
Graphics[DataPlot,Frame->True,Axes->False,AspectRatio->Automatic,FrameLabel->{"Re(z)","Im(z)"},RotateLabel->False,PlotLabel-> TraditionalForm["Julia set for  f(z) = (z^2)+c, c = "<>ToString[N[c]]],ImageSize->400]