# What is wrong with my approach to solving a heat transfer PDE?

I wanna solve the following heat transfer PDE using Mathematica.

$$\qquad u_{xx}=u_{t}$$

with following conditions:

$$\qquad \begin{cases}u(x,0)=sin(x) &00\\u_{x}(0,t)=1\\u_{x}(\pi,t)=-1\end{cases}$$

I used the following codes:

heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == Sin[x];
bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1};
sol = First @ DSolve[{heqn, ic, bc}, u[x, t], {x, t}]
sol = u[x, t] /. sol /. {K[1] -> n, Infinity -> 20}


But Mathematica says:

ReplaceAll::rmix: Elements of {(u^(0,1))[x,t]==(u^(2,0))[x,t],u[x,0]==Sin[x],{(u^(1,0))[0,t]==1,(u^(1,0))[[Pi],t]==-1}} are a mixture of lists and nonlists.

Where is the problem?

• use bc = Sequence[Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1] to fix the error indicated by the message.
– kglr
Commented Jan 17, 2019 at 19:33
• @Nasser sorry i didn't understand what do you mean. Dsolve just gives the PDE in your image. Where is the solution? Commented Jan 17, 2019 at 19:36
• @kglr No difference. Commented Jan 17, 2019 at 19:40
• Kevin, do you still get the same error message?
– kglr
Commented Jan 17, 2019 at 19:45
• @Nasser How about using numerical method?If it does not have a solution, so why when I use NDSolve and plot the results, it draws a plot? Sorry I'm not pro at this topic. Commented Jan 17, 2019 at 19:59

I solved this by hand to confirm Maple solution.

Since the boundary conditions are not homogeneous, we can't use separation of variables. Let the solution be

