We have the following system that describes the heat conduction in a rectangular region:
$$\begin{cases} u_{xx}+u_{yy}+S=u_t \\
u(a,y,t)=0 \\
u_x(x,b,t)=0 \\
u_y(0,y,t)=0 \\
u(x,0,t) = 0 \\
u(x,y,0) = f(x,y)
\end{cases}
$$
$S$ is source term (independent of time and space ) placed at the origin.
How can I solve this system analytically using separation of variables? or is there another method to solve it(analytically)?
EDIT: I am looking for the solution when $S$ is present because that's what makes the challenge .
EDIT This is as complete an answer as can be given with the ill-specified problem. Please see the bottom for a new discussion of the forcing term.
Interpreting your $S$ term as a source located at the origin $(S=S_0\delta(x)\delta(y))$, we have: $$ u_{xx} + u_{yy} = u_{t} + S_0\delta(x)\delta(y) $$ So, for all points not precisely the origin, we have the homogeneous equation $$ u_{xx} + u_{yy} = u_{t} $$ We don't yet get to what to do exactly at the origin. For now, let us try the separation of variables approach, guessing a form: $$ u(x,y,t) = X(x)Y(y)T(t) $$
Plugging into the original equation (and dropping the variable dependencies for notational simplicity) we have: $$ X''YT + XY''T = XYT' $$ Dividing by $XYT$ everywhere gives $$ \frac{X''}{X} + \frac{Y''}{Y} = \frac{T'}{T} $$ By the usual argument, the left hand side is a function of only the spatial coordinates and the right hand side is only a function of only the time coordinate, yet they are equal to each other. Therefore, each side must be equal to some constant, which we will denote as $-\lambda$. Let us consider the time equation first: $$ \frac{T'}{T} = -\lambda \;\;\Rightarrow\;\; T' + \lambda T = 0 $$ The general solution to this equation (up to a multiplicative constant we'll account for later) is: $$ T(t) = e^{-\lambda t} $$ Now consider the spatial equation $$ \frac{X''}{X} + \frac{Y''}{Y} = -\lambda $$ Moving terms gives $$ \frac{X''}{X} = -\left(\lambda+\frac{Y''}{Y}\right) $$ By the same argument as above, each side must also be a constant, which I will denote as $-k_x^2$. So for the $x$ equation we have: $$ \frac{X''}{X} = -k_x^2 \;\;\Rightarrow\;\; X'' + k_x^2 X = 0 $$ The general solution here is a sinusoid with spatial frequency $k_x$: $$ X(x) = A\cos(k_x x) + B\sin(k_x x) $$ Now we apply the boundary conditions, which, as written in the original question, are under-specified. I'm going to assume you meant two insulating walls (where the temperature gradient normal to the wall vanishes) and two zero temperature walls (where the temperature itself vanishes). For the $x$ dependence, these are: $$ u_x(0,y,t) = 0 \\ u(a,y,t) = 0 $$ Plugging in the separated forms reveals these correspond to $$ X'(0) = 0 \;\;\Rightarrow\;\; -Ak_x\sin(0) + Bk_x\cos(0) = 0 \;\;\Rightarrow\;\; B=0 $$ and $$ X(a) = 0 \;\;\Rightarrow\;\; A\cos(k_x a) = 0 \;\;\Rightarrow\;\; k_x=\frac{(n-1/2)\pi}{a}\;\forall n\in\mathbb{Z} $$ Now we turn to the $Y$ dependence. We have $$ -\left(\lambda+\frac{Y''}{Y}\right) = -k_x^2 \;\Rightarrow\; Y'' + (\lambda-k_x^2)Y = 0 $$ Now define the constant $k_y^2 = (\lambda-k_x^2)$ so that the solution to the above equation is $$ Y(y) = C\cos(k_y y) + D\sin(k_y y) $$ Now applying the assumed boundary conditions $$ u(x,0,t) = 0 \\ u_y(x,b,t) = 0 $$ to the separated form gives $$ Y(0)=0\;\Rightarrow\; C\cos(0) + D\sin(0) = 0 \;\Rightarrow\; C=0 $$ and $$ Y'(b) = 0 \;\Rightarrow\; Dk_y\cos(k_y b) = 0 \;\Rightarrow k_y=\frac{(m-1/2)\pi}{b}\;\forall m\in\mathbb{Z} $$ Now note that $k_x$ and $k_y$ can only take on discrete values as given above, and that $\lambda = k_x^2+k_y^2$ by the definition above as well. So for any two indices $(m,n)$, the following is a solution of the equation $$ u_{m,n} = C_{m,n} e^{-(k_x^2+k_y^2) t}\cos(k_x x)\sin(k_y y) $$ And now by linearity, a superposition of those is also a solution $$ u(x,y,t) = \sum_n^\infty\sum_m^\infty C_{m,n} e^{-(k_x^2+k_y^2) t}\cos(k_x x)\sin(k_y y) $$ There's two final things to consider, the delta function forcing and the initial condition. First, the $t=0$ initial condition $$ u(x,y,0) = \sum_n^\infty\sum_m^\infty C_{m,n} \cos(k_x x)\sin(k_y y) = f(x,y) $$ This tells me the $C_{n,m}$ must be chosen to be like the fourier coefficients of the $f(x,y)$ function, something like $$ C_{n,m} = \int_0^a\int_0^b f(x,y) \cos(k_x x)\sin(k_y y) dx dy $$ So now we have a full solution to the homogeneous problem with $S=0$, or equivalently, one that is valid precisely everywhere except $x=y=0$. Now we need to know how to handle that singularity at the origin. The way we normally handle this is to force a derivative discontinuity into the solution by flipping signs of the solution at a point in such a way that taking laplacian results in an impulse. However, I claim that the boundary conditions given are incompatible with a forcing at the origin. To see why, consider that $u_x(0,y,t) = 0$ says that there is no heat flux through the left wall. However, the source is touching that wall, meaning that the source is touching a perfectly insulating wall; how can heat from a source flow through an insulating wall? It cannot! The source has no effect in this problem because of the boundary conditions given. Therefore, the homogeneous solution given is the solution to the problem.
If you can clarify the boundary conditions and source specification and post in another question, you might get a better answer that helps you know how to solve these problems. For what it's worth, I can see a path forward with a point source if the boundary conditions were such that the modes had the form $\sin(k_x x)\sin(k_y y)$. As it is, we had to guess at what boundary conditions you meant because the ones you gave don't lead to a well-posed problem.