How to encrypt image using my own chaotic map instead of baker's map?

Professor, i am using a choatic map (in code it is with name a squence, this sequence is gotten from a matrix S2, that S2 (minima of a chaotic system, i have already solved that chatic system no need to solve it again) ) for image encrotion by using algorithms given in this paper (Image encryption using chaos and block cipher, with DOI No10.5539/cis.v4n1p172). But code is not running corecctly. Please checkt it. I have attached the all files. It is an image encryption scheme based on combination of pixel shuffling and new modified version of simplified AES. Chaotic map is used for shuffling and improving S-AES efficiency through S-box design. Chaos is used to expand diffusion and confusion in the image.

'' in the paper he use a baker's map "

(FOLLOWIGN IS THE CODE FOR MY MAP i,e SEQUENCE instead OF BAKER's MAP)

%% (1)
% Function for Shuffling Pixels Using Chaotic Sequence (Replacing Baker's Map by my sequence)
function shuffled_img = shuffle_pixels_with_sequence(img, sequence)
[rows, cols] = size(img);
num_points = size(sequence, 1);
% Generate initial grid coordinates
[X, Y] = meshgrid(1:cols, 1:rows);
X = X(:) / cols;
Y = Y(:) / rows;
% Iterate to compute the shuffled positions using the chaotic sequence
for i = 1:num_points
    X_next = mod(X + sequence(i, 1), 1);
    Y_next = mod(Y + sequence(i, 2), 1);
    X = X_next;
    Y = Y_next;
end
% Ensure X and Y values are within the valid range
X = ceil(X * cols);
Y = ceil(Y * rows);
% Make sure X and Y are within the range [1, rows] and [1, cols]
X = max(min(X, cols), 1);
Y = max(min(Y, rows), 1);
% Rearrange pixels based on shuffled indices
indices = sub2ind([rows, cols], Y, X);
shuffled_img = img(indices);
shuffled_img = reshape(shuffled_img, rows, cols);
% **Step 1: Vertical Permutation**
Pmap = randperm(rows);
shuffled_img = shuffled_img(Pmap, :);
% **Step 2: Horizontal Permutation**
Pmap = randperm(cols);
shuffled_img = shuffled_img(:, Pmap);
% **Step 3: Apply Circular Shift for Diagonal Permutation**
for it = 2:rows
    Shiftsize = mod(num_points - it + 1, cols);  % Prevents exceeding bounds
    shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, Shiftsize]);
end
% **Step 4: Apply Vertical Permutation Again**
shuffled_img = shuffled_img(Pmap, :);
% **Step 5: Apply Final Diagonal Permutation**
for it = 2:rows
    Shiftsize = mod(num_points - it + 1, cols);
    shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, -Shiftsize]);
end
end
%% (2)
function encrypted_img = s_aes_encrypt(img, sbox)
[rows, cols] = size(img);
% Ensure the S-box is expanded to match the image size
if size(sbox, 1) < rows || size(sbox, 2) < cols
    sbox_expanded = repmat(sbox, ceil(rows / 256), ceil(cols / 256));
    sbox_expanded = sbox_expanded(1:rows, 1:cols);
else
    sbox_expanded = sbox(1:rows, 1:cols);
end
% Convert image to double for computations
img = double(img);
% Round 1
img = add_round_key(img, sbox_expanded); % Add Round Key
img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles
img = shift_rows(img); % Shift Rows
img = mix_columns(img); % Mix Columns
img = add_round_key(img, sbox_expanded); % Add Round Key
% Round 2
img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles
img = shift_rows(img); % Shift Rows
img = add_round_key(img, sbox_expanded); % Add Round Key
% Convert back to uint8 for image display
encrypted_img = uint8(img);
end
%% (3)
function state = substitute_nibbles(state, sbox)
[rows, cols] = size(state);

% Substitute each nibble using the S-box
for i = 1:rows
    for j = 1:cols
        % Ensure the value is within the range of the S-box
        nibble_value = mod(state(i, j), 256) + 1; % Map to 1-256
        state(i, j) = sbox(nibble_value);
    end
