The problem is :
Let $$HS = \{(x,y):0\leq x\leq 1, y\geq 0\}$$ and consider the steady state heat problem on $HS$ with boundary conditions $u(0,y)=u(1,y)=0$ and $u(x,0) = x-x^2$.
Use separation of variables to obtain a family of basic solutions.
so I let $$u(x,y) = X(x)Y(y)$$
Since $u$ is the solution of heat equation,
$$u_{xx} + u_{yy} = X''Y + XY'' = 0$$
Manipulate this..
$$\frac{X''}{X} = -\frac{Y''}{Y}$$
But I don't know how to continue
Assume a solution of the form $$ u(x,y)=X(x)Y(y) $$ and as you did, find $$ \frac{X''(x)}{X(x)}=-\frac{Y''(y)}{Y} $$ which then implies both quantities are some constant, say $\lambda$. Solve the two second order linear ODE's $$ Y''(y)+\lambda Y=0\\ X''(x)-\lambda X=0\\ \implies Y(y)=ae^{\sqrt{-\lambda}y}+ce^{-\sqrt{-\lambda}y}\\ X(x)=b_1e^{\sqrt{\lambda}x}+b_2e^{-\sqrt{\lambda}x} $$ making your solution the product of the two solutions, $$ u(x,y)=X(x)Y(y)=(b_1e^{\sqrt{\lambda}x}+b_2e^{-\sqrt{\lambda}x})(ae^{\sqrt{-\lambda}y}+ce^{-\sqrt{-\lambda}y}) $$ Now use your boundary condition to figure out all the constants, $$ u(0,y)=(b_1+b_2)(ae^{\sqrt{-\lambda}y}+ce^{-\sqrt{-\lambda}y})=0\\ u(1,y)=(b_1e^{\sqrt{\lambda}}+b_2e^{-\sqrt{\lambda}})(ae^{\sqrt{-\lambda}y}+ce^{-\sqrt{-\lambda}y})=0 $$ in order for these to vanish identically for all $y$, we need $$ b_1=-b_2\\ b_1e^{2\sqrt{\lambda}}=-b_2\implies e^{2\sqrt{\lambda}}=1 $$ unless $b_1=0=b_2$ and we have the zero solution.
We could pick $\lambda=0$, but this would again give us the zero solution. So we take $$ 2\sqrt{\lambda}=2\pi i k\implies \lambda=-\pi^2k^2 $$ for any integer valued $k$.
Where has this left us, now we have $$ u(x,y)=X(x)Y(y)=(b_1e^{\pi i kx}-b_1e^{-\pi i kx})(a_1e^{\pi |k|y}+a_2e^{-\pi|k|y})=2b_1i\sin(\pi kx)(ae^{\pi |k|y}+ce^{-\pi|k|y}) $$ for any integer $k$. Then using superposition, we may sum over all integer $k$ to build the solution, with coefficients depending on $k$, $$ u(x,y)=\sum_{k\in \mathbb{Z}}B_k\sin(\pi kx)(a_ke^{\pi |k|y}+c_ke^{-\pi|k|y}) $$
Let us now impose your final condition, $$ u(x,0)=x-x^2=\sum_{k\in \mathbb{Z}}B_k\sin(\pi kx)(a_k+c_k) $$ where we can integrate away our troubles. Using orthogonality of trig functions, we integrate against $\sin(nx)$ over $[-1,1]$, with $\frac{1}{2}$ a normalizing constant. This yields, $$ \frac{1}{2}\int_{-1}^1(x-x^2)\sin(\pi n x)\mathrm dx=B_n(a_n+c_n)\\ \stackrel{\text{symmetry}}{\implies} \int_{0}^1x\sin(\pi n x)\mathrm dx=\frac{-1}{\pi n}\cos(\pi n)=B_n(a_n+c_n) $$ Here I think we are stuck unless we impose a condition at $y=\infty$, namely to insure that $$ \lim_{y\to \infty}|u(x,y)| $$ exists. This requires $a_k=0$ for any $k\in \mathbb{Z}$. Yielding for our coefficient $$ \frac{-1}{\pi n}\cos(\pi n)=B_nc_n $$ and our solution $$ u(x,y)=\sum_{k\in \mathbb{Z}}\frac{-1}{\pi k}\cos(\pi k)\sin(\pi kx)e^{-\pi|k|y} $$
We have $X''(x)/X(x) = -Y''(y)/Y(y)$ for all $x \in [0,1]$ and all $y \ge 0$.
Then there is a constant $c$ such that $X''(x)/X(x) =c$ for all $x \in [0,1]$ and $Y''(y)/Y(y)=-c$ for all $y \ge 0$.
Hence we get the differential equations $X''=cX$ and $Y''=-cY$.
Can you take it from here ?
