Thank you. I will try this method and check whether it is faster than mine. The 3 things I found and work well are
reduce the resolution of my matrices (for example 41x41x41 works quite well) because the interpolating functions will smooth everything anyway;
plot less points in my graph (with a smaller range of x,y,z) and use the command 'PlotPoints' in the following way
ContourPlot3D[{ReI[x,y,z] == 0, ImI[x,y,z]== 0},
{x, 1, 41}, {y, 1, 41}, {z, 1,41}, ContourStyle -> Opacity[0.1],
Mesh -> None,BoundaryStyle -> {{1,2} -> {{Green,Tube[0.5]}}}, AxesLabel->{x,y,z},PlotPoints -> 8]
- decrease the precision (in this case I use the mesh function method)
ContourPlot3D[{ReI[x,y,z] == 0, ImI[x,y,z]== 0},
{x, 1, 41}, {y, 1, 41}, {z, 1,41}, MeshFunctions -> {Function[{x, y, z}, ReI[x,y,z] - ImI[x,y,z]]},
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}}, ContourStyle -> Opacity[0.1],
AxesLabel->{x,y,z},
WorkingPrecision -> 2]