We have y=y₁+y₂ where the whole solution y is split into two parts: y₁ is the auxiliary, characteristic or complementary component, and y₂ is the particular solution.
We treat D as a variable in the quadratic D²-3D+2=0, so (D-1)(D-2)=0. This tells us that y₁=C₁e^x+C₂e^2x.
From this let Y₁=e^x and Y₂=e^2x. We use these to build y₂=u₁Y₁+u₂Y₂ where u₁ and u₂ are functions of x that have to be found. Now we need to work out some Wronskian determinants: W, W₁ and W₂ defined as:
W=Y₁Y₂'-Y₂Y₁'=(e^x)(2e^2x)-(e^2x)(e^x)=2e^3x-e^3x=e^3x.
W₁=(0)Y₂'-(1/(1+e^x))Y₂=-e^2x/(1+e^x).
W₂=Y₁(1/(1+e^x))-(0)Y₁'=e^x/(1+e^x).
u₁'=W₁/W=-e^-x/(1+e^x), u₂'=W₂/W=e^-2x/(1+e^x).
u₁=∫-e^-x/(1+e^x)dx=
∫(e^-x)(-e^-xdx)/(1+e^-x) and
u₂=∫e^-2x/(1+e^x)dx=
-∫(e^-2x)(-e^-xdx)/(1+e^-x).
Let z=1+e^-x then e^-x=z-1 and dz=-e^-xdx.
So u₁=∫(z-1)dz/z=∫(1-1/z)dz=z-ln(z)=1+e^-x-ln(1+e^-x).
Y₁u₁=e^x+1-e^xln(1+e^-x).
u₂=-∫((z-1)²/z)dz=
-∫(z-2+1/z)dz=-z²/2+2z-ln(z)=
-(1+2e^-x+e^-2x)/2+2(1+e^-x)-ln(1+e^-x),
u₂=-½-e^-x-e^-2x+2+2e^-x-ln(1+e^-x)=
3/2+e^-x-e^-2x-ln(1+e^-x).
Y₂u₂=3e^2x+e^x-1-e^2xln(1+e^-x).
Therefore y₂=e^x+1-e^xln(1+e^-x)+3e^2x+e^x-1-e^2xln(1+e^-x).
y₂=2e^x-e^xln(1+e^-x)+3e^2x-e^2xln(1+e^-x).
y=C₁e^x+C₂e^2x+2e^x-e^xln(1+e^-x)+3e^2x-e^2xln(1+e^-x).
Constants can be combined:
y=c₁e^x+c₂e^2x-e^xln(1+e^-x)-e^2xln(1+e^-x).
y=c₁e^x+c₂e^2x+e^xln(e^x/(1+e^x))+e^2xln(e^x/(1+e^x)).
y=c₁e^x+c₂e^2x-xe^xln(1+e^x)-xe^2xln(1+e^x).