
clear

r = linspace(0,5,31);
th = linspace(0,2*pi,81);
[R TH] = meshgrid(r,th);
T = linspace(0,6,51);
n0 = pi/4;
for i=1:length(T)
    n1 = erf(2*T(i))*pi^(1-5)*T(i)/2+pi*exp(-4*T(i)^2)/4+T(i)*pi^(1.5)/2;   % int(2*Pi*r*exp(-4*(r-t)^2),r=0..infinity);
    z(:,i) = exp(-4*(r-T(i)).^2)*n0/n1;
end
x = R.*cos(TH);
y = R.*sin(TH);

figure(1)
for i=1:length(T)
    mesh(x,y,ones(size(th'))*z(:,i)');
    axis tight
    set(gca,'fontsize',18,'zlim',[0 1])
    M(i) = getframe;
end

movie(M)

figure(2)
for i=1:length(T)
    plot([-r(end:-1:1) r(1:end)],z([end:-1:1 1:end],i),'k','linewidth',1)
    axis tight
    set(gca,'fontsize',18,'xlim',[-5 5],'ylim',[0 1])
    xlabel('$x$','interpreter','latex')
    ylabel('$u$','interpreter','latex')
    N(i) = getframe;
end

movie(N)
