Here's my (biophysical) system:
https://i.stack.imgur.com/Kchib.jpg
First, the intuition: the disk is like a cell (microorganism) that can swim around inside the (2D) box. $R(x,y,t)$ is the concentration of a resource that the cell consumes. The resource dynamics are purely diffusive and the cell absorbs resource in proportion to its local concentration, thereby acting as a resource sink. The cell also tries to swim towards regions of higher resource concentration by sensing the resource distribution around its circumference and computing its 'centre of mass' (really, centre of concentration),($x_{cr}(t),y_{cr}(t)$). It then applies a thrust force along the vector that runs from the cell's geometric centre ($x_c(t),y_c(t)$) to the 'centre of concentration'. This thrust is also balanced by a drag force. The resource diffusion and cell dynamics are therefore coupled and I want to solve for $R(x,y,t)$, $x_c(t)$ and $y_c(t)$.
The diffusion dynamics are: \begin{align} \frac{\partial R(x,y,t)}{\partial t}&=D\left(\frac{\partial^2 R(x,y,t)}{\partial x^2}+\frac{\partial^2 R(x,y,t)}{\partial y^2}\right)+\frac{\partial R(x,y,t)}{\partial t}\Big\rvert_C \end{align} Here, the cell sink term is assumed to be \begin{align} \frac{\partial R(x,y,t)}{\partial t}\Big\rvert_C= q R(x,y,t)\int_{-a}^a \int_{-\sqrt{a^2-\chi^2}}^{\sqrt{a^2-\chi^2}}\delta(x-x_c(t)-\chi)\delta(y-y_c(t)-\psi) d\psi d\chi\label{sink} \end{align} where $q$ is the resource absorptivity of the cell, $a$ is its radius, and $\delta$ is the Dirac delta function. Intuitively, the cell is a continuous disk of point sinks, the sink rates of which depend on the local resource concentration and the resource absorptivity. The integral evaluates to a unit step function that equals 1 inside the cell and 0 outside.
The cell accelerates per Newton's 2nd law, \begin{align} \frac{d^2x_c(t)}{d t^2}&=\frac{1}{m}\left(k\left(x_{cr}(R(x,y,t),x_c(t),y_c(t))-x_c(t)\right)-\gamma\frac{d x_c(t)}{d t}\right)\\ \frac{d^2y_c(t)}{d t^2}&=\frac{1}{m}\left(k\left(y_{cr}(R(x,y,t),x_c(t),y_c(t))-y_c(t)\right)-\gamma\frac{d y_c(t)}{d t}\right). \end{align} Here, $m$ is the cell mass, and $k$ and $\gamma$ respectively thrust and drag coefficients. The 'centre of resource' coordinates, $(x_{cr},y_{cr})$ are calculated analogously to the centre of mass of a circular ring of variable mass density:
\begin{align} x_{cr}(R(x,y,t),x_c(t),y_c(t))&=\frac{\int_0^{2\pi}R(x_c(t)+a\cos{\phi},y_c(t)+a\sin{\phi},t)\left(x_c(t)+a \cos \phi\right)d\phi}{\int_0^{2\pi}R(x_c(t)+a\cos{\phi},y_c(t)+a\sin{\phi},t)d\phi}\\ y_{cr}(R(x,y,t),x_c(t),y_c(t))&=\frac{\int_0^{2\pi}R(x_c(t)+a\cos{\phi},y_c(t)+a\sin{\phi},t)\left(y_c(t)+a \sin \phi\right)d\phi}{\int_0^{2\pi}R(x_c(t)+a\cos{\phi},y_c(t)+a\sin{\phi},t)d\phi}. \end{align}
These coordinates may also reasonably be approximated by replacing the integrals with sums over a finite set of $n$ points around the cell's circumference. Assuming these to be spaced at equal intervals, $\Delta\phi$, the $i$th point's coordinates are $(x_c(t)+a \cos\phi_i,y_c(t)+a\sin\phi_i)$, where $\phi_i=((i-1)\Delta\phi),~i=1,2,..n$. The `centre of resource' coordinates are then
\begin{align} x_{cr}(R(x,y,t),x_c(t),y_c(t))&=\frac{\sum_i^nR(x_c(t)+a\cos{\phi_i},y_c(t)+a\sin{\phi_i},t)\left(x_c(t)+a \cos \phi_i\right)}{\sum_i^nR(x_c(t)+a\cos{\phi_i},y_c(t)+a\sin{\phi_i},t)}\\ y_{cr}(R(x,y,t),x_c(t),y_c(t))&=\frac{\sum_i^n R(x_c(t)+a\cos{\phi_i},y_c(t)+a\sin{\phi_i},t)\left(y_c(t)+a \sin \phi_i\right)}{\sum_i^nR(x_c(t)+a\cos{\phi_i},y_c(t)+a\sin{\phi_i},t)}. \end{align}
I've tried solving the system, both numerically and quasianalytically, in Mathematica. It doesn't seem to be able to handle either the exact or approximate formulations. I'd really appreciate any tips on how to attack this thing, either computationally or analytically (or, presumably, both). Thanks!
I don't think you can get hope for an analytic solution.
Numerically, probably the most obvious line of attack is to discretise the box into small squares/rectangles. Then approximate $R$ by piecewise-linear on those small rectangles, or equivalently, linear interpolation between grid-points that are the vertices of your small square/rectangles. From this $R$ you can calculate $x_{cr},y_{cr}$. In the time direction there are lots of choice you can make. For propagating $R$ in time: the diffusion term is obvious, and for the resource consumption you could compute on the grid point according to how much the circle overlap with the "region of control" by the grid point.
There are also other numerical methods you could try. For example, Adomian method or more shiny homotopy analysis method. But they require more setting up.