$$u=v\left( x,t\right) +r\left( x\right)$$ Where $$v\left( x,t\right)$$ is the solution to $$v_{t}=v_{xx}$$ and homogenous B.C. $$v_{x}\left( 0,t\right) =0,v_{x}\left( \pi,t\right) =0$$ and $$r\left( x\right)$$ is any reference solution which only needs to satisfy the nonhomogeneous boundary conditions: $$r^{\prime}\left( 0\right) =1,r^{\prime }\left( \pi\right) =-1$$. By guessing, let $$r\left( x\right) =Ax+Bx^{2}$$. Let see if this satisfies the boundary conditions. $$r^{\prime}=A+2Bx$$. At $$x=0$$ this implies $$1=A$$. Hence $$r=x+Bx^{2}$$. Now $$r^{\prime}=1+2Bx$$. At $$x=\pi$$ this gives $$-1=1+2B\pi$$ or $$B=-\frac{1}{\pi}$$. Therefore $$r\left( x\right) =x-\frac{1}{\pi}x^{2}$$ Substituting $$u=v\left( x,t\right) +r\left( x\right)$$ into the PDE\ $$u_{t}=u_{xx}$$ and noting that $$r^{\prime\prime}\left( x\right) =-\frac{2}{\pi}$$ gives $$$$v_{t}=v_{xx}-\frac{2}{\pi}\tag{1}$$$$ PDE (1) is now solved using eigenfunction expansion. We need to find eigenfunctions and eigenvalues of $$v_{t}=v_{xx}$$ with $$v_{x}\left( 0,t\right) =0,v_{x}\left( \pi,t\right) =0$$. This is known PDE and have eigenfunctions and eigenvalues as follows. For zero eigenvalue, the eigenfunction is an arbitrary constant. Say $$\beta$$. let $$\beta=1$$ since scale is not important. $$\Phi_{0}\left( x\right) =1$$ And for $$n=1,2,3,\cdots$$ \begin{align*} \Phi_{n}\left( x\right) & =\cos\left( \sqrt{\lambda_{n}}x\right) \\ & =\cos\left( nx\right) \end{align*} with eigenvalues $$\lambda_{n}=n^{2}$$ for $$n=1,2,3,\cdots$$. Now we can eigenfunction expansion and assume the solution to (1)\ is $$$$v\left( x,t\right) =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \tag{2}$$$$ Plugging this into the PDE (1) gives $$\sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi_{n}^{\prime\prime}\left( x\right) -\frac{2}{\pi}$$ But $$\Phi_{n}^{\prime\prime}\left( x\right) =-\lambda_{n}\Phi_{n}\left( x\right)$$ and the above simplifies to $$\sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) =-\sum_{n=0}^{\infty}A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) -\frac{2}{\pi}$$ Since eigenfunctions are complete, we can expand $$\frac{2}{\pi}$$ using them and the above becomes \begin{align} \sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) & =-\sum_{n=0}^{\infty}A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) -\sum_{n=0}^{\infty}C_{n}\Phi_{n}\left( x\right) \nonumber\\ A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) +A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) & =-C_{n}\Phi_{n}\left( x\right) \nonumber\\ A_{n}^{\prime}\left( t\right) +A_{n}\left( t\right) \lambda_{n} & =-C_{n}\tag{3} \end{align} To find $$C_{n}$$ $$\sum_{n=0}^{\infty}C_{n}\Phi_{n}\left( x\right) =\frac{2}{\pi}$$ For $$n=0$$ $$C_{0}\Phi_{0}\left( x\right) =\frac{2}{\pi}$$ But $$\Phi_{0}\left( x\right) =1$$, hence $$C_{0}=\frac{2}{\pi}$$ All other $$C_{m}\,\$$\ for $$m>0$$ are zero. Hence (3) becomes, for $$n=0$$ (since $$\lambda_{0}=0$$) \begin{align*} A_{0}^{\prime}\left( t\right) & =-\frac{2}{\pi}\\ A_{0}\left( t\right) & =-\frac{2}{\pi}t+B_{0} \end{align*} Where $$B_{0}$$ is integration constant. For $$n>0$$ (3) becomes $$A_{n}^{\prime}\left( t\right) +A_{n}\left( t\right) n^{2}=0$$ This has the solution $$A_{n}\left( t\right) =B_{n}e^{-n^{2}t}$$ Where $$B_{n}$$ is constant of integration. Hence from (2) \begin{align*} v\left( x,t\right) & =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \\ & =A_{0}\left( t\right) +\sum_{n=1}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \\ & =-\frac{2}{\pi}t+B_{0}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \end{align*} Since $$u=v\left( x,t\right) +r\left( x\right)$$ then the solution becomes $$$$u\left( x,t\right) =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t+B_{0}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \tag{4}$$$$ At $$t=0$$ $$$$\sin\left( x\right) =\left( x-\frac{1}{\pi}x^{2}\right) +B_{0}+\sum _{n=1}^{\infty}B_{n}\cos\left( nx\right) \tag{5}$$$$ case $$n=0$$ $$\int_{0}^{\pi}\sin\left( x\right) \cos\left( \sqrt{\lambda_{0}}x\right) dx=\int_{0}^{\pi}\left( x-\frac{1}{\pi}x^{2}\right) \cos\left( \sqrt{\lambda_{0}}x\right) dx+\int_{0}^{\pi}B_{0}\cos\left( \sqrt {\lambda_{0}}x\right) dx$$ But $$\lambda_{0}=0$$ hence \begin{align*} \int_{0}^{\pi}\sin\left( x\right) dx & =\int_{0}^{\pi}\left( x-\frac {1}{\pi}x^{2}\right) dx+\int_{0}^{\pi}B_{0}dx\\ 2 & =\frac{\pi^{2}}{6}+B_{0}\pi\\ B_{0} & =\frac{2}{\pi}-\frac{\pi}{6} \end{align*} For $$n>0$$, Multiplying both sides of (5) by $$\cos\left( mx\right)$$ and integrating $$\int_{0}^{\pi}\sin\left( x\right) \cos\left( mx\right) dx=\int_{0}^{\pi }\left( x-\frac{1}{\pi}x^{2}\right) \cos\left( mx\right) dx+\sum _{n=1}^{\infty}B_{n}\int_{0}^{\pi}\cos\left( nx\right) \cos\left( mx\right) dx$$ For $$m=1$$ \begin{align*} 0 & =0+B_{1}\frac{\pi}{2}\\ B_{1} & =0 \end{align*} For $$m>1$$ \begin{align*} -\frac{1+\left( -1\right) ^{m}}{m^{2}\left( -1+m^{2}\right) } & =\frac{\pi}{2}B_{m}\\ B_{m} & =\frac{-2}{\pi}\left( \frac{1}{m^{2}}\frac{\left( -1\right) ^{m}+1}{m^{2}-1}\right) \end{align*} Hence solution (4) becomes \begin{align*} u\left( x,t\right) & =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t-\frac{\pi}{6}+\frac{2}{\pi}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \\ u\left( x,t\right) & =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t-\frac{\pi}{6}+\frac{2}{\pi}+\sum_{n=2}^{\infty}\frac{-2}{\pi}\left( \frac{1}{n^{2}}\frac{\left( -1\right) ^{n}+1}{n^{2}-1}\right) e^{-n^{2} t}\cos\left( nx\right) \end{align*}

