I'm not sure if Mathematica is the right tool to make these plots, but well we can try. First we can compile and optimize your code a bit:
ClearAll[cf,cfg]
cf=Compile[{{x, _Real},{y, _Real},{z, _Real}},
Module[{o=1,n=8,newx=x,newy=y,newz=z,rsq=x x+y y+z z,theta,phi},
While[o<400&&rsq<=4,
If[newz==0.0,newz=0.000001];
If[newz==0.0,newx=0.000001];
If[newy==0.0,newy=0.000001];
theta=n ArcTan[Sqrt[newx*newx+newy*newy],newz];
phi=n ArcTan[newy,newx];
rsq=rsq^(n/2);
newx=rsq Sin[theta]Cos[phi]+x;
newy=rsq Sin[theta]Sin[phi]+y;
newz=rsq Cos[theta]+z;
rsq=newx newx+newy newy+newz newz;
o=o+1
];
o
]
]
cfg[x_?NumberQ,y_?NumberQ,z_?NumberQ]:=cf[x,y,z]
Compared to your uncompiled version it is 11x faster:
RepeatedTiming[Mandelbulb[0.1, 0.1, 0.1], 2]
RepeatedTiming[cf[0.1, 0.1, 0.1], 2]
RepeatedTiming[cfg[0.1, 0.1, 0.1], 2]
0.0056
0.00052
0.00053
on my machine...
n = 45;
pts = Subdivide[-1.5, 1.5, n];
AbsoluteTiming[
vals = Table[cf[x, y, z], {x, pts}, {y, pts}, {z, pts}];]
Times @@ Dimensions[vals]
ListContourPlot3D[vals, Contours -> {20}, Mesh -> None]
Now we use ContourPlot rather than DensityPlot to get a surface:
With already a much better result, and only in 4 seconds (my laptop). Again, I'm not sure if Mathematica is the right tool, but with some optimisations it might be possible.
