I have the follwoing asymptotic exansion to solve a system of diffusion partial differential equations in a two-medium system (medium 1 and medium 2). Medium 1 and medium 2 are connected and have an interface between them at {x1 = 1 and x2 = 0}.
Concentration in the two media is defined here as:

C_1 = C_10 + (1/K)C_11 + (1/K^2)C_12 + ..., 0 < x1 < 1  
C_2 = (1/k)C_21 + (1/K^2)C_22 + ..., 0 < x2 < 1 

I have the following two governing PDE equations:

del(C_1)/del(t_1) = del^2(C_1)/(del(x_1))^2 - alpha_1*C_1, 0 < x_1 < 1  
del(C_2)/del(t_2) = del^2(C_21)/(del(x_2))^2 - alpha_2*C_21, 0 < x_2 < 1

Initial conditions:

C_1 = 0, t_1 = 0   
C_2 = 0, t_2 = 0

Boundary conditions:

C_1 = 1, x_1 = 0
del(C_2)/del(x_2) = 0, x_2 = 1

Interface conditions:

del(C_1)/del(x_1) = diffLen*del(C_2)/del(x_2), x_1 = 1 and x_2 = 0
C_1 = Parti*C_2, x_1 = 1 and x_2 = 0

Note: alpha1, alpha2, and diffLen are constants.

C1 is a concentration function of x1 and t1, and C2 is a concentration function of x2 and t2

I created the follwoing mathematica script to iteratively solve C10, then C21, then C11. I am getting the same value for C10 and C_11. Is my Mathematica implementation correct? Thanks.

wolframscript

ClearAll["Global`*"]

(* Define symbolic constants *)
alphOne = Symbol["alphOne"];
alphTwo = Symbol["alphTwo"];
diffLen = Symbol["diffLen"]; (* diffusivity&length factor: (D2/D1)*(L/gamm) *)


(****** C10 ******)

(* Define the PDE *)
pdeOneZero = D[ConeZero[x, t], t] == D[D[ConeZero[x, t], x], x] - alphOne*ConeZero[x, t];

(* Define the initial condition *)
initialConditionOneZero = ConeZero[x, 0] == 0;

(* Define the boundary conditions *)
boundaryConditionsOneZero = {
  ConeZero[x, t] == 1 /. x -> 0,
  D[ConeZero[x, t], x] == 0 /. x -> 1
};

(****** END C10 ******)


(****** C21 ******)

(* Define the PDE *)
pdeTwoOne = D[CtwoOne[x, t], t] == D[D[CtwoOne[x, t], x], x] - alphTwo*CtwoOne[x, t];

(* Define the initial condition *)
initialConditionTwoOne = CtwoOne[x, 0] == 0;

(* Define the boundary conditions *)
boundaryConditionsTwoOne = {
  CtwoOne[x, t] == ConeZero[1, t] /. x -> 0,
  D[CtwoOne[x, t], x] == 0 /. x -> 1
};

(****** END C21 ******)


(****** C11 ******)

(* Define the PDE *)
pdeOneOne = D[ConeOne[x, t], t] == D[D[ConeOne[x, t], x], x] -
alphOne*(ConeOne[x, t]);

(* Define the initial condition *)
initialConditionOneOne = ConeOne[x, 0] == 0;

(* Define the boundary conditions *)
boundaryConditionsOneOne = {
  ConeOne[x, t] == 1 /. x -> 0,
  D[ConeOne[x, t], x] == diffLen*(D[CtwoOne[0, t],x]) /. x -> 1
};

(****** END C11 ******)

(* Solve the pdeOneZero with initial and boundary conditions *)
solutionOneZero = DSolve[{pdeOneZero, initialConditionOneZero, boundaryConditionsOneZero}, ConeZero[x, t], {x, t}];

(* Solve the pdeTwoOne with initial and boundary conditions *)
solutionTwoOne = DSolve[{pdeTwoOne, initialConditionTwoOne, boundaryConditionsTwoOne}, CtwoOne[x, t], {x, t}];

(* Solve the pdeTwoOne with initial and boundary conditions *)
solutionOneOne = DSolve[{pdeOneOne, initialConditionOneOne, boundaryConditionsOneOne}, ConeOne[x, t], {x, t}];

(* Print solutions *)
Print["Solution for ConeZero[x, t]: ", ConeZero[x, t] /. solutionOneZero]
Print["Solution for CtwoOne[x, t]: ", CtwoOne[x, t] /. solutionTwoOne]
Print["Solution for ConeOne[x, t]: ", ConeOne[x, t] /. solutionOneOne]