A Mathematica notebook works much better than a jpg image, for the purpose at hand. Just posting in Mathematica InputForm would be better still.

(1) is a common (ie thermal / slab) PDE that any ODE book discusses. conditions upon it reduce the outer equations solutions.

(2) (3) are a hypothesis of a "systems of PDE". neither ODE book I have shows proofs that prove solutions to systems of PDE (which are more variate) being the same as systems of ODE (which can always be reduce in order). though they solve many slab or wave equations.

While I'm not much studied, I think you would be lucky if any of the PDE solvers gave you all solutions because PDE have many more degrees of freedom in solving than do ODE. You'd need a specific trace on a specific surface for an x,y,z solution: but PDE leave certain trace directions as "variable" and that's bound to trip up solvers ability to show all (relevant) sol'n.

Solvers which are quite complicated code - and Wolfram does show this in the Help you can access: this post suggested by Neil helped me understand the overall "problem" with asking "what can NDSolve do and not do".

Advanced NDSolve tutorial

http://reference.wolfram.com/language/tutorial/NDSolveOverview.html

http://community.wolfram.com/groups/-/m/t/1396110?p_p_auth=lf8ySf9M

I checked: https://www.math.psu.edu/tseng/class/Math251/Notes-PDE%20pt1.pdf, but it has no proof or methods than that of inspection. I understand that http://www.it.uu.se/edu/course/homepage/finmet2/vt14/material/Lab3.pdf modeling of 3space directions upon a surface of gradients (and thus possible anti-derivatives of integrals of these) can be done. It uses computer modeling to show possible directions involved in the parent equation of the partial differential, a backtrace analysis.

There is https://web.math.princeton.edu/~seri/homepage/papers/gws-2006-3.pdf

Wikipedia says Cauchy-Riemann equations and "impact upon imaginary" are used in the form of Laplace's equation - showing even a simple single PDE is a touchy thing. Laplace's equation in cylindrical coordinates looks "odd" and it is: it is formed not by simple translation of coord by conditions of preserving IMZ in such ways as Laplace needed (the handling of singularities needed for his solutions). It requires Cauchy-Riemann considerations and "sacrifices" of some areas to arrive at the form.

https://web.math.princeton.edu/~seri/homepage/papers/gws-2006-3.pdf seems to show principles and perhaps parts of proofs for systems of PDE: but I think in a limited way (a few forms where they can be solved).

http://farside.ph.utexas.edu/teaching/jk1/lectures/Electromagnetism.html shows many uses of PDE

i'd say PDE can be reduced to several ODE, because you can introduce trace equations so each is no longer partial. What i mean is simple. If one PDE is a gradient on a surface any (directional derivative allowed) is possible. It's easy to turn a partial equation into a (set of) ordinary differential equation knowing the original equation. But not knowing the equation(s) being solved for: my head is scratching. Analysis of all possible directions will likely result in finding sets of equations. and I do wonder if some proof in some book shows how that is done.

If anyone can suggest to me a free online source that explains PDE and contains proofs (including proofs for system of PDE reduction), I would appreciate it :).

Do you have the equations in a safer format? I don't open Microsoft formatted files. Sometimes they contain viruses.

Kindly see the attached file.

Attachments:

Equations.nb

Ok, can you call $Version? I checked the code on 11.01 and 11.3. It works. But on 10.3 and earlier versions it does not work.

The problem is not solvable in the combination of boundary conditions that you want to use. I will show an example of a solution of a problem with Neumann conditions on two boundaries y = 0, 1:

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] == 
   NeumannValue[1, y == 1] , 
  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] == NeumannValue[1, y == 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] == 0, 
     Subscript[\[Theta], f][x, y] == 1} , y == 1], 
   DirichletCondition[ {Subscript[\[Theta], c][x, y] == 0, 
     Subscript[\[Theta], s][x, y] == 0} , 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", ""}]}

Can you explain to me the physics of this task? What does this system of equations describe? I'll try to come up with a solution method using Mathematica.

In this problem, I have two parallel plates having two regions one is porous medium and the other is clear fluid both are combined. The upper plate is insulated (Y=1) and the bottom plate is applied constant heat flux(Y=0). The interface between two regions is Yp having heat flux which is equal to gamma. The first equation is representing the clear fluid region and the other two equations are for solid and fluid phases in a porous medium under local thermal non-equilibrium condition. I want to solve for temperature profile for both regions.

OK! But this is a completely different task than I thought. The first equation is defined in one area, and two others in another area. The solutions of these equations are conjugate for $y = y_p$.

There must be another model. Here we must solve the non-stationary problem, since the heat flux at $y = y_p$ is unknown beforehand. The solution tends to stationary for a long time. The question of the velocity of a fluid - is it the quadratic Poiseuille profile $A1Y^2+A2Y+A3$? The flow rate in a porous medium is constant $U_p$?

It remains for us to find out what to set for the velocity on the boundary $y = 1$ - the condition $u=0$ or is it a free surface?

Good work. You did a great job, highly appreciated.