Solve the Laplace Equation $u_{xx} + u_{yy} = 0$ inside the square $0 \leq x \leq 1$, $0 \leq y \leq 1$ 
Solve the Laplace Equation $u_{xx} + u_{yy} = 0$ inside the square $0 \leq x \leq 1$, $0 \leq y \leq 1$ subject to the boundary conditions $u_x(0, y) = 0$ and $u_x(1, y) = 0$ for all $y \in [0, 1]$ and with initial condition $u(x, 0) = 0$ and $u(x, 1) = 1-x$

So far what I've done is the following. I've assumed that the solution is given by $u(x, y) = X(x)Y(y)$, then seperating variables and applying $u_{xx} + u_{yy} = 0$ I get $$\frac{X''(x)}{X(x)} = \frac{Y''(y)}{Y(y)} = - \lambda^2$$ 
This gives me two ODE's $X''(x) = \lambda^2 X(x)$ and  $Y''(y) = \lambda^2 Y(y)$ and hence I get solutions to these ODE's as $X(x) = A \sin \lambda x + B \cos \lambda x$ and $Y(y) = C_1e^{\lambda y} +  C_2e^{-\lambda y}$.
Applying boundary conditions I get $A = 0$ and $\lambda = n\pi$ and $X(x) = \cos (n\pi x)$
Applying initial conditions (specifically $u(x, 0) = 0$) I get a system of two equations $$C_1 = \frac{-(1-x)e^y}{1-e^{2y}}$$ and $$C_2 = (1-x)e^y+ \frac{(1-x)e^3y}{1-e^2y}$$ 
Is what I've done so far correct? If so how can I proceed and solve the Laplace equation?
 A: Your answer is correct up until the $x$ part. Note that there are infinitely many values of $\lambda$ that satisfy the equation, so the product solution
$$u_n(x,y) = \big[c_1e^{n\pi y} + c_2e^{-n\pi y}\big]\cos(n\pi x) $$
Is just one of many. Now since all of these solutions are linearly independent, we can employ the law of superposition and write the full solution as
$$ u(x,y) = \sum_{n=0}^\infty k_n u_n(x,y) = A_0 + B_0y + \sum_{n=1}^\infty \big[A_n e^{n\pi y} + B_n e^{-n\pi y}\big]\cos(n\pi x) $$
Where I've renamed the constants $A_n \mapsto k_nc_1$ and $B_n \mapsto k_nc_2$. Also note the special case $\lambda = 0$ where the $y$ part has a slightly different form.
Applying the remain B.C.s gives
\begin{align}
u(x,0) &= A_0 + \sum_{n=1}^\infty (A_n + B_n)\cos(n\pi x) = 0 \\
u(x,1) &= A_0 + B_0 +\sum_{n=1}^\infty (A_ne^{n\pi} + B_ne^{-n\pi})\cos(n\pi x) = 1-x
\end{align}
You can apply the Fourier series expansion to find equations determining $A_n,B_n$
\begin{align} 
A_0 &= 0, & B_0 &= \int_0^1 (1-x)\ dx \\
A_n+B_n &= 0, & A_ne^{n\pi} + B_ne^{-n\pi} &= 2\int_0^1(1-x)\cos(n\pi x)\ dx
\end{align}

Edit: The second equation is really the one you care about (the first should be trivial). The basic idea is to rewrite the boundary function as a Fourier series
$$ f(x) = 1-x = f_0 + \sum_{n=1}^\infty f_n \cos(n\pi x) $$
where $f_i$ are all constants. Values $n<0$ are excluded due to symmetry. 
We say that the functions $\{\cos(n\pi x)\}$ from a complete orthogonal basis in the space $x\in [0,1]$. 
Orthogonal means the inner product
$$ \int_0^1 \cos(n\pi x)\cos(m\pi x) dx $$
Is equal to zero when $n \ne m$, and a non-zero constant when $n=m$. Specifically
$$ \int_0^1 \cos^2 (n\pi x) dx = \begin{cases}
1, & n=1\\
\frac12, & n \ne 1
\end{cases} $$
Complete means every continuous, bounded function in the interval $[0,1]$ can be expressed as a series of this form. The series converges to an even, periodic extension of $f(x)$
Using these facts, we can integrate across the sum and prove that
$$  \int_0^1 f(x) \cos(m\pi x)\ dx = \sum_{n=0}^\infty f_n \int_0^1 \cos (n\pi x)\cos(m\pi x) \ dx $$
I've skipped over the justification for interchanging integration and summation, but the key point here is: every term in this series will vanish except when $n=m$, leaving
$$  \int_0^1 f(x) \cos(m\pi x)\ dx = f_m\int_0^1 \cos^2 (m\pi x)\ dx $$
For a more general formulation including sine functions and an arbitrary interval, check out Wikipedia's Fourier Series page.
A: To find the separated solutions $u(x,y)=X(x)Y(y)$, substitute into the equation and impose the homogeneous boundary conditions in order to obtain
$$
                X''(x)Y(y) + X(x)Y''(y) = 0, \\
                -\frac{X''}{X} = \lambda = \frac{Y''}{Y}, \\
                X'(0)=0=X'(1),\;\; Y(0)=0.
$$
To normalize the solutions $Y$, further assume $Y'(0)=1$, which gives
$$
         Y(y)=\frac{\sinh(\sqrt{\lambda}y)}{\sqrt{\lambda}}.
$$
This works for $\lambda=0$ in the limit as $\lambda\rightarrow 0$, where $Y(y)=y$. The non-trivial solutions $X$ are
$$
            X_n(x) = \cos(n\pi x),\;\; n=0,1,2,3,\cdots,
$$
with corresponding $\lambda_n = n^2\pi^2$. The full solution $u$ has the form
$$
         u(x,y)= A_0 y + \sum_{n=1}^{\infty}A_n\cos(n\pi x)\sinh(n\pi y)
$$
The constants $A_n$ must be chosen so that $u(x,1)=1-x$, leading to
$$
         1-x = A_0 + \sum_{n=1}^{\infty}A_n \cos(n\pi x)\sinh(n\pi).
$$
Now use the mutual orthogonality of the functions $\cos(n\pi x)$ to isolate the constants $A_n$. To do this, multiply both sides by $\cos(k\pi x)$ and integrate over $[0,\pi]$ to obtain the equations for $A_k$:
$$
       \int_{0}^{1}(1-x)dx = A_0\int_{0}^{1} dx \\
     \int_{0}^{1}(1-x)\cos(k\pi x)dx = A_k\int_{0}^{1}\cos^2(k\pi x)dx\sinh(k\pi).
$$
