function plot2DTE(a, b, M, phi,T,eps)
    %%%%%%%%%%%%%save data
    filename = sprintf('T%.2feps%.2f.mat',T,eps);  
    filename = strrep(filename, '.', 'p');
    save(filename, 'phi', 'T', 'a', 'b', 'M', 'eps');
    fprintf('The data has been saved as: %s\n', filename);
    %%%%%%%%%%%%%%%%% Create the "figure" folder (if it does not exist)
    script_dir = fileparts(mfilename('fullpath'));
    figure_dir = fullfile(script_dir, 'outputimages');
    if ~exist(figure_dir, 'dir')
        mkdir(figure_dir);
    end
    % plot2D - Plot the modulus square of a two-dimensional complex wave function
    % Input:
    %   a, b   - Endpoints of the spatial interval
    %   M      - Number of grids
    %   phi    - Complex wave function (2 × M × M dimension)

    % Spatial grid
    h = (b - a) / M;
    t = a : h : b - h;
    [x, y] = meshgrid(t);

    % Calculate the probability density of the two components
    rho = abs(squeeze(phi(1, :, :))).^2+abs(squeeze(phi(2, :, :))).^2;

    % Set the graphics window
    figure('Position', [100, 100, 600, 600]);  

    % Public Chart Settings
    cmap = 'turbo';
    axlims = [a, b, a, b];  
    mesh(x, y, rho);
    view(2);
    title(sprintf('$T = %d\\ , \\ \\varepsilon = %.2f$', T, eps), ... 
      'Interpreter', 'latex', 'FontSize', 14);
    colormap(cmap);
    colorbar;
    axis equal tight; 
    xlabel('$x$', 'Interpreter', 'latex', 'FontSize', 18);
    ylabel('$y$', 'Interpreter', 'latex', 'FontSize', 18);
    set(gca, 'FontSize', 18);
    %%%%Save to a subfolder in the current directory
    exportgraphics(gcf, fullfile('outputimages', sprintf('T=%d_eps=%.2f.png', T, eps)), ...
               'Resolution', 300);
              
end