So the semi-infinite strip means that we are working in the region where $(x,y)\in [0,\infty)\times [0,H]$. We have Laplace's equation, which simply states:
$$\nabla^{2} \Phi = 0 \iff \frac{\partial^{2} \Phi}{\partial x^{2}}+\frac{\partial^{2} \Phi}{\partial y^{2}}=0$$
And we have $\Phi(0,y)=f(y)$ and $\lim_{x\to\infty}\Phi(x,y) = 0$ (from physical conditions) as our boundary conditions, along with:
$$\left.\frac{\partial \Phi}{\partial y}\right|_{(x,y)=(x,0)} = \left.\frac{\partial \Phi}{\partial y}\right|_{(x,y)=(x,H)}=0$$
We note that using the method of characteristics we can assume that: $\Phi(x,y)=X(x)Y(y)$, and therefore we can rewrite the Laplace equation above:
$$Y(y)\frac{\partial^{2} X}{\partial x^{2}}+X(x)\frac{\partial^{2} Y}{\partial y^{2}}=0 \implies \frac{1}{X(x)}\frac{\partial^{2} X}{\partial x^{2}}=-\frac{1}{Y(y)}\frac{\partial^{2} Y}{\partial y^{2}}$$
But $x$ and $y$ are independent, and therefore we must have that:
$$\frac{1}{X(x)}\frac{\partial^{2} X}{\partial x^{2}}=k^{2} \land \frac{1}{Y(y)}\frac{\partial^{2} Y}{\partial y^{2}}=-k^{2}$$
These are ODEs with known solutions:
$$X(x)=\alpha \exp(kx)+\beta\exp(-kx) \land Y(y) = \gamma \sin(ky) + \delta \cos(ky)$$
Where $\exists \alpha,\beta,\gamma,\delta\in\mathbb{R}$. We now wish to find the value of these coefficients based on the initial conditions. We can see that our second boundary condition, $\lim_{x\to\infty}\Phi(x,y)=0$ means that $X(x)$ must have a decaying form, thus $\alpha = 0$. We can now use the differential boundary conditions:
$$\left.\frac{\partial \Phi}{\partial y}\right|_{(x,y)=(x,0)}=0 \implies \gamma k \cos(0) - \delta k \sin(0) = 0$$
Therefore we have $\gamma k = 0$, and therefore $\gamma = 0$ (if we assume $k \neq 0$ for non-static solution). We can now look at the second differential boundary condition:
$$\left.\frac{\partial \Phi}{\partial y}\right|_{(x,y)=(x,H)}=0 \implies \delta k \sin(k H)=0$$
But assuming $\delta \neq 0$ (for non-trivial solution), we must have $k = \frac{n \pi}{H}$, $\exists n \in \mathbb{N}_{0}$.
We finally turn to our last boundary condition $\Phi(0,y)=f(y)$, where $f(y)$ is an arbitrary function. We can write this as:
$$\sum_{n = 0}^{\infty} \delta_{n} \beta_{n} \sin(k_{n} y) = f(y)$$
But WLOG we can assume $\beta_{n} = 1$, $\forall n \in \mathbb{N}_{0}$ therefore we have:
$$\sum_{n=0}^{\infty}\delta_{n} \sin(k_{n} y) = f(y)$$
We can recognize this as a Fourier sine series, thus giving us coefficients:
$$\delta_{n} = \frac{2}{H}\int_{0}^{H} f(y) \sin\left(\frac{n \pi y}{H}\right)\:\mathrm{d}y\tag{1}$$
So we have in general:
$$\Phi(x,y) = \sum_{n=0}^{\infty}\delta_{n}\exp\left(-\frac{n\pi x}{H}\right)\sin\left(\frac{n \pi y}{H}\right)$$
Where $\delta_{n}$ is defined by $(1)$.
I hope this helps clear things up!
This is because you are considering the equation on the half-line. Your domain is essentially $\Omega=[0,\infty)\times [0,y]$
If you look at the heat equation and wave equation on a bounded domain, you will see the solutions are typically given as Fourier Series, but on unbounded domains, the solutions look a lot different. You cannot use boundary methods and separation of variables to solve the PDE here.
Generally speaking, we want solutions $u(x,y)$ which "vanish at infinity" (tend towards $0$ as $||u(x,y)||_{L^2}\rightarrow \infty)$. Functions which are smooth and have compact support (denoted by $C_{c}^{\infty}(\Omega)$) include those that vanish at infinity $C_{0}(\Omega)$. These conditions imposed on the solutions are desirable because they make physical sense. The temperature/potential in a rod is not going to shoot off to infinity, it's going to decrease with time until it reaches a steady-state.
Now to actually solve the PDE, I don't think there is a slick and easy way to do it. Solutions to Laplace's equation in domains which aren't on a disk $D^n$ tend to be ad-hoc and very specific to the domain. Polynomials wouldn't work because they do not vanish at infinity. I would try to look for solutions which use sines and cosines. Alternatively, you could look for a harmonic function on the domain $[0,\infty)\times [0,y]$ (you don't necessarily need to work in the complex plane for this). Start by looking for qualitative solutions and then try to get a concrete, quantitative one.
Best Answer
You want homogenous conditions $u(0,y)=0=u(4,y)$ instead, which can achieved by subtracting $T$ from the $u$ in order to obtain a new equation for $v=u-T$ $$ v_{xx}+v_{yy} = 0, \;\; v(0,y)=0=v(4,y),\;\; v(x,0)=-T. $$ Now when you separate variables $v(x,y)=X(x)Y(y)$, you obtain $$ -\frac{X''}{X} = \lambda = \frac{Y''}{Y}. $$ The negative can go in either place, but I know how it works out, so I chose the above so that $\lambda$ is positive instead of negative. Then $X$ must satisfy $$ X''+\lambda X = 0,\;\; X(0)=X(4)=0. $$ The solutions are $$ \lambda_n=\frac{n^2\pi^2}{16},\;\; X_n(x)=\sin\left(\frac{n\pi}{4}x\right) $$ The corresponding solutions in $Y$ are $$ Y_n = A_n\exp\left(-\frac{n\pi}{4}y\right). $$ (I've excluded solutions with positive exponent in order to have a bounded solution.) The constants $A_n$ must be chosen so that $$ v(x,y) = \sum_{n=1}^{\infty}A_n\sin\left(\frac{n\pi}{4}x\right)\exp\left(-\frac{n\pi}{4}y\right),\;\;\; v(x,0)=-T. \\ \implies \sum_{n=1}^{\infty}A_n\sin\left(\frac{n\pi}{4}x\right)=-T \\ $$ Multiplying both sides by $\sin\left(\frac{n\pi}{4}x\right)$, integrating in $x$ over $[0,4]$, and using the orthogonality of the eigenfunctions determines the constants $A_n$ by the equations $$ A_n\int_{0}^{4}\sin^{2}\left(\frac{n\pi}{4}x\right)dx=-T\int_{0}^{4}\sin\left(\frac{n\pi}{4}x\right)dx. $$ The solution of the original problem is $$ u(x,y) = T+v(x,y) = T+\sum_{n=1}^{\infty}A_n\sin\left(\frac{n\pi}{4}x\right)\exp\left(-\frac{n\pi}{4}y\right). $$