1. There's no docstring for fd_solve. The post contains a description that would make a good start.

2. The docstring has zs[2:, 1:-1, 1:-1] = 0 but this doesn't match the code zs[1:,1:-1,1:-1] = 0.

3. eqn_parse does not seem to parse anything. I think you could pick a better name, for example pde_kernel.

4. eqn_parse might return None in some circumstances. Is this the expected behaviour? It would be better to raise an exception if the input is incorrect.

5. When building a dictionary, it is often clearer to call dict instead of using a dictionary literal. When calling dict you don't have to quote the keys, which reduces the amount of syntactic clutter. For example:

    sliders_dict = dict(active=0,
                        currentvalue=dict(font=dict(size=20),
                                          prefix='t: ',
                                          xanchor='right'),
                        steps=steps,
                        transition=dict(duration=100))
   

6. Instead of:

    ts[:,None,None], xs[None,:,None], ys[None,None,:]
   

I suggest using numpy.meshgrid to make it clear that you are constructing a grid of coordinates:

    np.meshgrid(ts, xs, ys, sparse=True, indexing='ij')

1. Instead of (t-t)+1, I suggest using np.ones_like(t) to make it clear that you are constructing an array of ones with the same shape and dtype as t. (See numpy.ones_like.) But since this is constant, why not just write zs[0,1:-1,1:-1] = 1 and avoid the need for ic?

2. The slice 1:-1 appears in many places; I suggest making a global constant, for example INTERIOR = slice(1, -1).

3. The boundary condition computation is wasteful as most of the values computed will be overwritten in the subsequent steps. It might be better to start with

    zs = np.zeros(ts.shape + xs.shape + ys.shape)
   

and then fill in the values only in the boundary planes:

    zs[0, INTERIOR, INTERIOR] = 1
    zs[:, (0, -1), :] = bc(*np.meshgrid(ts, (xs[0], xs[-1]), ys, sparse=True, indexing='ij'))
    zs[:, INTERIOR, (0, -1)] = bc(*np.meshgrid(ts, xs[INTERIOR], (ys[0], ys[-1]), sparse=True, indexing='ij'))