I'm trying to solve this system of PDEs. $$\begin{cases}\dfrac AxG^2=\partial_x F-k\cos\theta\\\dfrac AxG(\partial_\theta G)=-\dfrac 1x\partial_\theta F-k\sin\theta\\h(\theta)=\displaystyle\frac1{x_0l}\int_{x_0}^{x_0+h(\theta)}Gdx\end{cases}$$ Here, $x\in(x_0,x_0+h(\theta))$ and $\theta\in[0,2\pi)$ are the variables (they describe the positions of a particle in the physical system I'm interested in), $F,G$ are functions of $x$ and $\theta$, and $h$ is a function of $\theta$ but not of $x$. $x_0$ is a particular value of $x$ and is constant, and $A$, $k$ and $l$ are also constants. Furthermore, the functions $F,G$ are subject to the constraint that $F(x_0+h(\theta),\theta)=F_0$ is constant for all $\theta$, and $G(x_0,\theta)=x_0l$ for all $\theta$. The objective is to solve for $F,G$.
For context, this came up in a physics problem where the PDEs $F,G$ represent certain physical quantities in a physical system I'm trying to model mathematically. $(x,\theta)$ gives the position vector of the particle in question, and $A,k,l$ are relevant physical parameters which come up naturally in the system.
Obviously, this system of equations is very complicated, so most easy ways of analysing PDEs won't work. In case an analytical solution isn't possible, I would also be greatly appreciative of a good numerical method to obtain numerical values for $F,G$ given certain initial conditions. Any thoughts are appreciated!