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
The PDE is linear so the general solution is the sum of a homogeneous and particular solution, $u = u_h + u_p$.
The particular solution satisfies
$$ \frac{\partial u_p}{\partial t} - \nu \left(\frac{\partial^2 u_p}{\partial x_1^2} + \frac{\partial^2 u_p}{\partial x_2^2}\right) = f $$
Making the ansatz $u_p(x_1,x_2,t) = T(t)\cos (\sqrt{a}x_1)\cos (\sqrt{b}x_2)$ and subsituting into the PDE, we get
$$T'(t)\cos (\sqrt{a}x_1)\cos (\sqrt{b}x_2) + \nu T(t)(a+b) \cos (\sqrt{a}x_1)\cos (\sqrt{b}x_2) = e^{-\nu(a+b) t}\cos (\sqrt{a}x_1)\cos (\sqrt{b}x_2), $$
and a solution is obtained when
$$T'(t) + \nu(a+b)T(t) = e^{-\nu(a+b) t}$$
Solving the ODE with an integrating factor we get $T(t) = (C+t)e^{-\nu (a+b)t}$, and the particular solution is
$$u_p = (C+t)e^{-\nu (a+b)t}\cos (\sqrt{a}x_1)\cos (\sqrt{b}x_2)$$
As you have not specified an initial condition, nothing can be said yet about the constant $C$ appearing in the solution.
It remains to solve the homogeneous part,
$$ \frac{\partial u_h}{\partial t} - \nu \left(\frac{\partial^2 u_h}{\partial x_1^2} + \frac{\partial^2 u_h}{\partial x_2^2}\right) = 0$$
As you want to solve on a "bounded domain", once you supply the initial and boundary conditions you can apply separation of variables seeking a solution of the form
$$u_h(x_1,x_2,t) = X_1(x_1)X_2(x_2)T(t)$$
Best Answer
Use $$\frac{\partial u}{\partial s}=\frac{\partial u}{\partial t}\frac{\partial t}{\partial s}$$ Here, $\frac{\partial t}{\partial s}=1/k$, so $$\frac{\partial u}{\partial s} = \frac{1}{k}\frac{\partial u}{\partial t}$$ Substitute for $\frac{\partial u}{\partial t}$, and the answer follows.