You can derive two times the identity y=C1e−x+C2xe−x, take the two derivatives for solving C1 and C2, and then replace these constants in the initial identity, that is:
var("C1, C2")
y = function("y")(x)
f = C1*e^(-x) + C2*x*e^(-x)
sol = solve([diff(y==f,x), diff(y==f,x,2)], C1, C2)
ode = (y-f==0).subs(sol).full_simplify()
show(ode)
The output is
y(x)+2∂∂xy(x)+∂2(∂x)2y(x)=0
You can prettify it by using, for example, the following code inspired by @slelievre's answer to this question, i.e.:
y1 = function("y1", latex_name="y'")(x)
y2 = function("y2", latex_name="y''")(x)
ode_pretty = ode.subs({diff(y(x), x):y1, diff(y(x), x, x): y2})
show(ode_pretty)
which yields
y(x)+2y′(x)+y″