Not sure if you still need an answer, but here goes;

Define g(x) = f(x) - x*y, this makes the given equation homogeneous;

d^2/dx^2 d/dy g(x, y) = g(x, y)

And now assume g(x, y) = A(x)B(y). That lets you factorize the equation into two parts;

d^2/xx^2 A(x) = a A(x), d/dy B(y) = (1/a) B(y), with some positive constant a.

This can be solved as
A(x) = C1 exp(+-sqrt(a) x), B(y) = C2 exp(y / a)

And by multiplying them back,

g(x, y) = (C1 * C2) exp(+-sqrt(a) x + y / a) = C exp(+-sqrt(a) x + y / a)

Therefore, f(x, y) = C exp(+-sqrt(a) x + y / a) + x*y.

General solution is linear combination of these solution with different values of C and a.

But I might have missed something there, because this solution contains only two arbitrary parameters whereas the original equation is a third order differential equation.

Is this something that even has a symbolic solution? There is some numeric stuff that should be usable with PDEs, but I'm not as familiar with that.

True. It really depends on what you mean by solve... it can be a symbolic solution, an explicit (convergent) Taylor expansion at the origin, a numerical approximation, etc.