POV Ray 3.7 code
//=================================================================
// Newton's method 3D set ($f(z) = z^{4} + 1$) generated randomly
// (POV-Ray cylindrical symmetry)
//=================================================================

#include "colors.inc"  // Standard Color definitions
#include "glass.inc"   //  Glass pigment effect
#include "metals.inc"   // Metal pigment effect
global_settings {assumed_gamma 1.0}   // set display gamma
background { color rgb <00,0.0,0.0> }    // set black background

camera {							// set view point (camera) location
    location <-20, 20, -450>
    look_at  <-35, 1,  0>
}

  light_source    // Area light source (creates soft shadows)
{
  0*x // light's position (translated below)
  color Silver  // light's color
  // <widthVector> <heightVector> nLightsWide mLightsHigh
  area_light
  <8, 0, 0> <0, 0, 8> // lights spread out across this distance (x * z)
  2, 3                // total number of lights in grid (4x*4z = 16 lights)
  adaptive 0          // 0,1,2,3...
  jitter              // adds random softening of light
  translate <40, 80, -40>   // <x y z> position of light
}

#declare R = 0.9;     // set radius value
#declare L = 1.10;  // set square area side
#declare X = 0.01;  // set co-ordinate x initial value
#declare Y = 0.01;  // set co-ordinate y initial value
#declare Z = 0;  // set co-ordinate z initial value

#declare n = 200;    // set number of pixels per area side
#declare N = 30;   // set number of cycles
#declare Nr = 50;   // set number of sections
#declare Th = -60;
#declare Ph = 20;

#declare P = array[2*n+1][n+1];    //  P matrix
#declare Q = array[2*n+1][n+1];    //  Q matrix

union{ #for (p, -n, n, 1)
  #declare IncX = p*L/n;
  #for (q, 0, n, 1)
      #declare IncY = q*L/n;
      #declare X = IncX;
      #declare Y = IncY;
      #declare P[p+n][q] = 0;
      #declare Q[p+n][q] = 0;

       #for (i,1,N)
       #declare XX = 3*X/4 - X*(X*X - 3*Y*Y)/(X*X + Y*Y)/(X*X + Y*Y)/
       (X*X + Y*Y)/4;
       #declare YY = 3*Y/4 - Y*(Y*Y - 3*X*X)/(X*X + Y*Y)/(X*X + Y*Y)/
       (X*X + Y*Y)/4;
       #declare X = XX;
       #declare Y = YY;
       #if (X*X+Y*Y > R)
       	#if (X*X+Y*Y < R+0.02)
       	 #declare P[p+n][q] = p+n;
       	#declare Q[p+n][q] = q;
#end  // end if
#end  // end if
#end  // end for i
#end  // end for q
#end  // end for p

#declare Rnd_1 = seed (1153);
#declare Rnd_2 = seed (553);
#declare Rnd_3 = seed (758);

#for (j,0,n*n*Nr)      // partial values n*n/64, n*n/16, n*n/4
#declare p = int(2*n*rand(Rnd_1));
#declare q = int(n*rand(Rnd_2));
#declare k = int(Nr*rand(Rnd_3));
 	 	sphere {
  		  < P[p][q], Q[p][q]*cos(2*pi*k/Nr), Q[p][q]*sin(2*pi*k/Nr) >, 1
  		  hollow
		  radiosity { importance 1.0 }
   			 texture {
    			  pigment { color Col_Glass_Bluish }
      		 }
      		 finish { ambient rgb <0.3,0.1,0.1>
     			  diffuse .3
     			 reflection .3
     			 specular 1 }
   		}

#end  // end for j
rotate < 0, Th, Ph >
translate <-150,-50,0> }