Given a continuous function f(x) on the interval [0, 1] that satisfies the boundary conditions f(0) = f(1) = 0, use the method of separation of variables to solve the partial differential equation:
    ∂u/∂t = k ∂²u/∂x²
Provide step-by-step reasoning and explain your method.
To solve the given partial differential equation (PDE) using separation of variables, we assume that the solution u(x,t) can be written as a product of two functions, one depending only on x and the other depending only on t. That is,

u(x,t) = X(x)T(t)

Now let's substitute this into the PDE:

∂u/∂t = k ∂²u/∂x²

Differentiating with respect to t and x respectively, we get:

X(x)dT/dt = k T(t)d²X/dx²

Now, we separate the variables by dividing both sides by XT:

(1/T)dT/dt = (k/X)d²X/dx²

Notice that the left side depends only on t, while the right side depends only on x. For these expressions to be equal for all x and t, they must both be equal to a constant value, say -λ^2:

(1/T)dT/dt = -λ^2
(k/X)d²X/dx² = -λ^2

We now have two ordinary differential equations (ODEs):

dT/dt = -λ^2T          (1)
d²X/dx² = -(λ^2/k)X    (2)

Let's solve each ODE separately.

For Equation (1), it's a first-order linear ODE:

dT/dt + λ^2T = 0

The general solution is:

T(t) = Ce^(-λ^2t), where C is an arbitrary constant.

For Equation (2), it's a second-order linear ODE with boundary conditions X(0)=X(1)=0:

d²X/dx² + (λ^2/k)X = 0

This is an eigenvalue problem. The general solution depends on the value of λ^2/k:

1. If λ^2/k > 0, the general solution is X(x) = A*cos(√(λ^2/k)x) + B*sin(√(λ^2/k)x). However, this case doesn't satisfy the boundary conditions.

2. If λ^2/k < 0, the general solution is X(x) = A*cosh(√(-λ^2/k)x) + B*sinh(√(-λ^2/k)x). This case also doesn't satisfy the boundary conditions.

3. If λ^2/k = 0, the general solution is X(x) = Ax + B. But this only satisfies the boundary condition if A=B=0, which leads to a trivial solution (u=0).

Therefore, we need to consider other forms of solutions that can satisfy the given boundary conditions. One such form is the Fourier sine series:

X(x) = sin(nπx), where n is an integer.

Substituting this into Equation (2):

d²X/dx² + (n²π²/k)X = 0

This satisfies the ODE and the boundary conditions for any positive integer value of n. The corresponding time-dependent part is T(t) = Ce^(-n²π²kt).

Thus, the complete solution u(x,t) is a linear combination of these product functions:

u(x,t) = Σ [C_n * sin(nπx) * e^(-n²π²kt)]

where C_n are constants determined by the initial condition or other constraints on the problem.