how to plot the solution to this PDE?

Question

Maple generates a strange solution form for this PDE. I am having hard time plotting the solution.

The solution is in terms of infinite series. I set the number of terms to say 20, and then set the time to plot the solution at t=2 seconds. Then want to plot the solution now for x=0..1. But the plot comes out empty.

When I sample the solution, and use listplot, I get correct solution plot.

Here is MWE

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x$2)+x;
bc:=u(0,t)=0,u(1,t)=0;
ic:=u(x,0)=x*(1-x);
sol:=pdsolve({pde,ic,bc},u(x,t)):
sol:=value(sol);

Now set the number of terms to 20 and set t=2

sol2:=subs(t=2,sol):
sol2:=subs(infinity=20,sol2);

The above is what I want to plot.

plot(rhs(sol2),x=0..1);

I get empty plot

So had to manually sample it and use listplot

f:=x->rhs(sol2);
data:=[seq([x,f(x)],x=0..1,.01)]:
plots:-listplot(data);

Solution looks correct, when I compare it to Mathematica's result. But Mathematica result is simpler as it does not have those integrals in the sum.

pde=D[u[x,t],t]==D[u[x,t],{x,2}]+x;
bc={u[0,t]==0,u[1,t]==0};
ic=u[x,0]==x(1-x);
DSolve[{pde,ic,bc},u[x,t],x,t];
%/.K[1]->n;
%/.Infinity->20;
%/.t->2;

And the plot is

Question is: How to plot Maple solution without manually sampling it?

Answer 1

Short answer seems to be that it is a regression in Maple 2017.3.

For me, your code works directly in Maple 2017.2 and Maple 2016.2 (without any unevaluated integrals). I will submit a bug report against the regression.

[edited] Let me know if any of these four ways work for your version (presumably Maple 2017.3).

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x$2)+x;
bc:=u(0,t)=0,u(1,t)=0;
ic:=u(x,0)=x*(1-x);
sol:=pdsolve({pde,ic,bc},u(x,t)):
sol:=value(sol);

sol5:=value(combine(subs([sum=Sum,t=2,infinity=20],sol))):
plot(rhs(sol5),x=0..1);

sol4:=combine(subs([sum=Sum,t=2,infinity=20],sol)):
(UseHardwareFloats,oldUHF):=false,UseHardwareFloats:
plot(rhs(sol4),x=0..1);
UseHardwareFloats:=oldUHF: # re-instate

sol2:=subs([sum=Sum,int=Int,t=2],sol):
# Switch integration and summation in second summand of rhs(sol).
sol3:=subsop(2=Sum(int(op([2,1,1],rhs(sol2)),op([2,2],rhs(sol2))),
                   op([2,1,2],rhs(sol2))),rhs(sol2)):
# Rename dummy index and combine summations.
sol3:=Sum(subs(n1=n,op([1,1],sol3))+op([2,1],sol3),
          subs(n1=n,op([1,2],sol3))):
# Curtail to first 20 terms.
sol3:=lhs(sol2)=subs(infinity=20,simplify(sol3));
plot(rhs(sol3),x=0..1);

F:=unapply(subs([Sum='add'],rhs(sol3)),x):
plot(F,0..1);

[edited] Here is yet another way, working for me in Maple 2017.3 on 64bit Linux.

It produces the plot quickly, and doesn't involve curtailing any sum at 20 terms. Note that it does not do your earlier step of sol:=value(sol); since it does active int rather than Int before hitting any Sum with value. It also uses an assumption on x corresponding to the plotting range.

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x$2)+x:
bc:=u(0,t)=0,u(1,t)=0:
ic:=u(x,0)=x*(1-x):
sol:=pdsolve({pde,ic,bc},u(x,t)):

solA:=subs(sum=Sum,value(eval(eval(sol,t=2),Int=int))) assuming x>0, x<1;

plot(rhs(solA),x=0..1) assuming x>0, x<1;