Acá defino el potencial correspondiente a una caja de largo, ancho y alto 1, con todas sus tapas a potencial cero excepto la tapa de arriba que está a potencial 1. Elijo un sistema de coordenadas en la esquina de la caja. Dado que mi solución es una doble serie, tengo un parámetro N, que me controla el límite superior de la suma. Dado que son dos sumas y ambas arrancan desde cero, para un dado N, la siguiente definicion suma (N+1)x(N+1) términos.

f[x_,y_,z_,N_]:=NSum[(16/(π^2))((Sin[(2n+1)π  x]Sin[(2m+1)π  y])/((2n+1)(2m+1)))((Sinh[π (Sqrt[(((2n+1))^2)+(((2m+1))^2)])z])/(Sinh[π (Sqrt[(((2n+1))^2)+(((2m+1))^2)])])),{n,0,N},{m,0,N}]

Como se puede ver en los siguientes dos ejemplos, con N entre 4 y 5 sobra para tener una precisión de 7 dígitos. Observar que los valores x, y, z los elijo entre 0 y 1. Esto se debe a que como mi caja tiene lado 1, el interior de la misma son valores entre 0 y 1. Observar que a medida que me acerco a alguna de las paredes, el N que necesito para obtener precisión aumenta.

Manipulate[f[0.2,0.3,0.5,N],{N,0,20,1}]

Manipulate[f[0.8,0.33,0.18,N],{N,0,20,1}]

Acá ploteo para valores de y, z fijos, moviendome en x. Como es de esperar el potencial va de 0 a 0. La línea más alta y la más baja son las más cercanas a la pared y por lo tanto las que requieren más precisión.

Plot[{f[x,0.1,0.1,4],f[x,0.3,0.3,4],f[x,0.5,0.5,4],f[x,0.7,0.7,4],f[x,0.9,0.9,6]},{x,0,1}]

Acá ploteo para valores de x e y fijos, moviendome en z. Como es de esperar el potencial va de 0 (tierra) a 1 (potencial de referencia), pero nuevamente cerca de las paredes hay cierto error y debería considerar N más grande (costoso computacionalmente).

Plot[{f[0.5,0.5,z,4],f[0.6,0.6,z,5],f[0.7,0.7,z,5],f[0.8,0.8,z,6],f[0.9,0.9,z,6]},{z,0,1}]