How to solve this second order linear pde? I have the following pde for $f(t,x,y)$:
$a x^2 f_{xx} + bxf_x + f_t - bxy + c = 0$ subject to $f(T,x,y)=0$ for all positive $x,y$, where  $a,b$ and $c$ are constants.
The equation seems to be simple enough (second order linear pde, with polynomial coeficients of simple form and low degree) but non of the standard methods I now seems to work out well. 
Any ideas? Can we guarantee at least that there is existence and uniqueness of a regular solution?
 A: We can get ridd of the source term by defining 
$$f(x,y,z,t) = xy - ct + u(x,y,z,t)\tag{1}$$
so that the PDE for $u$ becomes
$$ax^2u_{xx} + bx u_x + u_t = 0\tag{2}$$
We can then perform a change of variables $z = \log(x)$ to get
$$au_{zz} + (b-a) u_z + u_t = 0\tag{3}$$
which is linear second order homogenous PDE with constant coefficients and you can for example use separation of variables to solve it. Note that the last step is not really needed if you intend to use separation of variables as this can be applied directly to $(2)$ (but you might need to perform a similar change variables on the resulting ODE to solve it).
A: Let us write f as a multivariate polynomial or power series
$$f(t,x,y) = \sum_{{k,l,m}\in{\mathbb N}^3}c_{klm}x^ky^lt^m$$
The partial differential operators will perform (sum omitted): 
$$\frac{\partial f} {\partial x} = kc_{klm}x^{k-1}y^lt^m \hspace{1cm} \frac{\partial f} {\partial y} = lc_{klm}x^{k}y^{l-1}t^m \hspace{1cm}\frac{\partial f} {\partial t} = mc_{klm}x^{k}y^lt^{m-1}$$
Multiplication with a monomial in x : $x^n$ will give $c_{klm}x^{k+n}y^lt^m$, so all terms will together create a linear equations system for the coefficients $c_{klm}$. 
Now remains the boundary conditions. I will assume that $T$ is a constant value for $t$. Then we see that $$f(T,x,y) = \sum_{k,l,m\in\mathbb N^3}c_{klm}x^ky^lT^m = \sum_{k,l,m\in\mathbb N^3}c_{klm}T^mx^ky^l$$ And we see that setting equal to zero and treating $T^m$ as constants just gives us one more set of linear equations involving all terms with $x^ky^lt^0$.
Now the problem that remains is the book-keeping to build this equation system. How you would want to do this in practice depends on what computer language(s) you are familiar with.
For existance I suppose we can analyze the rank of our equation system for a finite approximation. There will definitely exist at least some norm-minimizing solution for any approximation of finite dimensions due to the existence of pseudo-inverse and as the space of functions is so that each new is linearly independent from the others... However it is possible that we find infinite number of solutions and in that case what we need to do is to introduce some kind of a regularizer to encourage more interesting solutions and/or simpler representation over others.

EDIT
Just to show the reasonability (response to Winther's comment):
If we disregard the boundary condition for a second we can find the (trivial) solution $f = bxy - ct$ by hand. If we exchange the $-bxy$ term to $-bxy^k$ still solves it (still trivial since $f_{xx}$ zero, right?). Here is an example with $-bxy^6$,$a=4$,$b=1$, $c=10$.
 sparse(reshape(v_0,[9,9,9])(:,:,1))
        (2, 7) ->  1.00000
 sparse(reshape(v_0,[9,9,9])(:,:,2))
        (1, 1) -> -10.0000

So the indexes are exponents+1 ( Languages like Matlab and Octave start indexing at 1 instead of 0 ). So (2,7) first line is the coefficient for $x^1y^6t^0$ and the second line is for $x^0y^0t^1$. I don't have any other "ground truth" to compare to, so I guess this will have to do so far. As you can see the candidate polynomial includes coefficients for all up to exponents 8 for all k,l,m. A small regularizer pushes all of the rest coefficients to $0$ or close to the machine constant $10^{-15}$ for our example. 

