np = 500;
L = 50;
alpha2 = 1;
x = linspace(0, 1, np);
nmax = 100;
T0 = 1;
Tl = 0.2*T0;
Tr = 0.8*T0;


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

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

for t = tt
    u = termal(x, t, alpha2, L, Tl, Tr, nmax);
    plot(x, u);
    xlabel('x/L')        % x-axis label
    ylabel('T / T_0') % y-axis label
    ylim([-.1 1.05*T0]);
    teff = alpha2*pi^2*t/L^2;
    title(['\alpha^2 \pi^2 t/L^2 = ',num2str(teff,'%10.6f\n')])
    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, Tl, Tr, nmax)
    np = length(x);
    u = zeros(1, np);
    w = zeros(1, np);
    
    sig = 1;
    for m = 1:nmax
        u = u + 2*L*sig*sin(pi*m*x)*exp(-m^2*alpha2*pi^2*t/L^2)/(m*pi);
        w = w + 2*L*sin((2*m-1)*pi*x)*exp(-(2*m-1)^2*alpha2*pi^2*t/L^2)/(pi*(2*m-1));
        sig = -sig;
    end
    u = -u*(Tr - Tl)/L - w*Tl + Tl + (Tr - Tl)*x;
end
