function makeMovie(folder_path, run_index) 
    
    set(0,'defaulttextInterpreter','latex')
    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
    
    format long
    persistent yMaxNormFixed
    persistent yMaxEnergyFixed

    % Load data
    FileDir = dir(sprintf(horzcat(folder_path, '/Run_%03i/*.mat'),run_index));
    NumFiles = numel(FileDir);
    
    Data   = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_%i.mat'),run_index,2),'Params','Transf');
    Params = Data.Params;
    Transf = Data.Transf; 
    
    w0 = Params.w0;

    x  = Transf.x; y   = Transf.y; z   = Transf.z;
    dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
    
    mkdir(sprintf(horzcat(folder_path, '/Run_%03i/Figures'),run_index))
    outputVideo           = VideoWriter(fullfile(horzcat(folder_path, '/Movie')), 'MPEG-4');
    outputVideo.Quality   = 100; % Set quality to maximum (0–100)
    outputVideo.FrameRate = 10;
    open(outputVideo)
    
    figure(1);
    clf
    set(gcf,'Position', [100, 100, 1600, 900])
    t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid
    
    for ii = 2:(NumFiles-1)
        % Load data
        Data   = load(sprintf('./Results/Data_3D/RealTimeDynamics/Run_%03i/psi_%i.mat',run_index,ii),'psi','muchem','Observ','t_idx');
        
        if isgpuarray(Data.psi), psi = gather(Data.psi); end
        if isgpuarray(Data.Observ.tVecPlot)
            Observ.tVecPlot = gather(Data.Observ.tVecPlot);
        else
            Observ.tVecPlot = Data.Observ.tVecPlot;
        end
        if isgpuarray(Data.Observ.NormVec)
            Observ.NormVec = gather(Data.Observ.NormVec);
        else
            Observ.NormVec = Data.Observ.NormVec;
        end
        if isgpuarray(Data.Observ.PCVec)
            Observ.PCVec = gather(Data.Observ.PCVec);
        else
            Observ.PCVec = Data.Observ.PCVec;
        end
        if isgpuarray(Data.Observ.EVec)
            Observ.EVec = gather(Data.Observ.EVec);
        else
            Observ.EVec = Data.Observ.EVec;
        end
    
        muchem = Data.muchem;

        if isempty(yMaxNormFixed)
            yMaxNormFixed   = 1.05 * Observ.NormVec;  % Set once in the first iteration
        end

         if isempty(yMaxEnergyFixed)
            yMaxEnergyFixed = 1.05 * Observ.EVec;  % Set once in the first iteration
        end
    
        % Compute probability density |psi|^2
        n = abs(psi).^2;
    
        %Plotting
        nexttile;
        nxz = squeeze(trapz(n*dy,2));
        nyz = squeeze(trapz(n*dx,1));
        nxy = squeeze(trapz(n*dz,3));
        plotxz = pcolor(x,z,nxz');
        set(plotxz, 'EdgeColor', 'none');
        cbar1 = colorbar;
        cbar1.Label.Interpreter = 'latex';
        % cbar1.Ticks = []; % Disable the ticks
        colormap(gca, Helper.Colormaps.plasma())
        xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        title('$|\Psi(x,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
    
        nexttile;
        plotyz = pcolor(y,z,nyz');
        set(plotyz, 'EdgeColor', 'none');
        cbar1 = colorbar;
        cbar1.Label.Interpreter = 'latex';
        % cbar1.Ticks = []; % Disable the ticks
        colormap(gca, Helper.Colormaps.plasma())
        xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        title('$|\Psi(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
    
        nexttile;
        plotxy = pcolor(x,y,nxy');
        set(plotxy, 'EdgeColor', 'none');
        cbar1 = colorbar;
        cbar1.Label.Interpreter = 'latex';
        % cbar1.Ticks = []; % Disable the ticks
        colormap(gca, Helper.Colormaps.plasma())
        xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
        title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
    
        nexttile;
        plot(Observ.tVecPlot*1000/w0,Observ.NormVec,'-b')
        ylim([0, yMaxNormFixed]);
        ylabel('$N$', 'FontSize', 14); xlabel('$t$ [$m$s]', 'FontSize', 14);
        title('Normalization', 'FontSize', 14);
        grid on
    
        nexttile;
        plot(Observ.tVecPlot*1000/w0,1-2*Observ.PCVec/pi,'-b')
        ylim([0,1])
        ylabel('$C$', 'FontSize', 14); xlabel('$t$ [$m$s]', 'FontSize', 14);
        title('Coherence', 'FontSize', 14);
        grid on
    
        nexttile;
        plot(Observ.tVecPlot*1000/w0,Observ.EVec,'-b')
        ylim([0, yMaxEnergyFixed]);
        ylabel('$E_{tot}$', 'FontSize', 14); xlabel('$t$ [$m$s]', 'FontSize', 14);
        title('Energy', 'FontSize', 14);
        grid on
    
        tVal = Observ.tVecPlot(end)*1000/w0;
        sgtitle(sprintf('$\\mu =%.3f \\hbar\\omega_0$, $t=%.1f$ms',muchem,tVal),'Interpreter', 'latex','FontSize', 14)
        drawnow
        saveas(gcf,sprintf('./Results/Data_3D/RealTimeDynamics/Run_%03i/Figures/Image_%i.jpg',run_index,ii))
        img = imread(sprintf('./Results/Data_3D/RealTimeDynamics/Run_%03i/Figures/Image_%i.jpg',run_index,ii));
        writeVideo(outputVideo,img)
        clf
    end
    
    close(outputVideo)
    close(figure(1))
    delete(sprintf(horzcat(folder_path, '/Run_%03i/Figures/*.jpg'),run_index)) % deleting images after movie is made
end