np = 500;
L = 50;
alpha2 = 1;
x = linspace(0, 1, np);
nmax = 100;
T0 = 1;


h = figure;
axis tight manual % this ensures that getframe() returns a consistent size
filename = 'cooling4.gif';

nt = 300;
logtmin = -1;  % 10^(-1)
logtmax = 4;   % 10^4
tt = logspace(logtmin, logtmax, nt);

for t = tt
    u = termal(x, t, alpha2, L, nmax);
    plot(x, u);
    xlabel('x/L')        % x-axis label
    ylabel('T') % y-axis label
    teff = alpha2*pi^2*t/L^2;
    title(['\alpha^2 \pi^2 t/L^2 = ',num2str(teff,'%10.6f\n')])
    ylim([0 1.05*L]);
    drawnow();
    pause(0.1);
    
    % Capture the plot as an image 
    frame = getframe(h); 
    im = frame2im(frame); 
    [imind, cm] = rgb2ind(im, 256); 
    
    % Write to the GIF File 
    if t == tt(1) 
        imwrite(imind, cm, filename, 'gif', 'Loopcount', inf); 
    else 
        imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append'); 
    end 
end

function u = termal(x, t, alpha2, L, nmax)
    np = length(x);
    u = L*ones(1, np)/2;
    for m = 0:nmax
        u = u - 4*L/pi^2*cos(pi*(2*m+1)*x)*exp(-(2*m+1)^2*alpha2*pi^2*t/L^2)/(2*m+1)^2;
    end
end
