This is very involved problem. To solve, you'll need to break the solution into successively smaller pieces. First, separate the steady-state and transient solutions, then split up the boundary conditions in order to use separation of variables.
To start off, I'm going to label
$$ T(x,y,t) = u(x,y,t) + v(x,y) $$
where $v(x,y)$ is the time-independent, steady-state solution, and $u(x,y,t)$ is the decaying, transient solution.
The steady-state solution should satisfy $\nabla^2v = 0$ and all the boundary conditions as listed. Since all boundaries are inhomogeneous, we need so split it up further
$$ v(x,y) = v_1(x,y) + v_2(x,y) $$
such that
\begin{matrix}
\begin{cases}
\nabla^2 v_1 = 0 \\ \\
v_1(0,y) = T_1,\ v_1(L,y) = T_0 \\ \\
\dfrac{\partial}{\partial y}v_1(x,0) = \dfrac{\partial}{\partial y}v_1(x,L) = 0
\end{cases}
&&&
\begin{cases}
\nabla^2v_2 = 0 \\ \\
v_2(0,y) = v_2(L,y) = 0 \\ \\
\dfrac{\partial}{\partial y}v_2(x,0) = \dfrac{\partial}{\partial y}v_2(x,L) = a
\end{cases}
\end{matrix}
The homogeneous boundaries allow us to use separation of variables to solve each individual problem
- The first one is easier since you can intuitively guess that it's constant in $y$ and linear in $x$. This turns out to be
$$ v_1(x,y) = T_0\frac{x}{L} + T_1\frac{L-x}{L} = \frac{(T_0-T_1)x}{L} + T_1 $$
- For the second problem, the homogeneous boundary condition on $x$ returns a series solution of the form
$$ v_2(x,y) = \sum_{n=1}^\infty \sin\left(\frac{n\pi x}{L}\right)\left[A_n\cosh\left(\frac{n\pi y}{L}\right) + B_n \cosh\left(\frac{n\pi(L-y)}{L}\right)\right] $$
Then, applying the remaining boundary conditions on $y$ will give
$$ A_n = -B_n = \begin{cases} \dfrac{4aL}{n^2\pi^2\sinh(n\pi)}, & n \text{ odd} \\ 0, & n \text{ even} \end{cases} $$
For a more thorough explanation of why I'm using hyperbolic functions instead of exponentials, check out my answer for this similar problem. Short answer: It makes the math easier
The remaining transient solution is homogeneous on all boundaries and has initial conditions that cancel out the steady-state:
\begin{cases}
\dfrac{\partial u}{\partial t} = \nabla^2 u \\ \\
u(0,y,t) = u(L,y,t) = \dfrac{\partial}{\partial y}u(x,0,t) = \dfrac{\partial}{\partial y}u(x,L,t) = 0 \\ \\
u(x,y,0) = T_0 - v(x,y)
\end{cases}
Applying separation of variables once more and matching the homogeneous boundaries, we obtain
$$ u(x,y,t) = \sum_{n,m} c_{n,m} \sin\left(\frac{n\pi x}{L}\right)\cos\left(\frac{m\pi y}{L}\right)\exp \left[-(n^2+m^2)\frac{\pi^2}{L^2}t\right] $$
The initial condition for this last piece will have you solve a double Fourier series in $x$, $y$. You can use linearity to simplify calculations a bit, i.e.
\begin{align}
u_1(x,y,0) &= T_0 - T_1 - \frac{T_0-T_1}{L}x = \sum_{n,m} c_{1(n,m)} \sin\left(\frac{n\pi x}{L}\right)\cos\left(\frac{m\pi y}{L}\right) \\
u_2(x,y,0) &= -\sum_{n=2k+1} \frac{4aL}{n^2\pi^2\cosh(n\pi)} \sin\left(\frac{n\pi x}{L}\right)\left[\cosh\left(\frac{n\pi y}{L}\right) - \cosh\left(\frac{n\pi(L-y)}{L}\right)\right] \\ &\qquad = \sum_{n,m} c_{2(n,m)} \sin\left(\frac{n\pi x}{L}\right)\cos\left(\frac{m\pi y}{L}\right)
\end{align}
The first boundary function is constant in $y$ and the second is already a partial Fourier series, so they simplify to
\begin{align}
T_0 - T_1 - \frac{T_0-T_1}{L}x &= \sum_{n=1}^\infty c_{1(n,0)} \sin\left(\frac{n\pi x}{L}\right) \\
-\frac{4aL}{n^2\pi^2\sinh(n\pi)} \left[\cosh\left(\frac{n\pi y}{L}\right) - \cosh\left(\frac{n\pi(L-y)}{L}\right)\right] &= c_{2(n,0)} + \sum_{m=1}^\infty c_{2(n,m)}\cos\left(\frac{m\pi y}{L}\right)
\end{align}
The last round of integration finishes the job
Best Answer
You are correct to use separation of variables to break up the Partial Differential Equation (PDE) into two independent Ordinary Differential Equations (ODEs). The next step is to solve these equations, then we will find the constants.
First, you should use $ \frac {X''}{X} = \frac{4T'}{T} = - \lambda^2$ as this gives you a negative only eigenvalue. We want negative because we cannot match the boundary conditions with a zero or positive eigenvalue (check the case of $0$ to convince yourself of this, the proof for $\lambda^2$ shown below). You do in general need to check these possibilities, but for this PDE with these boundary conditions the eigenvalues are strictly negative. This will give us a non-zero imaginary number in the ODE solving for $X$.
The independent ODEs are
$ X'' = - \lambda^2 X$
and
$T'= - \frac 14\lambda^2T$
This will give you the following,
$X(x) = A\sin( \lambda x)+B\cos( \lambda x)$
$T(t) = Ce^{\frac {-\lambda^2 t} 4}$
After which you plug in your boundaries of $U(0,t)=0$, $U(2,t)=0$ giving you
$X(0) = A\sin(0)+B\cos(0) = B =0$
And $X(2) = A\sin( 2\lambda ) = 0$ (already removed $B$)
From the second equation we can see that $\lambda$ must be a multiple of $\frac \pi 2$, in fact any integer multiple of $\frac \pi 2$ must be considered. You have a countably infinite number of solutions since every single integer multiple is a valid answer for the question, we will denote these solutions $X_n(x)=C_n \sin(n\frac\pi2x)$ with $\lambda_n =n\frac\pi2$. Since this is a linear differential equation we know the linear combination of solutions is also a solution so we must sum up all possible answers.
By summing all possible solutions together, and multiplying the $T(t)$ term back in we get
$U(x,t) = \sum_{n=1}^\infty T_n(t)X_n(x)$ (remember, $T(t)$ is also a function of $\lambda_n$) Which we simplify to
$U(x,t) = \sum_{n=1}^\infty C_ne^{\frac {-\pi^2 n^2t} {16}}\sin(n\frac\pi2x)$
This is the general form of the answer, all that is left is determining the constants.
In order to do this, we will use Fourier's trick (see link for more details) which uses the initial condition $U(x,0)=2\sin(\frac\pi2x)-\sin(\pi x) + 4\sin(3\pi x)$ and boundary $ 0 \le x \le 2$ to find our undetermined constants.
First we find $U(x,0)$ given our formula as $U(x,0) = \sum_{n=1}^\infty C_n \sin(n\frac\pi2x)$ because $e^0=1$
we then set this equal to $2\sin(\frac\pi2x)-\sin(\pi x) + 4\sin(3\pi x)$ on the boundary $ 0 \le x \le 2$ which we know from the boundary conditions.
$U(x,0) = \sum_{n=1}^\infty C_n \sin(n\frac\pi2x) = 2\sin(\frac\pi2x)-\sin(\pi x) + 4\sin(3\pi x)$
This is where the trick comes in. We will first multiply both sides by $\sin (n'\frac\pi2x)$ where $n'$ is an integer. After that, we integrate in $x$ over the boundary $ 0 \le x \le 2$
This gives us \begin{align} &\sum_{n=1}^\infty C_n \int_0^2 \sin( n\frac\pi2x)\sin (n'\frac\pi2x) dx = \\ &\qquad=\int_0^2 2\sin(\frac\pi2x)\sin (n'\frac\pi2x)-\sin(\pi x)\sin (n'\frac\pi2x) + 4\sin(3\pi x)\sin (n'\frac\pi2x) dx \end{align} (notice we pulled the summation out of the integral)
We will take advantage of the following property.
if $n = n'$ then $\int_0^2 \sin( n\frac\pi2x)\sin (n'\frac\pi2x) dx = \frac 2 2 = 1$
if $n \neq n'$ then $\int_0^2 \sin(n\frac\pi2x)\sin (n'\frac\pi2x) dx = 0$
in essence $\int_0^2 \sin(n\frac\pi2x)\sin (n'\frac\pi2x) dx = \begin{cases} 1, & n = n' \\ 0, & n \neq n'\end{cases}$
Therefore \begin{align} &\sum_{n=1}^\infty C_n \int_0^2 \sin( n\pi x )\sin (n' \pi x) dx =\\ &\qquad=C_{n} = \int_0^2 2\sin(\frac\pi2x)\sin (n\frac\pi2x)-\sin(\pi x)\sin (n\frac\pi2x) + 4\sin(3\pi x)\sin (n\frac\pi2x) dx \end{align}
giving us
$C_{n} = \int_0^2 2\sin(\frac\pi2x)\sin (n\frac\pi2x) dx -\int_0^2 \sin(\pi x)\sin (n\frac\pi2x)dx + 4\int_0^2\sin(3\pi x)\sin (n\frac\pi2x) dx = \int_0^2 2\sin((1)\frac\pi2x)\sin (n\frac\pi2x) dx -\int_0^2 \sin((2)\frac\pi2x)\sin (n\frac\pi2x)dx + 4\int_0^2\sin((6)\frac\pi2x)\sin (n\frac\pi2x) dx $
(Can you see why the sum drops out? Hint: the sum runs through the entire index of integers, but if $n$ and $n'$ are not equal then the integral is $0$)
Simplifying to $C_{n} = \begin{cases} 2, & n = 1 \\ 0, & n \neq 1\end{cases} +\begin{cases} -1, & n = 2 \\ 0, & n \neq 2\end{cases} + \begin{cases} 4, & n = 6 \\ 0, & n \neq 6\end{cases} $
We actually have a finite number of terms in the sum that are not $0$, using $U(x,t) = \sum_{n=1}^\infty C_ne^{\frac {-\pi^2 n^2t} {16}}\sin(n\frac\pi2x)$, and plugging in those constants.
$U(x,t) = 2e^{\frac {-\pi^2 1^2t} {16}}\sin(1\frac\pi2x)-e^{\frac {-\pi^2 2^2t} {16}}\sin(2\frac\pi2x)+4e^{\frac {-\pi^2 6^2t} {16}}\sin(6\frac\pi2x)$
Which simplifies to
$U(x,t) = 2e^{\frac {-\pi^2t} {16}}\sin(\frac\pi2x)-e^{\frac {-\pi^2t} {4}}\sin(\pi x)+4e^{\frac {-9\pi^2 t} {4}}\sin(3\pi x)$
Here is the graph, which I plotted between $x=0$ and $x=2$ with $t$ ranging from $0$ to $1$. you can see that near $t=0$ there are strong fluctuations in $x$ which quickly smooth out. The physical representation of this is large temperature gradients which would come into thermal equilibrium as heat flows using thermodynamic laws (well, as much as a system with two points that are held at temperature zero can be called a thermodynamic system). Eventually the function would be zero across the entire boundary as the limit of time tends to infinity.
Edit Attaching the determination that our independent ODEs cannot equal $\lambda^2$ (I.E. any positive number)
Now are equations are
$ X'' = \lambda^2 X$
$T'= \frac 14\lambda^2T$
giving us
$X(x) = Ae^{ \sqrt\lambda x}$
$T(t) = Ce^{\frac {\lambda^2 t} 4}$
Let us try to apply our boundary conditions of $U(0,t)=0$, $U(2,t)=0$ into $X(x)$
$U(0,t)= Ae^0 = A = 0$
This shows that the constant $A$ must be $0$ so no equations of this form are included. If more than one case, for example positive and negative eigenvalues, are valid you must sum them together with their respective constants, complicating the problem.