I'll assume that $X$ and $Y$ are non-negative random variables. Let $F_X$ be the cumulative distribution function of $X$ (that is $F_X(t)=\mathbb P(X\le t)$) and $F_Y$ be the cumulative distribution function of $Y$.
In your notation, probably $F_X(t)=\int_0^t p(s)\,ds$ and $F_Y(y)=\int_0^t q(t)\,dt$.
Now define two functions on $[0,1]$: $g_X(x)=\sup\lbrace t\colon \mathbb P(X\le t)\le x\rbrace $
and similarly $g_Y(x)=\sup\lbrace t\colon \mathbb P(Y\le t)\le x\rbrace$. These functions are the increasing rearrangements of $X$ and $Y$. That is these are non-decreasing functions with the property that $m\lbrace x\colon g_X(x)\le t\rbrace =\mathbb P(X\le t)$ and $m\lbrace x\colon g_Y(x)\le t\rbrace = \mathbb P(Y\le t)$.
Now the largest possible value of $\mathbb E XY$ given the distributions is $\int_0^1 g_X(t)g_Y(t)\ dt$. Intuitively the reason for this is that the largest value for the expectation is obtained when the largest values of $X$ are multiplied by the largest values of $Y$. Slightly more precisely imagine you've arranged the $X$ values from largest to smallest. Think of these as "weights" for the $Y$ values. Obviously you get the biggest integral if you weight the big $Y$ values with the biggest weights.