Solving differential equation with two variables

Using the definition of implicit differentiation for $f(x,y)$ one can write $$\frac{dy}{dx}=-\frac{f_x}{f_y}$$ And from your equation: $$\frac{a x}{f_x}=-\frac{b y+f_y}{f_y}$$ Therefore: $$\begin{align} \frac{dy}{dx}=\frac{a x}{b y+f_y}&\implies\int a x\,dx=\int(b y+f_y)dy\\ &\implies\dfrac 12 ax^2+C=\dfrac 12 by^2+f(x,y)\\ &\implies f(x,y)=\dfrac 12 ax^2-\dfrac 12 by^2+C \end{align}$$ Update:
The boundary condition is $f(x_0,y_0)=c$ which means $$\dfrac 12 (ax_0^2-by_0^2)+C=c\implies C=c-\dfrac 12 (ax_0^2-by_0^2)$$ and finally $$f(x,y)=\dfrac 12 a(x^2-x_0^2)-\dfrac 12 b(y^2-y_0^2)+c$$ But it seems the answer must be $\frac a2(x^2-x_0^2)-\frac b{\color{red}{4}}(y^2-y_0^2)+c$. Don't know why I ended up with that.

$$\dfrac{ax}{\dfrac{\partial{f(x,y)}}{\partial x}} + \dfrac{by}{\dfrac{\partial{f(x,y)}}{\partial y}}=-1 \quad \text{ with condition} \quad f(x_0,y_0)=c$$ This is a non-linear PDE : The given condition isn't sufficient to determine a unique solution.

For example, the method of separation of variables leads to : $$f(x,y)=\frac{a(x^2-x_0^2)}{2\lambda}-\frac{b(y^2-y_0^2)}{2(\lambda+1)}+c$$ which is a family of solutions of the PDE fitting with the given condition.

This is true any constant $\lambda$, except $0$ and $-1$. So, they are an infinity of solutions.

Note : The solution $\frac a2(x^2-x_0^2)-\frac b{\color{red}{4}}(y^2-y_0^2)+c$ cited by "polfosol" is the particular solution in case of $\lambda=1$.