# Solving matrix system of ODE coupled to a non linear PDE (Reynolds)

I am trying to solve the following problem for $$p$$ (the pressure) between a piston skirt and a cylinder lining. This depends on the secondary motion of top and bottom of the piston skirt $$e_t$$ and $$e_b$$ respectively.

These values are given by $$\begin{pmatrix} m_{pin}(1-a/l)+m_{pis}(1-b/l) && m_{pin}a/l+m_{pis}b/l\\I_{pis}/l+m_{pis}(a-b)(1-b/l) && m_{pis}(a-b)(b/l)-I_{pis}/l \end{pmatrix}\begin{pmatrix}\ddot{e_t}\\\ddot{e_b}\end{pmatrix}=\begin{pmatrix}F+F_s+F_ftan(\phi)\\M+M_s+M_f\end{pmatrix}$$

However, $$F$$,$$F_f$$, $$M$$ and $$M_f$$ depend on the pressure, which depends on $$e_t$$ and $$e_b$$. For example $$F$$ is given by the integral of the pressure over the surface of the skirt that is in contact with the oil ($$A$$= {$$0\leq y\leq L, \quad2\pi-\theta_1/2\leq\theta\leq\theta_1/2$$ }:

$$F=R\int\int_A p(\theta,y)\cos(\theta)d\theta dy$$

where the pressure is calculated from the reynolds lubrication equation:

$$\frac{\partial}{\partial x}\left(\phi_x(h)h^3\frac{\partial p}{\partial x}\right)+\frac{\partial}{\partial y}\left(\phi_y(h)h^3\frac{\partial p}{\partial y}\right)= 6\mu U\left(\phi_c(h)\frac{\partial h}{\partial x}+\sigma\frac{\partial \phi_s(h)}{\partial x}\right)+12\mu\phi_c(h)\frac{\partial h}{\partial t}$$

With $$x=R\theta$$ and BC:$$p=0 \quad\theta_1\leq\theta\leq\theta_2\\\frac{\partial p}{\partial x}|_{\theta=0}=\frac{\partial p}{\partial x}|_{\theta=\pi}=0\\p(\theta,0)=p(\theta,L)=0$$

The coefficients $$\phi_i(h)$$ all depend on $$h$$ (complicated and long expression which we can call $$f(h)$$ for now. The coupling is that $$h(\theta,y,t)=C+e_t(t)\cos(\theta)+[e_b(t)-e_t(t)]\frac{y}{L}\cos(\theta)$$ and $$\frac{\partial h}{\partial t}=\cos(\theta)[\dot{e_t}(1-y/L)+\dot{e_b}y/L]$$

Is this possible to solve in Mathematica? And if so what is the best way to go about it? Dsolve or Nsolve? Apologies for the rather crude questions but I am very new to Mathematica, i only started using it yesterday as i had spent an awful lot of time trying to solve this problem in Matlab via finite difference method and an iterative scheme but with no success. I played around with Mathematica trying to set up the problem but i can't seem to figure out how to do so. Any help or tips would be greatly appreciated.