clc
animatesolution("stiff")
animatesolution("generic")

function animatesolution(method)

    % Parameters and initial conditions
    delta = 0.001;
    F = @(t, y) y^2 - y^3;
    opts = odeset('RelTol', 1.e-4);

    if method == "generic"
        integrator = @ode45;
        file = 'flame_animation_generic.gif';
    else
        integrator = @ode23s;
        file = 'flame_animation_stiff.gif';
    end

    % Solving ODE
    [T, Y] = integrator(F, [0 2/delta], delta, opts);

    % Set up first frame
    figure('Color', 'white')
    plot(T, Y, 'LineWidth', 1)

    hh1(1) = line(T(1), Y(1,1), 'Marker', '.', 'MarkerSize', 20, 'Color', 'r');
    xlabel('time (a.u.)')
    ylabel('radius (a.u.)')
    ht = title(sprintf('Time: %.0f', T(1)));

    % Get figure size
    pos = get(gcf, 'Position');
    width = pos(3);
    height = pos(4);

    % Preallocate data (for storing frame data)
    mov = zeros(height, width, 1, length(T), 'uint8');

    % Loop through by changing XData and YData
    for id = 1:length(T)
        % Update graphics data. This is more efficient than recreating plots.
        set(hh1(1), 'XData', T(id), 'YData', Y(id, 1))
        set(ht, 'String', sprintf('Time: %0.0f', T(id)))

        % Get frame as an image
        f = getframe(gcf);

        % Create a colormap for the first frame. For the rest of the frames,
        % use the same colormap
        if id == 1
            [mov(:,:,1,id), map] = rgb2ind(f.cdata, 256, 'nodither');
        else
            mov(:,:,1,id) = rgb2ind(f.cdata, map, 'nodither');
        end
    end

    % Create animated GIF
    imwrite(mov, map, file, 'DelayTime', 0, 'LoopCount', inf)

end
