Thank you for your reply. I appreciate that!
So I can use the the DEfactor or the distance after the mandelbox part to apply it's own modification of a mandelbulb. The thing is I don't know how since there were trigonometric functions. I found a method for the mandelbulb on the web of Iņigo Quilez that could be implemented by using the derivative which can be calculated in the iteration loop. But that gives me finally a distance and I can't apply this to the DEfactor or distance calculated before in the mandelbox part.
Where can I find your Delta DE method? I looked here in the forum and you mention it a lot but I couldn't find it (or I didn't recognise it). It would be nice if you could send me a link or something because I get really lost when you talk about smooth iteration value

a space of 5 iterations touching a space of 7 with no space of 6 in between
I don't understand that at all. I'm using the same number of iterations for the entire fractal to calculate the DE
You don't need to worry about the trig too much, in the case of polynomials (even using the odd Mandelbulb triplex maths) it turns out that the running DE is exactly the same as that for plain complex or quaternions etc.
i.e. If your current running DE value is d and you're applying z^p+c on a given iteration then the new running DE value is
d = p*d*z^(p1)
(That may need a +1 on the end for either Mandy's or Julias I forget which)
Of course the z^(p1) needs to be calculated using the correct number form i.e. using triplex math for the Mandelbulb matching the triplex math used for z^p+c.
This is "off the top of my head" and it's a while since I messed with it and my memory isn't what it was but...
My delta DE is simply done by calculating the smooth iteration value for the current point on the ray, for a point deltax up the x axis, a point deltay (magnitude same as deltax) up the y axis and a point deltaz (same magnitude) up the z axis  the delta being a suitably small value according to the current required distance threshold for "solid".
Of course this is an overhead of 3 full iteration loops per step along the ray compared to using analytical DE *but* of course the analytical DE calculation involves calculations done on every iteration whereas the smooth iteration value only involves calculation on bailout.
You then take the differences in the smooth iteration value i.e. (value for point offset in x  value for base point (at step on ray), value for point offset in y  value for base point, value for point offset in z  value for base point) and get the magnitude of this (as a 3D vector) and then adjust appropriately to give a DE value.
Buddhi's method I think is similar except if I remember correctly he iterates at the delta offsets always to the same iteration count as the base point (at current step on the ray) and uses the bailout values rather than the smooth iteration values.
With respect to the iteration count  I was talking about points along the ray before you reach solid  here the ray will pass through space of different iteration levels at bailout and for many fractals (like the Mandelbulb) although the shape of the boundaries from one iteration depth to the next are fractal and irregular they are always consistent i.e. there will never be a point of 5 iterations at bailout directly next (I mean in the infinitessimal limit) to a point of 7 iterations at bailout  there's always a point between of 6 iterations at bailout.
To a certain degree even analytical DE relies on this but some fractals, notably some fractals with conditionals in the formula and fractals with fractional powers, have adjacent points where the iteration count at bailout is discontinuous and such formulas always present rendering issues due to "incorrect" calculation of the distance estimate.