Desmos, 401 375 357 353 bytes

-26 bytes thanks to @aiden-chow

-18 bytes thanks to flexabrotnt#1409

H(u,v)=(u.xv.x-v.yu.y,u.xv.y+v.xu.y)
F(u,v)=(u.xv.x+u.yv.y,u.yv.x-u.xv.y)/(v.x²+v.y²)
E(z)=e^{z.x}(cos(z.y),sin(z.y))
C(z,k)=H((1/2,0),E(z)+kE(H((-1,0),z)))
N(z)=F(H(z,C(z,-1))-C(z,1)+(1,0),C(z,-1))
G(z)=max(z.x²,z.y²)
Y(z)=N(N(N(N(N(z)))))
Z(z)=Y(Y(Y(Y(z))))
l=[0...2]
O=hsv(l110,1,1)
P=Z((x,y))
X=round(P.y/(2π))
G((0,2Xπ)-P)<0.01\{\mod(X,3)=l\}

Try It On Desmos!

Before golfing

Rust + Num ASCII Art - 420 413 bytes

-6 bytes thanks to @LF

use{num::complex::*,std::f32::consts::*};(0..255).any(|y|(0..512).any(|x|if let Some(n)=(0..99).scan(Complex::new(x as f32/256.0-1.,y as f32/256.0-0.5+PI),|s,l|{*s=*s-(s.cosh()-1.)/s.sinh();Some((*s,l))}).filter(|s|s.0.re.abs().max(s.0.im-((s.0.im/(2.*PI)).round()*2.*PI).abs())<0.01).next(){print!("\x1b[{:0<2};{:0<2}m0",30+n.1%8,40+((n.0.im/PI/2. + 36.).round()%8.)as usize)}else{print!("0")}>())||print!("
")>());

Python + Pygame, 410 392 bytes (Example Submission)

-16 bytes thanks to wizzwizz4 and others

from pygame import*
from cmath import*
w=256;s=Surface((w,w));f=lambda p:p-(cosh(p)-1)/sinh(p);g=lambda t:round(t.imag/pi/2)*pi*2j;h=lambda p:max(abs(p.real),abs(p.imag))>0.01
for x in range(w):
 for y in range(w):
  t=x/w-0.5+(y/w-0.5-pi)*1j
  for i in range(99):
   try:t=f(t);1/h(g(t)-t)
   except:break
  s.set_at((x,y),[0,[w-1,5,5],[9,w-1,9]][(round(t.imag/pi/2)+i)%3])
image.save(s,"f")

Result will be saved as a TGA file named "f":

JavaScript (ES6), 271 260 253 251 bytes

c=>{z=c.getContext`2d`;for(i=384;i--;)for(j=384;i--;)with(Math){x=i/256-.75;y=j/256-.75+PI;for(k=10;k<30>(x*x<1e-4&abs(y-2*PI*(v=round(y/2/PI)))<.01);k++)w=cosh(x)+cos(y),x-=sinh(x)/w,y-=sin(y)/w;z.fillStyle=`#0${(k+v%10)*3%10}9`;z.fillRect(i,j,1,1)}}

Takes as parameter a 384×384 canvas on which to draw the result. Demo, with whitespace for readability and animation:

f = async c => {
  z = c.getContext`2d`;
  for (i = 384; i--; ) {
    for (j = 384; j--; ) with (Math) {
      x = i / 256 - .75;
      y = j / 256 - .75 + PI;
      for (k = 10; k < 30 > (x * x < 1e-4 & abs(y - 2 * PI * (v = round(y / 2 / PI))) < .01); k++)
        w = cosh(x) + cos(y),
        x -= sinh(x) / w,
        y -= sin(y) / w;
      z.fillStyle = `#0${(k+v%10)*3%10}9`;
      z.fillRect(i, j, 1, 1);
    }
    await new Promise(resolve => requestAnimationFrame(resolve));
  }
};
f(c);

<canvas width=384 height=384 id=c>

Edit: Saved 11 bytes thanks to @Ausername. Saved 7 bytes thanks to @Steffan. Saved 2 bytes thanks to @Dingus.

Pyxplot 0.8.4, ^{190 188} 186 bytes

se nu c
su g(x,y){z=x+y*i
fo n=1to20{z=z-coth(z)+csch(z)
if(Re(z)**2<1e-4)*cos(Im(z))>cos(.01){ret ceil(Im(z)/pi-.5);}
}
}
se sa gr 1e3x1e3
se colm hsb c2:c1:1
p[-1:1][2:4]g(x,y):n w col

Conveniently for code golf, most commands in Pyxplot may be abbreviated zealously. By default (empty .pyxplotrc config file), output is written to pyxplot.eps.

If the text and axis labels produced by default must be removed, replace the last line with

se ax xy in
p[-1:1][2:4]g(x,y):n t''w col

at a cost of 15 bytes. Add se s s for a square aspect ratio.

Some mathematical simplifications are used. Firstly, we write \$g\$ as \$g(z) = z-\coth(z)+\mathop{\rm csch}(z)\$. Secondly, we note that \$\max(\lvert\Delta x\rvert, \lvert\Delta y \rvert) < t\$ implies that both \$\lvert\Delta x\rvert < t\$ and \$\lvert\Delta y\rvert < t\$ and treat these two constraints separately:

• \$\lvert\Delta x\rvert = \lvert x-0\rvert = \lvert x\rvert\$, so the first constraint becomes \$x^2<t^2\$, which is shorter in Pyxplot than \$\lvert x\rvert < t\$ (which would be abs(Re(z))<.01).
• \$\lvert\Delta y\rvert = \lvert y-2k\pi\rvert\$ for some integer \$k\$. Observe that the attractors lie at points where the cosine function is maximal, i.e. \$\cos(2k\pi)=1\$. Thus the constraint on \$y\$ may be written as \$\cos(\lvert y-2k\pi\rvert) > \cos(t)\$, which simplifies to \$\cos(y) > \cos(t)\$ by the symmetry and periodicity of the cosine function.

For the plot, an HSB colourmap is used. At each point the hue is determined by the number of iterations (n in the code), automatically rescaled to span the range from 0 and 1 (i.e. the full range of hues). The saturation is determined by the value of \$2k\$ (the return value of g(x,y)) for the point's attractor, again rescaled to span the full range of saturations. The brightness is fixed at maximum.

An iteration limit of 20 is imposed. Red and orange hues become more predominant at the expense of green and blue as this limit is increased.

In theory, it is possible to omit [-1:1][2:4] in the last line of the code (saving 10 bytes), have Pyxplot apply its default axis ranges (\$-10\le x,y\le 10\$), and increase the grid resolution (1e3x1e3) to compensate. In practice, setting the resolution to 4e3x4e3 exhausted the available memory (16 GB RAM plus 20 GB swap) on my machine.

Ungolfed code

set numerics complex
subroutine g(x, y) {
  z = x + y * i
  for n = 1 to 20 {
    z = z - coth(z) + csch(z)
    if Re(z)**2 < 1e-4 and cos(Im(z)) > cos(1e-2) {
      return ceil(Im(z) / pi - 0.5)
    }
  }
}
set samples grid 1e3 x 1e3
set colourmap hsb c2:c1:1
plot [-1:1] [2:4] g(x, y):n with colourmap

Wolfram Language (Mathematica), 128 bytes

Image@Table[n=0;x+y*I//.z_/;Max@Abs@{Re@z,Mod[i=Im@z,Pi,n+=8t;-1]}>2t:>z-Coth@z+Csch@z;i~Hue~n,{y,-1+Pi,1+Pi,t=.005},{x,-1,1,t}]