end
end
%% (4)
function state = shift_rows(state)
% Shift rows as per S-AES
state(2, :) = circshift(state(2, :), -1); % Shift second row by 1 position
end
%% (5)
function state = mix_columns(state)
% Mix columns as per S-AES
temp = state;
state(1, 1) = bitxor(temp(1, 1), temp(2, 1));
state(1, 2) = bitxor(temp(1, 2), temp(2, 2));
state(2, 1) = bitxor(temp(1, 1), temp(2, 1));
state(2, 2) = bitxor(temp(1, 2), temp(2, 2));
end
%% (6)
function state = add_round_key(state, round_key)
% XOR the state with the round key
state = bitxor(state, round_key);
end
%% (7)
%% (5) generate_sbox - Fully Corrected Version
function sbox = generate_sbox(sequence)
NoIt = length(sequence);  
D = zeros(1, NoIt);
% **Step 1-5: Compute D(it) Using Iterative Transformation**
for it = 1:NoIt
    idx = mod(it - 1, length(sequence)) + 1;
    X = sequence(idx, 1);
    Y = sequence(idx, 2);
    % ? Corrected: Use iterative transformations
    X_next = mod(X + Y, 1);
    Y_next = mod(X + 2 * Y, 1);
    
    % Update X, Y for next iteration
    X = X_next;
    Y = Y_next;
    % ? Compute D(it) using correct transformation
    D(it) = sqrt(abs(X.^3 + Y.^3));
end
% **Step 6-10: Apply Sorting & Nonlinear Transformation**
D = sort(D, 'ascend');
M = max(D);
for it1 = 1:NoIt
    for it2 = 1:NoIt
        if D(it2) == M
            D(it2) = -16 + it1; 
        end
    end
end
% **Step 11-20: Normalize & Create the S-Box**
D = abs(D);  
max_D = max(D);
sbox = uint8(mod(round(D * 255 / max_D), 256));
% **Ensure the S-Box is exactly (256 × 256)**
sbox = repmat(sbox(:), ceil(256*256/numel(sbox)), 1);
sbox = reshape(sbox(1:256*256), 256, 256);
end
%%%  (8)
%% (8) encryption_main_subscript - Fully Matches Figure 5
clc; clear; close all;
%% **Step 1: Load and Preprocess Image**
filePath = 'C:\Users\HP EliteBook 840 G3\Desktop\pic1.jpeg';
img = imread(filePath);

Error using imread>get_full_filename (line 633)
Unable to find file "C:\Users\HP EliteBook 840 G3\Desktop\pic1.jpeg".

Error in imread (line 395)
fullname = get_full_filename(filename);

% Convert to grayscale if the image is RGB
if size(img, 3) == 3
    img = rgb2gray(img);
end
% Resize image to 256x256 (if not already)
img = imresize(img, [256, 256]);
% Display (a) Original Image
figure;
subplot(2, 3, 1);
imshow(img);
title('(a) Original Image');
%% **Step 2: Load Sequence from S2.xlsx for Pixel Shuffling**
sequenceFile = 'C:\Users\HP EliteBook 840 G3\Desktop\S2.xlsx';
data = xlsread(sequenceFile);
data = data(:); % Flatten matrix
% Scale the sequence between 0 and 1 (Normalization)
scaled_sequence = (data - min(data)) / (max(data) - min(data));
maxima_x = scaled_sequence(1:end-1);
maxima_y = scaled_sequence(2:end);
sequence = [maxima_x, maxima_y];
%% **Step 3: Apply Algorithm 1 (Pixel Shuffling)**
shuffledImg = shuffle_pixels_with_sequence(img, sequence);
% Display (b) Shuffled Image
subplot(2, 3, 2);
imshow(shuffledImg);
title('(b) Shuffled Image');
%% **Step 4: Apply Algorithm 2 (Generate Key-Dependent S-Box)**
sbox = generate_sbox(sequence); % Generate dynamic S-Box
%% **Step 5: Apply Simplified AES (S-AES) Encryption**
I = s_aes_encrypt(shuffledImg, sbox); % Pass S-Box explicitly
% Display (c) Ciphered Image
subplot(2, 3, 3);
imshow(I);
title('(c) Encrypted Image');
%% **Step 6: Compute and Display Histograms**
subplot(2, 3, 4); imhist(img); title('(d) Histogram of Original Image');
subplot(2, 3, 5); imhist(shuffledImg); title('(e) Histogram of Shuffled Image');
subplot(2, 3, 6); imhist(I); title('(f) Histogram of Encrypted Image');
% %% **Step 7: Save Encrypted Image**
% outputPath = 'C:\Users\HP EliteBook 840 G3\Desktop\encrypted_image.png';
% imwrite(cipheredImg, outputPath);
% disp(['Encrypted image saved to: ', outputPath]);

For S1; from the following code just consider MAXIMA and in above codes i have converted S1 in sequence.

The system this paper (DOI: 10.1017/S0022112003003835) Thank you.
%% funn
function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);
%persistent iter
%if isempty(iter)
%  iter = 0;
%end
%iter = iter + 1;
%if mod(iter,1000)==0
%  t
%  iter = 0;
%end
G = y(1:nx);
F = y(nx+1:2*nx);
% F
% Compute spatial derivatives
dFdx = zeros(nx,1);
d2Fdx2 = zeros(nx,1);
dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);
dFdx(nx) = Hdot(t);           % This corresponds to F'(1,t) = H'(t)
d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;
%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;
Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);
d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;
% Compute temporal derivatives
dFdt = zeros(nx,1);
dFdt(1) = F(1);
dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...
    F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);
