I want to solve this differential equation using Mathematica:

1/2*\[Sigma]X^2*V1''[x] + ((a0*\[Delta])/a1 - a1*x)*V1'[x] -
em*V1[x] + (2*x - 2*ep*(a0 + a1*x)* \[Lambda])/(em - ep) == 0


The general solution of the nonhomogeneous equation should be equal to the sum of the general solution of the homogeneous equation and a particular solution of the nonhomogeneous equation.

The general solution of the homogeneous equation can be found as follows:

SolGeneralHomo =
DSolve[1/2*\[Sigma]X^2*V1''[x] + ((a0*\[Delta])/a1 - a1*x)*V1'[x] -
em*V1[x] == 0, V1[x], x] // FullSimplify


A particular solution of the nonhomogeneous equation can be found using the Variation of Parameters method. More precisely, given the nonhomogeneous equation:

$$a(x) y′′(x) + b(x) y′(x) + c(x) y(x) = f(x)$$

a particular solution is given by:

$$y1(x) \left( \int \frac{y2(x) (-f(x))}{a(x) W(x)} dx \right) + y2(x) \left( \int \frac{y1(x) f(x)}{a(x) W(x)} dx \right)$$

where $$W(x)$$ is the Wronskian:

$$W(x) = y1(x) y2′(x) - y1′(x) y2(x)$$

Therefore, I do exactly this in Mathematica:

a[x_] := 1/2*\[Sigma]X^2;
y1[x_] :=
HermiteH[-(em/a1), (a1^2 x - a0 \[Delta])/(a1^(3/2) \[Sigma]X)];
y2[x_] :=
Hypergeometric1F1[em/(2 a1), 1/2, (a1^2 x - a0 \[Delta])^2/(
a1^3 \[Sigma]X^2)];
f[x_] := (2*x - 2* \[Lambda]*ep*(a0 + a1*x))/(em - ep);
W[x_] := Wronskian[{y1[x], y2[x]}, x];
SolParticularNonHomo = -y1[
x]*(\[Integral](y2[x]*f[x])/(a[x]*W[x]) \[DifferentialD]x) +
y2[x]*(\[Integral](y1[x]*f[x])/(a[x]*W[x]) \[DifferentialD]x)


It follows that the general solution of the homogeneous equation is equal to:

FinalSol = (V1[x] /. SolGeneralHomo[[1]]) + SolParticularNonHomo //
FullSimplify


If I test the solution, I can see that it does not satisfy the differential equation:

SolvedGeneralNonHomo =
1/2*\[Sigma]X^2*D[FinalSol, {x, 2}] + ((a0*\[Delta])/a1 - a1*x)*
D[FinalSol, x] - em*(FinalSol) + (
2*x - 2*ep*(a0 + a1*x)* \[Lambda])/(em - ep) // FullSimplify;
SolvedGeneralNonHomo == 0


The output is indeed:

(4 (x - ep (a0 + a1 x) \[Lambda]))/(em - ep) == 0


This is equal to f[x] + f[x], which makes me think that there is a sign error. I really can't see where I made a mistake. Could anyone help please?

You got wrong sign when you defined $$f$$. When moving it to RHS it picks a MINUS sign. You did not do that. You wrote

f[x_] := (2*x - 2* \[Lambda]*ep*(a0 + a1*x))/(em - ep);


It should be

f[x_] := - (2*x - 2* \[Lambda]*ep*(a0 + a1*x))/(em - ep);


Here is full code

    ClearAll["Global*"]
ode = 1/2*σ*X^2*V1''[x] + ((a0*δ)/a1 - a1*x)*V1'[x] -
em*V1[x] == -(2*x - 2*ep*(a0 + a1*x)*λ)/(em - ep)
SolGeneralHomo = DSolveValue[ode[[1]] == 0, V1[x], x] // FullSimplify
a = 1/2*σ*X^2;
y1 = First@Cases[SolGeneralHomo, C[1]*any_ :> any]
y2 = First@Cases[SolGeneralHomo, C[2]*any_ :> any]
f = ode[[2]]
W = Wronskian[{y1, y2}, x];
SolParticularNonHomo = (-y1*Integrate[y2*f/(a*W), x] +
y2*Integrate[y1*f/(a*W), x])
FinalSol = SolGeneralHomo + SolParticularNonHomo // FullSimplify
trialSol = V1 -> Function[{x}, Evaluate[FinalSol]]
ode /. trialSol // FullSimplify


• σ*X^2 should be (σX)^2 (i.e. it's a variable named with two letters) or, for ease of notation, let's rename it σ^2. If I make this adjustment the code no longer works. Commented Apr 20, 2023 at 18:35
• After several attempts I found out that replacing σ*X^2 with σX works but (σX)^2 does not work. Commented Apr 20, 2023 at 18:46
• @NC520 I think it is simply that Mathematica can't simplify it to True in the second case. It should not make difference if we replace occurences of σX^2 with (σX)^2. But FullSimplify for some reason does not simplify it to True` in that case. But the main point is that you had minus sign wrong on $f$. Btw, it is better not to use two letter for a single variable in Mathematics. This can be confusing for the reader. Better to use single letter. So you can replace σ*X with say $w$ or any other letter. Commented Apr 20, 2023 at 18:50
• Thank you very much for your answer. The big deal was the wrong sign. Thanks a lot ! Commented Apr 20, 2023 at 18:51