Maple solution verified:

heqn := diff(u(x, t), t) = diff(u(x, t), x$2): ic := u(x, 0) = sin(x): bc := eval(diff(u(x,t),x),x=0)=1, eval( diff(u(x,t),x),x=Pi)=-1: sol := pdsolve({heqn, ic, bc}, u(x, t))  Verification of hand solution (used 10 terms in sum, good enough): mysol=(x-1/Pi x^2)-2/Pi t-Pi/6+2/Pi-2/Pi Sum[((-1)^n+1)/(n^2(n^2-1)) Cos[n x]Exp[-n^2 t],{n,2,10}]; D[mysol,x]/.t->0/.x->0 (*1*) D[mysol,x]/.t->0/.x->Pi (*-1*) Plot[{mysol/.t->0},{x,0,Pi},AxesOrigin->{0,0}]  Animation of hand solution Manipulate[ Plot[Evaluate[mysol /. t -> time], {x, 0, Pi}, PlotRange -> {{0, Pi}, {-0.5, 1}}], {{time, 0, "time"}, 0, 1, .1}, TrackedSymbols :> {time}, Initialization :> {mysol = (x - 1/Pi x^2) - 2/Pi t - Pi/6 + 2/Pi - 2/Pi Sum[((-1)^n + 1)/(n^2 (n^2 - 1)) Cos[n x] Exp[-n^2 t], {n, 2, 10}]} ]  ref (1): boundary_values • 1++,Great answer . Commented Jan 18, 2019 at 5:48 • Changing the coefficient of t is not the same as adding any number to the solution. And indeed, changing that coefficient does not effect the bc's or ic's, but it breaks the pde. Commented Jan 18, 2019 at 7:36 • The premise that only having neumann boundary conditions makes the solution non-unique does not apply to the heat equation. There are numerous examples all over the web. One is found at ramanujan.math.trinity.edu/rdaileda/teach/s12/m3357/lectures/…. For some pde's the premise is correct, but different pde's have different bc requirements. By the way, separation of variables works very well in this case, but you do have to use plus separation as well as times separation to get all the pieces. Commented Jan 18, 2019 at 19:42 • The B1 must be zero. How it is 2/Pi?! Commented Jan 21, 2019 at 20:30 • @Nasser When using [0 , Pi] interval in integrating, there should be a 2 factor because we changed [-Pi, Pi] to [0, Pi] so the Mathematica is right. I don't know why you didn't used the 2 factor. Commented Jan 21, 2019 at 21:12  MAPLESOL = (x - 1/Pi x^2) - 2/Pi t - Pi/6 + 2/Pi - 2/Pi Inactivate[ Sum[((-1)^n + 1)/(n^2 (n^2 - 1)) Cos[n x] Exp[-n^2 t], {n, 2, Infinity}]]; Plot3D[MAPLESOL /. Infinity -> 20 // Activate, {x, 0, Pi}, {t, 0, 1}]  Check equation,initial and boundary conditions:  (D[MAPLESOL /. Infinity -> 20, t] == D[MAPLESOL /. Infinity -> 20, {x, 2}]) // Activate (*True*) D[MAPLESOL /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> 0 (* 1 *) D[MAPLESOL /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> Pi (* -1 *) Plot[{Sin[x], Evaluate[(MAPLESOL /. Infinity -> 20 // Activate) /. t -> 0]}, {x, 0, Pi}, PlotStyle -> {Red, {Dashed, Black}}] Plot3D[D[MAPLESOL /. Infinity -> 20 // Activate, x] // Evaluate, {x, 0, Pi}, {t, 0, 1}]  It's seems symbolic solution by Maple is OK. • good answer. Yes, this seems like a valid solution. But it is not unique solution as adding any number to the solution is also a solution. D[(10 + MAPLESOL) /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> 0 gives 1 as well. May be Mathematica did not solve it because of this? I do not know. Maple 2018 generally can solve more PDE's than Mathematica 11.3 Commented Jan 18, 2019 at 0:57 • Adding a constant to the stated solution will break the ic u[x,0] = Sin[x]. I think that the given ic is enough to make the solution unique.When you work the solution with separation of variables there are no extra constants lying around. Commented Jan 18, 2019 at 1:10 • @Nasser and Mariusz …This solution is incorrect, just compare it to the numeric solution. Commented Jan 18, 2019 at 10:30 This question is strongly related to, if not a duplicate of this question. It can be solved with the help of finite Fourier cosine transform and its inversion: heqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ic = u[x, 0] == Sin[x]; bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1}; help[index_] := Module[{tset = finiteFourierCosTransform[{heqn, ic}, {x, 0, Pi}, index] /. Rule @@@ bc /. HoldPattern@finiteFourierCosTransform[f_, __] :> f}, tsol = DSolve[tset, u[x, t], t][[1, 1, -1]]] tsolgeneral = help[n] // Simplify tsolzero = help[0] tsolfunc[n_] = Piecewise[{{tsolgeneral, n != 0}}, tsolzero] sol = inverseFiniteFourierCosTransform[tsolfunc[n], n, {x, 0, Pi}] sol = sol /. HoldForm@Sum[expr_, {n, C}] :> sum[Simplify[expr, n > 1], {n, 2, C}] /. sum -> (HoldForm@Sum@## &)  Notice finite Fourier cosine transform for $$n=0$$ has been calculated separately because currently finiteFourierCosTransform, which is built on Integrate, isn't strong enough to obtain the complete solution for general n. Finally let's compare the solution to the numeric one: solapprox = Compile[{t, x}, #] &[sol /. C -> 10 // ReleaseHold]; solnum = NDSolveValue[{heqn, ic, bc}, u, {t, 0, 1}, {x, 0, Pi}] Manipulate[Plot[{solnum[x, t], solapprox[t, x]}, {x, 0, Pi}], {t, 0, 1}]  # Addendum The solution found above looks different from the one found by method of separation of variables shown in other answers, but it's possible to prove they're identical i.e. sol= -((2 (-1 + t))/Pi) + (2 HoldForm[Sum[-(((1 + (-1)^n) (1 + E^(n^2 t) (-1 + n^2)) Cos[n x])/(E^(n^2 t) (n^2 (-1 + n^2)))), {n, 2, C}]])/Pi  is equivalent to mysol=(x-1/Pi x^2)-2/Pi t-Pi/6+2/Pi- 2/Pi Inactivate[Sum[((-1)^n+1)/(n^2(n^2-1)) Cos[n x]Exp[-n^2 t],{n,2,Infinity}]];  We first find the difference between the summand and calculate the sum: summandseparate = (((-1)^n + 1) Cos[n x] Exp[-n^2 t])/(n^2 (n^2 - 1)); summandtransform = ((1 + (-1)^n) E^(-n^2 t) (1 + E^(n^2 t) (-1 + n^2)) Cos[n x])/( n^2 (-1 + n^2)); diff = Sum[summandseparate - summandtransform, {n, 2, Infinity}] // FullSimplify (* 1/4 (-PolyLog[2, E^(-2 I x)] - PolyLog[2, E^(2 I x)]) *)  Remark At least in v11.2, you need to modify definition of diff to diff = Sum[summandseparate - summandtransform // Simplify // Evaluate, {n, 2, Infinity}]  or the calculation will be extremely slow. Then the problem boils down to proving the following identity: eq = (x - x^2/π) - (2 t)/π - π/6 + 2/π - 2/Pi diff2 == -(( 2 (-1 + t))/π); -4 diff == -4 diff2 /. First@Solve[eq, diff2] (* PolyLog[2, E^(-2 I x)] + PolyLog[2, E^(2 I x)] == 1/3 (π^2 - 6 π x + 6 x^2) *)  which seems to be beyond the reach of Mathematica, so I asked it in math.SE and get the solution by hand in 6 minutes. Please check this post for more details. This identity can be proved using ReducePiecewise: $Assumptions = 0 < x < Pi;

ReducePiecewise[
FullSimplify@
PiecewiseExpand@
PowerExpand@
FunctionExpand[PolyLog[2, E^(-2 I x)] + PolyLog[2, E^(2 I x)]], x, $Assumptions] (* 1/3 (π^2 - 6 π x + 6 x^2) *)$Assumptions =.;

• -4 Series[diff, {x, 0, 10}] // FullSimplify[#, 0 < x < \[Pi]] & gives you what you want. I guess your tails are harder to get rid of because your series converges much more slowly than the others because your Exp[-n^2 t] only damps out part of the sum while the sum damps out very rapidly in the separate variable case. I still don't understand why I can't prove your solution satisfies the pde. Probably need to know a closed form for the whole sum to do that. Commented Jan 20, 2019 at 0:53
• @Bill Yeah Mathematica can verify the identity in experimental mathematics style. We can prove the solution satisfy the PDE in the same manner as shown in Addendum, the key point is to sum part of the summand, the following code takes about 100 seconds: terms = List @@ sol /. HoldForm@Sum[a_, _] :> Sequence@a; sum = Sum[#, {n, 2, Infinity}] &; {Simplify@D[terms, t], {#, sum /@ Expand@#2} & @@ D[terms, x, x] // Simplify} // AbsoluteTiming Commented Jan 20, 2019 at 6:32
• See my update. We don't need the math guys for this. Commented Jan 20, 2019 at 22:23
• @BillWatts I find a simpler way to prove the identity, see my update. Commented Jul 19, 2020 at 14:02
• Nice, but I get the same result with PowerExpand@ FunctionExpand[ PolyLog[2, E^(-2 I x)] + PolyLog[2, E^(2 I x)]] // Simplify. I don't know what the Piecewise functions are doing. Commented Jul 26, 2020 at 7:29

And yet another way.

pde = D[u[x, t], t] == D[u[x, t], {x, 2}]

ic = u[x, 0] == Sin[x]

bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1}


Separate variables, first with times separation.

u[x_, t_] = X[x] T[t]

pde/u[x, t] // Expand
T'[t]/T[t] == X''[x]/X[x]


Each side must be equal to a constant. First try 0.

t0eq = T'[t]/T[t] == 0;

DSolve[t0eq, T[t], t] // Flatten
(*{T[t] -> C[1]}*)

t0 = 1


since we will combine with constants from the x equation.

x0eq = X''[x]/X[x] == 0

DSolve[x0eq, X[x], x] // Flatten
(*{X[x] -> C[2] x + C[1]}*)

x0 = X[x] /. % /. {C[1] -> c1, C[2] -> c2}


Now use a negative constant.

t1eq = T'[t]/T[t] == -α^2

DSolve[t1eq, T[t], t] // Flatten
(*{T[t] -> C[1] E^(α^2 (-t))}*)

t1 = T[t] /. % /. C[1] -> 1

x1eq = X''[x]/X[x] == -α^2

DSolve[x1eq, X[x], x] // Flatten
(*{X[x] -> C[2] Sin[α x] + C[1] Cos[α x]}*)

x1 = X[x] /. % /. {C[1] -> c3, C[2] -> c4}


We don't have all the necessary pieces yet, so do plus separation.

u[x_, t_] = X[x] + T[t]

pde
(*T'[t] == X''[x]*)


Again, each side must be equal to a constant. Call it δ.

xpeq = X''[x] == δ

DSolve[xpeq, X[x], x] // Flatten
(*{X[x] -> C[2] x + C[1] + (δ x^2)/2}*)

xp = X[x] /. % /. {C[1] -> c5, C[2] -> c6}

tpeq = T'[t] == δ

DSolve[tpeq, T[t], t] // Flatten
(*{T[t] -> C[1] + δ t}*)

tp = T[t] /. % /. C[1] -> 0


Put all the pieces together.

Clear[c1, c2, c3, c4, c5, c6, α, δ]

u[x_, t_] = x0 t0 + x1 t1 + xp + tp
(*c1 + c2 x + E^(α^2 (-t)) (c3 Cos[α x] + c4 Sin[α x]) + c5 +
c6 x + δ t + (δ x^2)/2*)


c1 and c5 can be combined, c2 and c6 can be combined.

c5 = 0;
c6 = 0;


Look at the first bc

bc[[1]]
(*c2 + α c4 E^(α^2 (-t)) == 1*)


from which we can say

c2 = 1;
c4 = 0;


and the next bc.

bc[[2]]
-(α c3 Sin[π α] E^(α^2 (-t))) + π δ + 1 == -1


Make the Sin 0

Solve[α π == n π, α] // Flatten

from which

α = n;
$Assumptions = n ∈ Integers && n >= 0  Use the rest of the equation to solve for δ c26eq2 = π δ + 1 == -1 δ = δ /. Solve[π δ + 1 == -1, δ][[1]]  So far we have u[x, t] c1 + c3 E^(-n^2 t) Cos[n x] - (2 t)/π - x^2/π + x  Now for the ic. ic (*c1 + c3 Cos[n x] - x^2/π + x == Sin[x]*)  We are going to have an infinite series for n and the n=0 will be a constant term. Since we don't need two of them we can throw out c1. c1 = 0 ic - (-(x^2/π) + x) (*c3 Cos[n x] == x^2/π - x + Sin[x]*)  To solve for c3, multiply by Cos[n x] and integrate. c3*Integrate[Cos[n*x]^2, {x, 0, Pi}] == (1/Pi)*Integrate[x^2*Cos[n*x], {x, 0, Pi}] - Integrate[x*Cos[n*x], {x, 0, Pi}] + Integrate[Sin[x]*Cos[n*x], {x, 0, Pi}] c3 = c3 /. Solve[%, c3][[1]] // Simplify (*-((2 ((-1)^n + 1))/(\[Pi] n^2 (n^2 - 1)))*)  The constant term for n = 0 needs to be done separately. c30*Integrate[1, {x, 0, Pi}] == (1/Pi)*Integrate[x^2*1, {x, 0, Pi}] - Integrate[x*1, {x, 0, Pi}] + Integrate[Sin[x]*1, {x, 0, Pi}] c30 = c30 /. Solve[%, c30][[1]] // Expand 2/π-π/6  Without the sum on n we have u[x_, t_] = u[x, t] + c30 -((2 ((-1)^n + 1) E^(-n^2 t) Cos[n x])/(π n^2 (n^2 - 1))) - ( 2 t)/π - x^2/π + x - π/6 + 2/π  Check against the pde. pde (*True*)  We can see that odd n the terms are zero, so change to. -((2 ((-1)^n + 1) E^(-n^2 t) Cos[n x])/(π n^2 (n^2 - 1))) /. n -> 2 n // Simplify (*-((E^(-4 n^2 t) Cos[2 n x])/(π n^2 (4 n^2 - 1)))*)  We have an infinite series of n starting from 1, but lets make a finite series. u[x_, t_, mm_] := 2/Pi - Pi/6 - (2*t)/Pi - x^2/Pi + x - (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), {n, 1, mm}]  Check against xzczd's numerical solution. gifs = Table[ Plot[{u[x, t, 20], solnum[x, t]}, {x, 0, π}, PlotRange -> {-1, 1}], {t, 0, 1, .01}]; ListAnimate[%]  This solution matches Maple's and Nasser's solution and the numerical solution. I could not get xzczd's analytic solution to satisfy the pde, which is I think is the cause for the little tails at the end points. Update on showing the equivalence in the separate variable solution and xzczd's fourier transform solution. The separate variable solution simplified by letting the original n go to 2n to eliminate the odd n zero's: usep[x_, t_] := -((2*(t - 1))/Pi) - Pi/6 - x^2/Pi + x - (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), {n, 1, Infinity}]  Rewrite as usep[x, t] := 2/Pi - (2*t)/Pi - (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), {n, 1, Infinity}] + Sum[An*Cos[2*n*x], {n, 1, Infinity}]  where we find a series such that for 0 <= x <= π x - x^2/Pi - Pi/6 == Sum[An*Cos[2*n*x], {n, 1, Infinity}];  Find An by multiplying by Cos[2 n x] and integrating. Integrate[(x - x^2/Pi - Pi/6)*Cos[2*n*x], {x, 0, Pi}] == An*Integrate[Cos[2*n*x]^2, {x, 0, Pi}] An = An /. Solve[%, An][[1]] (*-(1/(π n^2))*)  Add to the existing sum -(1/π) (E^(-4 n^2 t) Cos[2 n x])/(n^2 (4 n^2 - 1)) + An Cos[2 n x] // Simplify (*((E^(-4 n^2 t)/(1 - 4 n^2) - 1) Cos[2 n x])/(π n^2)*)  And we end up with the same solution as the fourier series solution. • There're little tails at the end points in the graphic of solapprox because I've truncated the summation at$n=10\$ (C -> 10). If e.g. C -> 60 is used, the tails will be almost invisible. And, it's possible to prove my solution is identical to the one found by method of separation of variables, see my update. Commented Jan 19, 2019 at 12:14