I have three PDE's and I tried to solve it through Finite Fourier Transform (FFT) method but could not get the required result. How can I solve it? Is there any other method to solve these equations?
\begin{align*} \frac{\partial\theta_c}{\partial X}-\frac{A_4}{\epsilon(Y^2A_1+YA_2+A_3)\,}\,\frac{\partial^2\theta_c}{\partial Y^2}&=0 \qquad (1) \\ \frac{\partial^2\theta_s}{\partial Y^2}-Bi(\theta_s-\theta_f)&=0 \qquad (2) \\ \kappa\,\frac{\partial^2\theta_f}{\partial Y^2}+Bi(\theta_s-\theta_f)-\kappa\,U_p\,\frac{\partial\theta_f}{\partial X}&=0 \qquad (3) \end{align*}
Boundary Conditions:
\begin{align*} \text{at}\quad Y&=0,\qquad \theta_s=\theta_f,\qquad 1=-\kappa\,\frac{\partial\theta_f}{\partial Y}-\frac{\partial\theta_s}{\partial Y}\\ \text{at}\quad Y&=Y_p, \qquad \theta_s=\theta_f=\theta_c,\qquad \gamma=-\frac{\kappa}{\epsilon}\,\frac{\partial\theta_c}{\partial Y}=-\kappa\,\frac{\partial\theta_f}{\partial Y}-\frac{\partial\theta_s}{\partial Y} \\ \text{at}\quad Y&=1,\qquad \frac{\partial\theta_c}{\partial Y}=0 \\ \text{at}\quad X&=0,\qquad \theta_s=\theta_f=\theta_c=0. \end{align*} Note that $\epsilon, \kappa, Bi, U_p, A_1, A_2, A_3, A_4$ are constants.
I want to solve for $\theta_s,\;\theta_f,$ and $\theta_c.$