restart:

bfn:=(n,x,a,b)->sqrt(2/(a+b))*sin(n*Pi*(x+b)/(a+b));

Tfn:=(n,t,a,b,sigma)->exp(-n^2*Pi^2*sigma^2*t/2/(a+b)^2);

'p(x,t)'=Sum(b[n](0)*b[n](x)*T[n](t),n=1..infinity);

pnfn:=(n,x,t,a,b,sigma)->bfn(n,0,a,b)*bfn(n,x,a,b)*Tfn(n,t,a,b,sigma): 'pnfn(x,t)'=pnfn(n,x,t,a,b,sigma);

Check to make sure the "Heat equation" is satisfied

Diff(pnfn(n,x,t,a,b,sigma),t)-(sigma^2/2)*Diff(pnfn(n,x,t,a,b,sigma),x,x);

diff(pnfn(n,x,t,a,b,sigma),t)-(sigma^2/2)*diff(pnfn(n,x,t,a,b,sigma),x,x);

N:=50: p:=sum(pnfn(n,x,t,a,b,sigma),n=1..N):

a:=1: b:=1: subs(t=0,sigma=1,p): plot(%,x=-b..a);