OK, I translated into the Wolfram Language and made a code for numerical integration. It is not possible to solve the problem analytically.

A1 = 1; A2 = 1; A3 = 1; A4 = 1; \[Epsilon] = 1; \[Kappa] = 1;

Bi = 1; Subscript[U, p] = 1; eq = {D[Subscript[\[Theta], c][x, y],

x] – A4/(\[Epsilon]*(y^2 A1 + y*A2 + A3))*

D[Subscript[\[Theta], c][x, y], y, y] == 0 ,

D[Subscript[\[Theta], s][x, y], y, y] –

Bi*(Subscript[\[Theta], s][x, y] –

Subscript[\[Theta], f][x, y]) ==

0 , \[Kappa]*D[Subscript[\[Theta], f][x, y], y, y] +

Bi*(Subscript[\[Theta], s][x, y] –

Subscript[\[Theta], f][x, y]) – \[Kappa]*Subscript[U, p]*

D[Subscript[\[Theta], f][x, y], x] == 0 };

ic = {Subscript[\[Theta], s][0, y] == 0,

Subscript[\[Theta], f][0, y] == 0,

Subscript[\[Theta], c][0, y] == 0};

bc = {DirichletCondition[ {Subscript[\[Theta], s][x, y] == 1,

Subscript[\[Theta], f][x, y] == 0,

Subscript[\[Theta], c][x, y] == 0} , y == 1],

DirichletCondition[ {Subscript[\[Theta], s][x, y] == 0,

Subscript[\[Theta], f][x, y] == 0,

Subscript[\[Theta], c][x, y] == 1} , y == 0] };

sol = NDSolve[{eq, ic, bc}, {Subscript[\[Theta], c],

Subscript[\[Theta], s], Subscript[\[Theta], f]}, {x, 0, 1}, {y, 0,

1}]

{Plot3D[Evaluate[Subscript[\[Theta], c][x, y] /. sol], {x, 0, 1}, {y,

0, 1}, Mesh -> None, ColorFunction -> Hue,

AxesLabel -> {“x”, “y”, “”}],

Plot3D[Evaluate[Subscript[\[Theta], s][x, y] /. sol], {x, 0, 1}, {y,

0, 1}, Mesh -> None, ColorFunction -> Hue,

AxesLabel -> {“x”, “y”, “”}],

Plot3D[Evaluate[Subscript[\[Theta], f][x, y] /. sol], {x, 0, 1}, {y,

0, 1}, Mesh -> None, ColorFunction -> Hue,

AxesLabel -> {“x”, “y”, “”}]}