dFdt(nx) = F(nx) + Hdot(t);     % This corresponds to F(1,t) =- H'(t)
% G
% Compute spatial derivatives
dGdx = zeros(nx,1);
d2Gdx2 = zeros(nx,1);
dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);
d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;
% Compute temporal derivatives
dGdt = zeros(nx,1);
dGdt(1) = G(1);
dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...
    1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...
    1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...
    2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...
    (d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...
    (H(t)^2*Reynolds);
dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);
%Taken from the article, should be the same as set
%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));
dydt = [dGdt;dFdt];
end
%%%   MAIN SUBSCRIPTT
format long;
clear all;
close all;
clc;
% Parameters
delta = 0.3;       % Amplitude of H
Reynolds = 900;   % Reynolds number
tstart = 0;        % Start time
tend = 206400;     % End time
nx = 101;          % Number of spatial points
xstart = 0.0;      % Spatial domain start
xend = 1.0;        % Spatial domain end
x = linspace(xstart, xend, nx).'; % Spatial grid
dx = (xend - xstart) / (nx - 1);  % Spatial step size
tspan = [tstart, tend];           % Time span
% numpoints = 1000;
% tspan = linspace(tstart, tend,numpoints);           % Time span
% Initial conditions
G0 = zeros(nx, 1); % Initial condition for G
F0 = zeros(nx, 1); % Initial condition for F
y0 = [G0; F0];     % Combine into one vector
% Define H(t) and its derivative
H = @(t) 1 + delta * cos(2 * t);
Hdot = @(t) -2 * delta * sin(2 * t);
% Mass matrix for DAE
M = [eye(nx), zeros(nx); zeros(nx, 2*nx)];
M(1, 1) = 0; % Boundary condition G(0,t) = 0
M(nx, nx) = 0; % Boundary condition G(1,t) = 0
options = odeset('Mass', M, 'RelTol', 1e-3, 'AbsTol', 1e-6);
% Solve the system
[T, Y] = ode23t(@(t, y) funn(t, y, x, nx, dx, delta, Reynolds, H, Hdot), tspan, y0, options);
% Extract G and F
G = Y(:, 1:nx);
F = Y(:, nx+1:2*nx);
% Evaluate G at eta = 1 using interpolation
eta_query = 1; % Change from 1/4 to 1
G_eta_1 = zeros(size(T)); % Preallocate
for i = 1:length(T)
    G_eta_1(i) = interp1(x, G(i, :), eta_query, 'spline'); % Interpolation
end
% Plot G(1, t) over time
figure(1);
plot(T, G_eta_1, 'LineWidth', 1.5);
xlabel('Time (t)');
ylabel('G(1, t)');
title('Time evolution of G(1, t)');
grid on;
[maxima_G_eta_1, locs_max_G_eta_1] = findpeaks(G_eta_1); % Maxima of G(1/4, t)
[minima_G_eta_1, locs_min_G_eta_1] = findpeaks(-G_eta_1); % Minima of G(1/4, t)
minima_G_eta_1 = -minima_G_eta_1; % Correct sign of minima
% Display counts of maxima and minima
sajjad1 = length(maxima_G_eta_1);
sajjad2 = length(minima_G_eta_1);
disp(['Number of maxima: ', num2str(sajjad1)]);
disp(['Number of minima: ', num2str(sajjad2)]);
% Multiply maxima and minima by 10^15 and apply mod 256
scaled_maxima = mod(round(maxima_G_eta_1 * 1e15), 256);
scaled_minima = mod(round(minima_G_eta_1 * 1e15), 256);
% Display scaled maxima and minima
disp('Scaled maxima (mod 256):');
disp(scaled_maxima);
disp('Scaled minima (mod 256):');
disp(scaled_minima);
% Return Map for Maxima
figure(2); % New figure
plot(maxima_G_eta_1(1:end-1), maxima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'r', 'MarkerSize', 4);
xlabel('Maxima_{n}');
ylabel('Maxima_{n+1}');
title('Return Map of Maxima');
grid on;
% Return Map for Minima
figure(3); % Another new figure
plot(minima_G_eta_1(1:end-1), minima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'g', 'MarkerSize', 4);
xlabel('Minima_{n}');
ylabel('Minima_{n+1}');
title('Return Map of Minima');
grid on;

Please feel to ask any other information about it.

Thank you.