Let me demonstrate how to solve the specific case $n=2$ using the method of characteristics. The case arbitrary $n$ is below. The PDE for $f(x_1,x_2)$ is

$$\tag{1} x_1\partial_1f+x_2\partial_2f=f+\lambda (x_1+x_2) $$

Together with initial condition $f(x_1,g(x_1))=\psi(x_1)$ along some curve in the $x_1,x_2$ plane which is the graph of $x_2=g(x_1)$. Parametrize the characteristics with $z$ and differentiate

$$\tag{2} \frac{df}{dz}=\frac{dx_1}{dz}\partial_1f+\frac{dx_2}{dz}\partial_2f:=\text{RHS of (1)} $$

On the characteristic curves we have the following coupled ODEs

$$\tag{3} \frac{dx_1}{dz}=x_1 \\ \frac{dx_2}{dz}=x_2 \\ \frac{df}{dz}=f+\lambda (x_1+x_2) $$

We choose $z=0$ along the initial data curve $f(a_1,g(a_1))=\psi(a_1)$. The first two equations in (3) may be integrated directly

$$\tag{4} x_1(z,a_1)=a_1e^z \\ x_2(z,a_1)=g(a_1)e^z $$

We can substitute these into the last equation in (3)

$$ \frac{df}{dz}-f=\lambda e^z(a_1+g(a_1)) $$

This is a linear, inhomogeneous, first order equation that we may solve directly (eg. using integrating factor or homogeneous plus particular solution), with the initial condition $f(z=0)=\psi(a_1)$

$$\tag{5} f(z,a_1)=\psi(a_1)e^z+\lambda z e^z (a_1+g(a_1)) $$

What remains is to invert the co-ordinate transformation defined in (4) back to $x_1,x_2$ from $z,a_1$. For example, if $g(a)=a^2$, then the solution is

$$ f(x_1,x_2)=\frac{x_1^2}{x_2}\psi\left(\frac{x_2}{x_1}\right)+\lambda(x_1+x_2)\ln\left(\frac{x_1^2}{x_2}\right) $$

Which you can check does solve (1) for arbitrary $\psi$. Now you can generate all kinds of specific solutions by choosing particular forms for $\psi$.

Looking over the calculations, the only change for arbitrary $n$ is the initial condition: $$ f(x_1,x_2,\dots ,x_{n-1},g)=\psi(x_1,x_2,\dots ,x_{n-1}) $$

Where $g=g(x_1,x_2,\dots ,x_{n-1})$. The co-ordinate transformation in (4) becomes

$$ x_i(z,a_1,\dots,a_{n-1})=a_ie^z \qquad, \qquad i\neq n \\ x_n(z,a_1,\dots,a_{n-1})=g e^z $$

And the solution is

$$\tag{6} f(z,a_1,\dots,a_{n-1})=\psi(a_1,\dots,a_{n-1})e^z+\lambda z e^z (a_1+\cdots+a_{n-1}+g) $$

Where again the solution is to be expressed in the original co-ordinates after inverting the co-ordinate transform. If we choose $g=a_1^2$ (for simplicity) then the inversion may be done neatly for arbitrary $n$

$$\tag{7} e^z=\frac{x_1^2}{x_n} \\ a_i=\frac{x_i x_n}{x_1^2} $$

Substituting (7) into (6)

$$ f(x)= \frac{x_1^2}{x_n} \psi(a)\big|_{a=a(x)}+\lambda (x_1+\cdots+x_n) \ln\left( \frac{x_1^2}{x_n}\right) $$