% Credits are indicated where needed. The general idea is based on a template by Vel (vel@LaTeXTemplates.com) and Frits Wenneker.

\documentclass[11pt, a4paper]{article} % General settings in the beginning (defines the document class of your paper)
% 11pt = is the font size
% A4 is the paper size
% “article” is your document class
%----------------------------------------------------------------------------------------
%	Packages
%----------------------------------------------------------------------------------------
% Necessary
\usepackage{tabularx}
\usepackage{float}
\usepackage[german,english]{babel} % English and German language 
\usepackage{booktabs} % Horizontal rules in tables 
% For generating tables, use “LaTeX” online generator (https://www.tablesgenerator.com)
\usepackage{comment} % Necessary to comment several paragraphs at once
\usepackage[utf8]{inputenc} % Required for international characters
\usepackage[T1]{fontenc} % Required for output font encoding for international characters

% Might be helpful
\usepackage{amsmath,amsfonts,amsthm} % Math packages which might be useful for equations
\usepackage{tikz} % For tikz figures (to draw arrow diagrams, see a guide how to use them)
\usepackage{tikz-cd}
\usetikzlibrary{positioning,arrows} % Adding libraries for arrows
\usetikzlibrary{decorations.pathreplacing} % Adding libraries for decorations and paths
\usepackage{tikzsymbols} % For amazing symbols ;) https://mirror.hmc.edu/ctan/graphics/pgf/contrib/tikzsymbols/tikzsymbols.
\usepackage{blindtext} % To add some blind text in your paper
\DeclareFontFamily{T1}{calligra}{}
\DeclareFontShape{T1}{calligra}{m}{n}{<->s*[1.44]callig15}{}
\DeclareMathAlphabet\mathcalligra   {T1}{calligra} {m} {n}
\DeclareMathAlphabet\mathzapf       {T1}{pzc} {mb} {it}
\usepackage[backref]{hyperref}
\usepackage{cite}
\usepackage[backend=biber,style=ieee,sorting=none]{biblatex} % Use biber as the backend
\addbibresource{References/Ref.bib} % Name of your .bib file
%\include{biblfe.bib}
\hypersetup{
    colorlinks=true,   % Enable colored text for links (no boxes)
    linkcolor=black,   % Links to equations
    citecolor=black,   % Links to references
    urlcolor=black     % Links to URLs
}
%\urlstyle{same}

\usepackage{graphicx}
%\usepackage{subfigure}
\usepackage{caption}
\usepackage{subcaption}

%\usepackage{authblk}
%---------------------------------------------------------------------------------
% Additional settings
%---------------------------------------------------------------------------------

%---------------------------------------------------------------------------------
\usepackage{titlesec}
\titleformat{\subsection}[block]
  {\normalfont\itshape}  % Format: normal font and italics (you can change this)
  {\thesubsection}{1em}{}
%\renewcommand{\subsection}[1]{\noindent\textbf{\normalfont\itshape #1}}

% Define your margins
\usepackage{geometry} % Necessary package for defining margins

\geometry{
	top=2cm, % Defines top margin
	bottom=2cm, % Defines bottom margin
	left=2.2cm, % Defines left margin
	right=2.2cm, % Defines right margin
	includehead, % Includes space for a header
	%includefoot, % Includes space for a footer
	%showframe, % Uncomment if you want to show how it looks on the page 
}

\setlength{\parindent}{0pt} % Adjust to set you indent globally 

%---------------------------------------------------------------------------------
% Define your spacing
\usepackage{setspace} % Required for spacing
% Two options:
\linespread{1.0}
%\onehalfspacing % one-half-spacing linespread

%---------------------------------------------------------------------------------------------------------------------------------------------------
% Define your fonts
\usepackage[T1]{fontenc} % Output font encoding for international characters
\usepackage[utf8]{inputenc} % Required for inputting international characters

\usepackage{XCharter} % Use the XCharter font


%---------------------------------------------------------------------------------
% Define your headers and footers

\usepackage{fancyhdr} % Package is needed to define header and footer
\pagestyle{fancy} % Allows you to customize the headers and footers
\usepackage{epsfig}
\usepackage{epstopdf}
%\renewcommand{\sectionmark}[1]{\markboth{#1}{}} % Removes the section number from the header when \leftmark is used

% Headers
\lhead{} % Define left header
\chead{\textit{}} % Define center header - e.g. add your paper title
\rhead{} % Define right header

% Footers
\lfoot{} % Define left footer
\cfoot{\footnotesize \thepage} % Define center footer
\rfoot{ } % Define right footer

%---------------------------------------------------------------------------------
%	Add information on bibliography
%\usepackage{natbib} % Use natbib for citing
%\usepackage{har2nat} % Allows to use harvard package with natbib https://mirror.reismil.ch/CTAN/macros/latex/contrib/har2nat/har2nat.pdf

% For citing with natbib, you may want to use this reference sheet: 
% http://merkel.texture.rocks/Latex/natbib.php

%---------------------------------------------------------------------------------
% Add field for signature (Reference: https://tex.stackexchange.com/questions/35942/how-to-create-a-signature-date-page)
\newcommand{\signature}[2][5cm]{%
  \begin{tabular}{@{}p{#1}@{}}
    #2 \\[2\normalbaselineskip] \hrule \\[0pt]
    {\small \textit{Signature}} \\[2\normalbaselineskip] \hrule \\[0pt]
    {\small \textit{Place, Date}}
  \end{tabular}
}

\renewcommand{\figureautorefname}{Fig.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                             %  AUTHORS BIO  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{ \textbf{Bifurcation Analysis and Soliton Structures of Davey-Stewartson Fokas System} }
\author{
Muhammad Hammad$^1$\footnote{hammad.22311039@math.qau.edu.pk }, 
Amjad Hussain$^1$\footnote{a.hussain@qau.edu.pk},
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %  BEGNINING  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\maketitle
% Print your title, author name and date;
\setcounter{page}{1} % Sets counter of page to 1
\begin{center}
$^1$Department of Mathematics, Quaid-i-Azam University 45320, Islamabad 44000,  Pakistan.
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %  ABSTRACT  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}
In this research, we studied the (2+1)-dimensional Davey-Stewartson Fokas (DS-Fokas) system, which serves as an optimal model for nonlinear pulse propagation in mono-mode optical fibers. We employ the Jacobi elliptic function approach to obtain the novel soliton solutions for the DS-Fokas system. The employed method is a very efficient and robust mathematical approach for solving non-linear models of various nonlinear Schrödinger’s equations (NLSEs) in mathematical physics and sciences. The obtained solutions are useful and significant in elucidating the DS-Fokas system's physical aspects, as they provide insights. Furthermore, we discuss these obtained solutions graphically using 3D and 2D graphs to gain a deep understanding and vision of the analytical results. We also looked at the unpredictable and changing behaviors of the system we studied by using phase portraits, quasi-periodic and chaotic portraits, Poincare maps, bifurcation diagrams, and sensitivity. The theory of planar dynamical systems looks at chaotic patterns in the systems under study when the disturbance term $\cos \omega t$ is added. Numerical simulations demonstrate how changes in frequency and amplitude impact the dynamics of the system.


\noindent\textbf{Keywords:} Davey-Stewartson Fokas System; Dynamical analysis;  Chaotic behavior; Bifurcation; Sensitivity analysis; Soliton Solutions.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                             %  INTRODUCTION  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}\label{S0}
Nonlinear evolution equations are frequently employed in several scientific domains, including fluid dynamics, marine science, medicinal sciences, optics, and plasma physics \autocite{Ref1,Refn1,Refr12,Ref2,Refr2,Ref3,Ref4,Ref5,mohan2024generalization}.  Researchers are increasingly interested in solving nonlinear evolution problems. To find exact solutions for solitary waves, different mathematical methods have been created. These include the tanh method \autocite{Ref6}, the Lie Symmetry Method \autocite{Refn} and Unified Approach \autocite{kumar2024exploring}, the Hirota bilinear method \autocite{Ref8}, the direct symbolic approach
 \autocite{kumar2024bilinearization}, the generalized exponential rational function method \autocite{Ref10}, the improved $\frac{G'}{G}$-expansion method \autocite{Ref11}, the extended Jacobi elliptic expansion \autocite{Ref12}, and the modified Sardar-sub equation method \autocite{Ref13}.
Besides, the dynamical properties of nonlinear planar dynamical systems have been the subject of numerous investigations. For example, Rani et al. have analyzed the dynamical aspects of the Complex Kraenkel–Manna–Merle System \autocite{rani2024dynamical}, while Kumar has conducted a qualitative analysis of integrable Kuralay equations and a bifurcation analysis of the (2+1)-dimensional Konopelchenko–Dubrovsky (KD) system \autocite{kumar2024dynamic,kumar2023dynamical}.

The main goal of this study is to look into the nonlinear (2+1)-dimensional Davey-Stewartson-Fokas system. The DS-Fokas system plays a crucial role in nonlinear optics, plasma physics, and fluid dynamics. Fiber-optic communication and optical signal processing also utilize this system. The DS-Fokas system is an extension of the nonlinear Schrödinger equation, derived by Fokas in 1994 (\autocite{Ref31}). The system is a complicated integrable model that takes the NLSE to higher dimensions. It explains the evolution of complex wave amplitudes and their corresponding mean flow fields in two spatial dimensions and over time. Equation \eqref{E1} represents the mathematical structure of this system, as given in \autocite{Ref14}:
\begin{equation}\label{E1}
   \begin{align*}
       \iota\ \phi_\tau + k_1\ \phi_{\mu \nu} + k_2\ \phi\ \psi &= 0, \\
        k_3\ \psi_\nu - k_4\ (\left |\phi \right |^2)_\mu &= 0.
   \end{align*} 
\end{equation}
In single-mode optical fibers, nonlinear pulse propagation is shown by the real function $\phi = \phi(\mu, \nu, \tau)$ and the imaginary function $\psi = \psi(\mu, \nu, \tau)$. The constants $k_j, \quad (j=1,2,3,4)$, are arbitrary to be assessed using the proposed computational framework.\\

The DS-Fokas system is an important model in optics that explains how light pulses move through single-mode optical fibers in a nonlinear way. Analytical study of this system yields novel perspectives on modern nonlinear science, which makes its importance and the basis for our research choice clear. Our focus is utilizing the JEF technique and detecting chaos in the DS-Fokas system. The Jacobian elliptic function method works well for the DS-Fokas system because it accurately describes complex wave patterns that repeat over time. This method helps to analyze equations that involve nonlinearity and turbulence, as seen in fluid and plasma physics. It gives answers as elliptic functions that connect spherical and hyperbolic functions, showing how remote wave solutions change into soliton solutions. To find soliton solutions, researchers have studied the DS-Fokas system using various analytical and numerical methods, like the modified auxiliary equation method and the extended $\frac{G′}{G^2}$-expansion method \autocite{Ref16}. Some authors employ the ansatz method and the projective Riccati expansion method to identify solitary wave solutions for the DS-Fokas system. The Hirota bilinear expansion method has also been used \autocite{Ref15} to find new analytical solutions for the same system. The Jacobi elliptic function method used in this study is superior to other methods currently used. This method provides a wide variety of soliton solutions and more general and flexible results that are suitable for the DS-Fokas system.

The uniqueness and significance of the obtained solutions are the main point of this research. We have discovered new solutions for the DS-Fokas system, previously unreported, using the Jacobi elliptic function method. In addition, we uncover the dynamical features of the same system using different tools. This work is a significant addition to the field of finding analytical solutions and dynamical features for complex systems like the DS-Fokas system and may be applicable to other nonlinear models in the future.
\section{Dynamical System of the Govening Model}\label{S3}
Let us consider the following wave transformation:
\begin{equation}\label{E2}
   \begin{align*}
      \phi(\mu,\nu,\tau) = e^{\iota \delta}\ \Phi(\epsilon), \\
      \psi(\mu,\nu,\tau) = \Psi(\epsilon),
   \end{align*} 
\end{equation}
with
\begin{equation}\label{E3}
    \delta = \lambda_{1}\ \mu + \lambda_{2}\ \nu + \lambda_{3}\ \tau, \quad \epsilon = \alpha\ \mu + \beta\ \nu - \gamma\ \tau.
\end{equation}
where $\lambda_{i}, i=1,2,3$, and $\alpha, \beta, \gamma$ are nonzero arbitrary constants. Using Eq. \eqref{E3} alongside Eq. \eqref{E2} in Eqs. \eqref{E1} we obtain the following system having the real part as:
\begin{equation}\label{E4}
   k_{{1}} {
\alpha}^{2}
\Phi'' + k_{{2}}\Phi\, \Psi -(\lambda_{{3}} + k_{{1}}{\lambda_{{1}}}^{2})\Phi
=0,
\end{equation}
\begin{equation}\label{E5}
    -\gamma\ k_{3}\ \Psi'' - 2 \alpha\ k_{4} \Phi\ \Phi'=0,
\end{equation}
and the imaginary part as:
\begin{equation}\label{E6}
    (\gamma - 2 \alpha\ k_{1}\ \lambda_{1})\ \Phi'=0.
\end{equation}
The imaginary part will vanish for $\gamma = 2 \alpha\ k_{1}\ \lambda_{1}$. Integrating Eq. \eqref{E5} gives
\begin{equation}\label{E7}
    \Psi = - \frac{k_{4}\ \alpha}{k_{3}\ \gamma}\ \Phi^2.
\end{equation}
Substituting Eq. \eqref{E7} in Eq. \eqref{E5}, we get
\begin{equation}\label{E8}
        k_{{1}}\ {\alpha}^{2} \Phi''
- \left( k_{{1}}{\lambda_{{1}}}^{2}+\lambda_{{3}} \right) \Phi-{\frac 
{k_{{2}}k_{{4}}\alpha}{k_{{3}}\gamma}}{\Phi}^{3}=0,
\end{equation}
\begin{equation}\label{E9}
        \Phi'' = \eta_{1} \Phi + \eta_{2} {\Phi}^{3}.
\end{equation}
where $\eta_{1}=\frac{\left( k_{{1}}{\lambda_{{1}}}^{2}+\lambda_{{3}} \right)}{k_{{1}}\ {\alpha}^{2}}$ and $\eta_2= {\frac {k_{{2}}k_{{4}}}{k_{1} k_{{3}}\alpha\gamma}}$ and $\left('\right)$ denotes $\frac{d}{d\epsilon}$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                          %  PHASE PORTRAITS  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Phase Portraits}\label{SS41}
A two-dimensional representation of phase space, a phase portrait, shows the interaction between state variables. Each point in the portrait corresponds to a fixed-point solution, while closed curves indicate periodic solutions. Chaotic solutions are represented by discrete, non-repeating trajectories. Fig. \ref{fig:subfig1a}-d provide examples of such portraits for the analyzed model.
This section will explore Eq. \eqref{E1} utilizing bifurcation theory and phase portrait analysis \autocite{Ref17,Ref18}. Changing parameters can cause qualitative changes in a dynamical system with a bifurcation as defined by Eq. \eqref{E1}. Bifurcation theory examines family-level bifurcations. This is performed by recognizing common bifurcation patterns. Using the Galilean transformation, the planar dynamical system for Eq. \eqref{E9} may be expressed as:
\begin{equation}\label{E10}
\left\{
    \begin{align*}
        \Phi' &= \chi, \\
        \chi' &= \eta_{1} \Phi + \eta_{2} {\Phi}^{3}.
    \end{align*}.
\right.
\end{equation} 
Equilibrium points are obtained by solving the following system:
\begin{equation}\label{E11}
\left\{
    \begin{align*}
         \chi = 0,\\
         \eta_{1} \Phi + \eta_{2} {\Phi}^{3} = 0.
    \end{align*}
\right.
\end{equation}
As a result, three equilibrium points are as under:
\begin{equation*}
    E_{1} = (0,0), \quad E_{2} = \left(\sqrt{-\frac{\eta_{1}}{\eta_{2}}} , 0 \right), \quad E_{3} = \left(-\sqrt{-\frac{\eta_{1}}{\eta_{2}} }, 0 \right).
\end{equation*}
    The determinant of the Jacobian matrix of transformation has the following form:
    \begin{equation}\label{E12}
        J(\Phi,\chi) = \begin{vmatrix}
        0 & 1 \\
        \eta_{1} + 3 \eta_2 \Phi^2 & 0
        \end{vmatrix} = -\left(\eta_{1} + 3 \eta_2 \Phi^2\right).
    \end{equation}

\textbf{Case 1:} \(\eta_{1} > 0, \quad \eta_{2} > 0.\)\\
Eq. \eqref{E11} produces single equilibrium point, \(E_1 = (0,0)\) which may be attained by using some positive values of \(\eta_{1}\) and \(\eta_{2}\). This is shown in Fig. \ref{fig:subfig1a}. It can be seen that \(E_1\) displays as a saddle point.

\textbf{Case 2:} \(\eta_{1} > 0, \quad \eta_{2} < 0\)

Eq. \eqref{E11} yields three equilibrium points, \(E_1 = (0,0)\) and \(E_{2,3} = \left(\pm\sqrt{-\frac{\eta_{1}}{\eta_{2}}} , 0 \right)\), which may be achieved by using the values of \(\eta_{1}\) and \(\eta_{2}\) which satisfy the conditions assigned in this case. This is shown in Fig. \ref{fig:subfig1b}. As can be seen, \(E_1\) is a saddle point and \(E_{2,3}\) displays as the centre point.

\textbf{Case 3:} \(\eta_{1} < 0, \quad \eta_{2} > 0\)

Eq. \eqref{E11} gives three equilibrium points, \(E_1 = (0,0)\) and \(E_{2,3} = \left(\pm\sqrt{-\frac{\eta_{1}}{\eta_{2}}} , 0 \right)\), which may be achieved by taking some negative values of \(\eta_{1}\) and positive values for \(\eta_{2}\). This is shown in Fig. \ref{fig:subfig1c}. As seen, \(E_1\) is a center point and \(E_{2,3}\) displays the saddle point.

\textbf{Case 4:} \(\eta_{1} < 0, \quad \eta_{2} < 0\)

Eq. \eqref{E11} yields single equilibrium point, \(E_1 = (0,0)\)  which may be obtained by using negative values for both \(\eta_{1}\) and \(\eta_{2}\). This is shown in Fig. \ref{fig:subfig1d}. As seen, \(E_1\) is a center point.


\begin{figure}[H]
    \centering
    % First row
    \begin{subfigure}[b]{0.45\textwidth}
        \centering
        \includegraphics[[width=3in,height=3in]{fig1a.eps}
        \caption{$\eta_{1} = 1, \eta_{2} = 1$}
        \label{fig:subfig1a}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.45\textwidth}
        \centering
        \includegraphics[width=3in,height=3in]{fig1b.eps}
        \caption{$\eta_{1} = 1, \eta_{2} = -1$}
        \label{fig:subfig1b}
    \end{subfigure}
    % Second row
    \begin{subfigure}[b]{0.45\textwidth}
        \centering
        \includegraphics[[width=3.1in,height=3.1in]{fig1c.eps}
        \caption{$\eta_{1} = -1, \eta_{2} = 1$}
        \label{fig:subfig1c}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.45\textwidth}
        \centering
        \includegraphics[[width=3.1in,height=3.1in]{fig1d.eps}
        \caption{$\eta_{1} = -1, \eta_{2} = -1$}
        \label{fig:subfig1d}
    \end{subfigure}
    \caption{Comparison of different cases of $\eta_{1}$ and $\eta_{2}$.}
    \label{fig:main}
\end{figure}
\subsection{Quasi-Periodic and Chaotic Behavior}\label{SubS2}
This subsection will examine the dynamical system \eqref{E13} by introducing a perturbation term $\cos{\omega t}$ in the system \eqref{E10}. The modified system \eqref{E13} will be utilized to detect the quasi-periodic and chaotic patterns:
\begin{equation}\label{E13}
\left\{
    \begin{align*}
        \Phi' &= \chi, \\
        \chi' &= \eta_{1} \Phi + \eta_{2} {\Phi}^{3} + A\cos{\omega t}.
    \end{align*}
\right.
\end{equation}
Periodic disruptions are rather prevalent in real-world systems and can cause havoc, occasionally leading to chaotic consequences; for more detailed insights, readers are suggested these articles \autocite{Ref25,Ref26,Ref27}. Our research examines how the systemstem \eqref{E13} responds to variations in the intensity and frequency of these disruptions. We set the major parameters, $\eta_{1}=1, \eta_{2}=-1$, and tested the periodic, quasi-periodic, and chaotic dynamics at several levels of disruption. When we adjust the amplitude A to zero, system \eqref{E10} appears as depicted in Fig \ref{fig:QP1}. We experimented with the system's course by adjusting the frequency and amplitude, as shown in Figs. \eqref{fig:QP2} to \eqref{fig:QP4}. We used 2D and 3D phase space plots and time series graphs to figure out how system \eqref{E13} behaved in a periodic and quasi-periodic way.

    %%% Including Figures for Chaotic Behavior %%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.37\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.1a_A-0.000_W-0.200.png}
        \caption{3D Portrait}
         \label{fig:QP1a}
    \end{subfigure}%
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.1b_A-0.000_W-0.200.png}  
        \caption{2D Portrait}
         \label{fig:QP1b}
    \end{subfigure}%
    \\  
    \vspace{0.25cm}% Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.1c_A-0.000_W-0.200.eps} 
        \caption{Time series Plot for $Q$}
        \label{fig:QP1c}
    \end{subfigure}
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.1d_A-0.000_W-0.200.eps}  
        \caption{Time series Plot for $R$}
         \label{fig:QP1d}
    \end{subfigure}%
    \caption{Phase Portraits showing Periodic Behavior for Dynamical System \eqref{E13}, with $A=0,w=0.2$}
    \label{fig:QP1}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.37\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.8a_A-0.090_W-1.700.png}
        \caption{3D Portrait}
         \label{fig:QP2a}
    \end{subfigure}%
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.8b_A-0.090_W-1.700.eps}  
        \caption{2D Portrait}
         \label{fig:QP2b}
    \end{subfigure}%
    \\ 
    \vspace{0.25cm}% Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.8c_A-0.090_W-1.700.eps} 
        \caption{Time series Plot for $Q$}
        \label{fig:QP2c}
    \end{subfigure}
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.8d_A-0.090_W-1.700.eps}  
        \caption{Time series Plot for $R$}
         \label{fig:QP2d}
    \end{subfigure}%
    \caption{Phase Portraits showing Quasi-Periodic Behavior for Dynamical System \eqref{E13}, with $A=0.09,w=1.7$}
    \label{fig:QP2}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Third Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.37\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.6a_A-1.100_W-1.700.png}
        \caption{3D Portrait}
         \label{fig:QP3a}
    \end{subfigure}%
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.6b_A-1.100_W-1.700.png}  
        \caption{2D Portrait}
         \label{fig:QP3b}
    \end{subfigure}%
    
    \vspace{0.15cm}% Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.6c_A-1.100_W-1.700.eps} 
        \caption{Time series Plot for $Q$}
        \label{fig:QP3c}
    \end{subfigure}
    \hspace{0.28cm} % Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.6d_A-1.100_W-1.700.eps}  
        \caption{Time series Plot for $R$}
         \label{fig:QP3d}
    \end{subfigure}%
    \caption{Phase Portraits showing Chaotic Behavior for Dynamical System \eqref{E13}, with $A=1.1,w=1.7$}
    \label{fig:QP3}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fourth Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.37\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.3a_A-0.500_W-0.200.png}
        \caption{3D Portrait}
         \label{fig:QP4a}
    \end{subfigure}%
    \hspace{0.7cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.35\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.3b_A-0.500_W-0.200.eps}  
        \caption{2D Portrait}
         \label{fig:QP4b}
    \end{subfigure}%
    \vspace{0.15cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.3c_A-0.500_W-0.200.eps} 
        \caption{Time series Plot for $Q$}
        \label{fig:QP4c}
    \end{subfigure}
    \hspace{0.7cm} % Horizontal space between figures
    \begin{subfigure}[b]{0.38\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Fig2.3d_A-0.500_W-0.200.eps}  
        \caption{Time series Plot for $R$}
         \label{fig:QP4d}
    \end{subfigure}%
    \caption{Phase Portraits showing Chaotic Behavior for Dynamical System \eqref{E13}, with $A=0.5,w=0.2$}
    \label{fig:QP4}
\end{figure}

\subsection{Bifurcation}\label{BifSec}
Particularly when things turn chaotic, a bifurcation diagram is a useful tool for understanding system behavior \autocite{Ref19,Ref20}. It demonstrates how a system's behavior changes as you change a certain parameter. The graph's x-axis records parameter changes, while the y-axis depicts the system's various states. Chaos frequently manifests as chaotic, twisted patterns, such as branches or startling leaps. For our model (13), the diagram's rough, unpredictable patterns indicate chaotic behavior. Figure 7 illustrates this type of behavior with the values $\eta_{1}=1, \eta_{2}=-1, w=1.7$.
\begin{figure}[H]
    \centering 
    \includegraphics{bif.eps}
    \vspace{-4.7em}
    \caption{Bifurcation Diagram for Parameter $A$ vs $\Phi$ and $\chi$}
    \label{fig:fig2a}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Lyapunov Exponent}
Lyapunov exponents are an excellent tool to determine whether a system is behaving chaotically or not \autocite{Ref21,Ref22}. They help us comprehend how a few discrepancies in the initial configuration may develop or diminish over time. If the Lyapunov exponent is positive, it indicates pandemonium, since neighboring routes are rapidly diverging. Negative exponents represent stability, with routes moving closer together, but zero exponents indicate that objects remain around the same distance apart. Fig. \ref{fig: fig3} displays a positive Lyapunov exponent, indicating that our model exhibits chaotic behavior. We created this figure with the following parameter values: $\eta_{1}=1, \eta_{2}=-1, w=1.7, A=0.09$.
\begin{figure}[H]
    \centering 
    \includegraphics{Lyapfig.eps}
    
    \caption{Lyapunov Exponent Diagram}
    \label{fig: fig3}
\end{figure}
\subsection{Poincaré Maps}

The Poincare map method is a brilliant technique to investigate large systems by breaking them down into smaller defined phases. Instead of constantly observing the system, we capture snapshots of its phase portrait at regular intervals. This improves our understanding of stability and simplifies the analysis of complex systems \autocite{Ref23,Ref24}. In a Poincare map, a stable point indicates that the system is repeating itself (periodic behavior), whereas closed curves or scattered points indicate quasi-periodic or chaotic behavior. Chaos appears as a series of random dots, but quasi-periodic behavior appears as smooth loops. To make a Poincare map, we select a surface and track where the system's route passes it repeatedly. The Poincare map consists of these crossing locations. One may view instances of quasi-periodic and chaotic behavior in Figs. \ref{fig:PC1}-\ref{fig:PC4}.

%%% Including Figures for Poincare Map %%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   \begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
        \centering
        % First image without rectangle border
        \begin{tikzpicture}
            \node[anchor=south west, inner sep=0] at (0, 0)
            {\includegraphics[width=2.7in]{PoiCareMap0_A-0.000_w-0.20.eps}};
        \end{tikzpicture}
        \caption{for $\alpha = 1 $, $A = 0.009$ and $\omega = 1.1$}
        \label{fig:PC1}
    \end{subfigure}%
    \hspace{0.32cm} % Horizontal space between figures
    \begin{subfigure}[b]{0.4\textwidth}
        \centering
        % Second image without rectangle border
        \begin{tikzpicture}
            \node[anchor=south west, inner sep=0] at (0, 0)
            {\includegraphics[width=2.7in]{PoiCareMap1_A-0.090_w-1.70.eps}};
        \end{tikzpicture}
        \caption{for $\alpha = 1 $, $A = 0.05$ and $\omega = 1.1$}
        \label{fig:PC2}
    \end{subfigure}%
    \\  
    \vspace{0.3cm} % Vertical space between rows
    \begin{subfigure}[b]{0.4\textwidth}
        \centering
        % Third image without rectangle border
        \begin{tikzpicture}
            \node[anchor=south west, inner sep=0] at (0, 0)
            {\includegraphics[width=2.7in]{PoiCareMap2_A-1.100_w-1.70.eps}};
        \end{tikzpicture}
        \caption{for $\alpha = 1 $, $A = 0.1$ and $\omega = 1.1$}
        \label{fig:PC3}
    \end{subfigure}
    \hspace{0.3cm}  % Horizontal space between figures
    \begin{subfigure}[b]{0.4\textwidth}
        \centering
        % Fourth image without rectangle border
        \begin{tikzpicture}
            \node[anchor=south west, inner sep=0] at (0, 0)
            {\includegraphics[width=2.7in]{PoiCareMap3_A-0.500_w-0.20.eps}};
        \end{tikzpicture}
        \caption{for $\alpha = 1.8 $, $A = 0.001$ and $\omega = 1.1$}
        \label{fig:PC4}
    \end{subfigure}%
    \caption{Poincare Maps for Dynamical System \eqref{E13}}
    \label{fig:PCM}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 %  Quasi-Periodic and Chaotic Behavior  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






\begin{figure}[H]
    \centering
    \begin{subfigure}[b]{0.8\linewidth}
        \centering
        \includegraphics[width=\linewidth,height=5cm]{Fig3_A-0.000_w-0.20.png}
        \caption{Sensitivity Plot for System \eqref{E13} with $A=0$}
        \label{fig:SA1}
    \end{subfigure}
    
    \begin{subfigure}[b]{0.8\linewidth}
        \centering
        \includegraphics[width=\linewidth]{Fig2_(0.72, 0.1, 0)_A-0.100_w-1.70.png}
        \caption{Sensitivity Plot for System \eqref{E13} with $A=0.1$ and $\omega=1.7$}
        \label{fig:SA2}
    \end{subfigure}
    
    \begin{subfigure}[b]{0.8\linewidth}
        \centering
        \includegraphics[width=\linewidth]{Fig5_A-0.400_w-1.70.png}
        \caption{Sensitivity Plot for System \eqref{E13} with $A=0.4$ and $\omega=1.7$}
        \label{fig:SA3}
    \end{subfigure}
    
    \begin{subfigure}[b]{0.8\linewidth}
        \centering
        \includegraphics[width=\linewidth]{Fig2_(0.72, 0.1, 0)_A-0.500_w-1.70.png}
        \caption{Sensitivity Plot for System \eqref{E13} with $A=0.5$ and $\omega=1.7$}
        \label{fig:SA4}
    \end{subfigure}

   
\label{fig:sensitivity1}
   % \caption{Sensitivity plots for different values of $A$. Each plot shows the system's response as per Equation \eqref{E13}.}
    \label{fig:sensitivity}
\end{figure}

\begin{figure}
    \ContinuedFloat
    \centering
    \begin{subfigure}[b]{0.8\linewidth}
        \centering
        \includegraphics[width=\linewidth]{Fig2_A-0.500_w-0.20.png}
        \caption{Sensitivity Plot for System \eqref{E13} with $A=0.5$ and $\omega=0.2$}
        \label{fig:SA5}
    \end{subfigure}


    \caption{Sensitivity plots for different values of $A$. Each plot shows the system's response as per Equation \eqref{E13}.}
    \label{fig:sensitivity2}
\end{figure}

\subsection{Sensitivity Analysis}
It might not find chaotic patterns right away, but sensitivity analysis is useful for learning how systems respond to small changes in their initial conditions (IC) \cite{Ref28}. When it comes to chaotic behaviors, sensitivity analysis can explain what causes these crazy dynamics, particularly the legendary sensitivity to IC—where even minor variations may send things in very unexpected routes. Using sensitivity analysis, we may investigate how a system reacts to modest changes, as seen in Fig. \ref{fig:SA1}. For $A=0$, Fig. \ref{fig:SA1} depicts three periodic trajectories that begin from near ICs and follow similar routes. However, when $A$ and $\omega$ alter, things become interesting. Fig. \ref{fig:SA1} shows that increasing $A$ causes trajectories from the identical ICs to diverge. Instead, all three split onto separate, non-periodic pathways. This shift demonstrates how even little adjustments in settings may significantly affect the system's behavior. Figures \ref{fig:SA1} through \ref{fig:SA5} illustrate the evolution of dynamics as $A$ and $\omega$ alter. This procedure not only displays the system's complexity, but it also emphasizes the significance of thoroughly examining how initial conditions and parameters impact outcomes.

\section{Soliton Solution of the (2+1)-dimensional nonlinear coupled Davey-Stewartson Fokas system Equation by Jacobi Elliptic Function Method}\label{Section2}

In this part, we get solitary wave solutions for the (2+1)-dimensional nonlinear coupled Davey-Stewartson Fokas system using the Jacobi Elliptic Function (JEF) Method. This approach, which converts nonlinear partial differential equations into algebraic equations, is a potent analytical tool that is frequently employed to solve them. Using the properties of Jacobi elliptic functions, the method makes it possible to get accurate waveforms, such as periodic and solitary waveforms. This approach makes it easier to investigate the dynamic behaviors and underlying physical phenomena of the coupled Davey-Stewartson Fokas system. The papers \autocite{Ref29,Ref30} provide a thorough description of the method. Remember the ODE \eqref{E9}, which can be written as:
\begin{equation}\label{E14}
        \Phi'' - \frac{\left( k_{{1}}{\lambda_{{1}}}^{2}+\lambda_{{3}} \right)}{k_{{1}}\ {\alpha}^{2}} \Phi - {\frac {k_{{2}}k_{{4}}}{k_{1} k_{{3}}\alpha\gamma}} {\Phi}^{3}=0.
\end{equation}
It is assumed that Eq. \eqref{E9} has the solution of the form.
\begin{equation}\label{E15}
    \Phi(\epsilon) = \sum_{j=0}^{K} m_j R^j(\epsilon).
\end{equation}
Using the balancing approach on the nonlinear $U^3$ terms and highest order derivative $U''$ in \eqref{E14}, we obtain $K = 1$, and the solution form \eqref{E15} can be expressed as:
\begin{equation}\label{E16}
    \Phi(\epsilon) = m_0 + m_1 R(\epsilon),
\end{equation}
where $m_0, m_1$ are the constants to be determined, and the function $R(\epsilon)$ satisfies the ansatz given by:
\begin{equation}\label{E17}
    R'(\epsilon) = \sqrt{c_1 + c_2 R(\epsilon)^2 + \frac{1}{2} c_3 R(\epsilon)^4}.
\end{equation}
The following algebraic system of equations is obtained by first using equation \eqref{E16} together with equation \eqref{E17} in \eqref{E14}:
\begin{align}
    m_{{1}}c_{{3}}-{\frac {k_{{2}}k_{{4}}{m_{{1}}}^{3}}{k_{{1}}\alpha\,k_{
{3}}\gamma}}=0,\\
-3\,{\frac {k_{{2}}k_{{4}}m_{{0}}{m_{{1}}}^{2}}{k_{{1}}\alpha\,k_{{3}}
\gamma}}=0,\\
m_{{1}}c_{{2}}-{\frac {m_{{1}}{\lambda_{{1}}}^{2}}{{\alpha}^{2}}}-{
\frac {m_{{1}}\lambda_{{3}}}{k_{{1}}{\alpha}^{2}}}-3\,{\frac {k_{{2}}k
_{{4}}{m_{{0}}}^{2}m_{{1}}}{k_{{1}}\alpha\,k_{{3}}\gamma}}=0,\\
-{\frac {m_{{0}}{\lambda_{{1}}}^{2}}{{\alpha}^{2}}}-{\frac {m_{{0}}
\lambda_{{3}}}{k_{{1}}{\alpha}^{2}}}-{\frac {k_{{2}}k_{{4}}{m_{{0}}}^{
3}}{k_{{1}}\alpha\,k_{{3}}\gamma}}=0.
\end{align}
Solving this system of equations with the computer program Maple yields the following parameters:
\begin{equation}\label{E23}
    c_1=c_1, \quad c_{{2}}={\frac {k_{{1}}{\lambda_{{1}}}^{2}+\lambda_{{3}}}{k_{{1}}{
\alpha}^{2}}}, \quad c_{{3}}={\frac {k_{{2}}k_{{4}}{m_{{1}}}^{2}}{k_{{1}}\alpha\,k_{{3}}
\gamma}}, \quad m_0 = 0, \quad m_1 = m_1.
\end{equation}
Subject to these results equation \eqref{E16} gives us:
\begin{equation}\label{E24}
    \Phi(\epsilon) = m_1 R(\epsilon).
\end{equation}

\textbf{Solution Set 1: } For the choice of parameters such as \(c_1 = -(1 + \eta^2)\), \(c_2 = 2\eta^2\), and \(c_3 = 1\), the wave profile is derived using the Jacobi elliptic sine function $R(\epsilon)=$\(\operatorname{sn}(\epsilon, \eta)\), producing:
\[
\Phi(\epsilon) = m_1 \operatorname{sn}(\epsilon, \eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the periodic wave profile transitions into a hyperbolic tangent shock wave:
\begin{equation}\label{E24}
    \Phi(\epsilon) = m_1 \tanh{\epsilon},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms as:
\begin{equation}\label{E25}
   \phi_1(\mu, \nu, \tau) = m_{{1}}\tanh \left( \alpha\,\mu+\beta\,\nu-\gamma\tau \right) {{\rm e}^{i
 \left( \mu\,\lambda_{{1}}+\nu\,\lambda_{{2}}+\tau\,\lambda_{{3}}
 \right) }},
\end{equation}
\begin{equation}\label{E26}
   \psi_1(\mu, \nu, \tau) = -{\frac {k_{{4}}\alpha\,{m_{{1}}}^{2} \left( \tanh \left( \alpha\,\mu+
\beta\,\nu-g\tau \right)  \right) ^{2}}{k_{{3}}\gamma}}.
\end{equation}

\textbf{Solution Set 2: } For the choice of parameters such as \(c_1 = 2\eta^2-1\), \(c_2 = 2\), and \(c_3 = -\eta^2 (1-\eta^2)\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\operatorname{ds}(\epsilon, \eta)\), yielding:
\[
\Phi(\epsilon) = m_1 \operatorname{ds}(\epsilon, \eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{dn}(\epsilon, \eta) \to \operatorname{csch}{\epsilon}$
\begin{equation}\label{E27}
    \Phi(\epsilon) = m_1 \operatorname{csch}{\epsilon},
\end{equation}
which on substituting back into Eq. \eqref{E2} gives us:
\begin{equation}\label{E28}
   \phi_2(\mu, \nu, \tau) = m_1 \operatorname{csch}{ (\alpha \mu + \beta \nu - \gamma \tau)} e^i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3),
\end{equation}
\begin{equation}\label{E29}
   \psi_2(\mu, \nu, \tau) = -\frac{k_4 \alpha m_1^2 \operatorname{csch}^2(\alpha \mu + \beta \nu - \gamma \tau)}{k_3 \gamma}.
\end{equation}

\textbf{Solution Set 3: } For the choice of parameters such as \(c_1 = 2-\eta^2\), \(c_2 = 2\),, and \(c_3 = 1-\eta^2\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\operatorname{cs}(\epsilon, \eta)\), producing:
\[
\Phi(\epsilon) = m_1 \operatorname{cs}(\epsilon, \eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{dn}(\epsilon, \eta) \to \operatorname{coth}{\epsilon}$
\begin{equation}\label{E30}
    \Phi(\epsilon) = m_1 \coth{\epsilon},
\end{equation}
substituting back this solution in Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms:
\begin{equation}\label{E31}
   \phi_3(\mu, \nu, \tau) = m_1 \text{csch }(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)},
\end{equation}
\begin{equation}\label{E32}
   \psi_3(\mu, \nu, \tau) = -\frac{k_4 \alpha m_1^2 \text{csch}^2(\alpha \mu + \beta \nu - \gamma \tau)}{k_3 \gamma}.
\end{equation}

\textbf{Solution Set 4: } For the choice of parameters such as \(c_1 = 2\eta^2-1\), \(c_2 = -2\eta^2\), and \(c_3 = 1-\eta^2\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\operatorname{cn}(\epsilon, \eta)\), yielding:
\[
\Phi(\epsilon) = m_1 \operatorname{cn}(\epsilon, \eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{\epsilon}$
\begin{equation}\label{E33}
    \Phi(\epsilon) = m_1 \operatorname{sech}{\epsilon},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms:
\begin{equation}\label{E34}
   \phi_4(\mu, \nu, \tau) = m_1 \text{sech}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)},
\end{equation}
\begin{equation}\label{E35}
   \psi_4(\mu, \nu, \tau) = -\frac{k_4 \alpha m_1^2 \text{sech}^2(\alpha \mu + \beta \nu - \gamma \tau)}{k_3 \gamma}.
\end{equation}

\textbf{Solution Set 5: } For the choice of parameters such as \(c_1 = 2-\eta^2\), \(c_2 = -2\), and \(c_3 = \eta^2-1\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\operatorname{dn}(\epsilon, \eta)\), yielding:
\[
\Phi(\epsilon) = m_1 \operatorname{dn}(\epsilon, \eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{\epsilon}$
\begin{equation}\label{E36}
    \Phi(\epsilon) = m_1 \operatorname{sech}{\epsilon},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms:
\begin{equation}\label{E37}
  \phi_5(\mu, \nu, \tau) =m_1 \text{sech}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)},
\end{equation}
\begin{equation}\label{E38}
   \psi_5(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 \text{sech}^2(\alpha \mu + \beta \nu - \gamma \tau)}{k_3 \gamma}.
\end{equation}
\textbf{Solution Set 6: } For the choice of parameters such as \(c_1 = \frac{\eta^2 -2}{2}\), \(c_2 = \frac{\eta^2}{2}\), and \(c_3 = \frac{1}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\frac{\operatorname{sn}(\epsilon, \eta)}{1 \pm \operatorname{dn}(\epsilon, \eta)}\), yielding:
\[
\Phi(\epsilon) = m_1 \frac{\operatorname{sn}(\epsilon, \eta)}{1 \pm \operatorname{dn}(\epsilon, \eta)}.
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{sn}(\epsilon, \eta) \to \operatorname{tanh}{\epsilon}$ and $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{(\epsilon)}$
\begin{equation}\label{E39}
    \Phi(\epsilon) = m_1 \frac{\operatorname{tanh}(\epsilon)}{1 \pm \operatorname{sech}(\epsilon)},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms:
\begin{equation}\label{E40}
  \phi_6(\mu, \nu, \tau) =\frac{m_1 \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{1 + \text{sech}(\alpha \mu + \beta \nu - \gamma \tau)},
\end{equation}
\begin{equation}\label{E41}
   \psi_6(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 \text{tanh}^2(\alpha \mu + \beta \nu - \gamma \tau)}{(1 + \text{sech}(\alpha \mu + \beta \nu - \gamma \tau))^2 k_3 \gamma}.
\end{equation}
\textbf{Solution Set 7: } For the choice of parameters such as \(c_1 = \frac{\eta^2 -2}{2}\), \(c_2 = \frac{\eta^2}{2}\), and \(c_3 = \frac{\eta^2}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=$\(\frac{\operatorname{dn}(\epsilon, \eta)}{(\eta^2+1)(\operatorname{sn}(\epsilon,\eta) \pm \operatorname{dn}(\epsilon, \eta))}\), yielding:
\[
\Phi(\epsilon) = m_1 \frac{\operatorname{dn}(\epsilon, \eta)}{(\eta^2+1)(\operatorname{sn}(\epsilon,\eta) \pm \operatorname{dn}(\epsilon, \eta))}.
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{sn}(\epsilon, \eta) \to \operatorname{tanh}{\epsilon}$ and $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{(\epsilon)}$
\begin{equation}\label{E42}
    \Phi(\epsilon) = m_1 \frac{\operatorname{sech}(\epsilon)}{2(\operatorname{tanh}(\epsilon) \pm \operatorname{sech}(\epsilon))},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in terms of multidimensional wave transforms:
\begin{equation}\label{E43}
  \phi_7(\mu, \nu, \tau) =\frac{m_1 \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{1 - \text{sech}(\alpha \mu + \beta \nu - \gamma \tau)},
\end{equation}
\begin{equation}\label{E44}
   \psi_7(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 \text{tanh}^2(\alpha \mu + \beta \nu - \gamma \tau)}{(1 - \text{sech}(\alpha \mu + \beta \nu - \gamma \tau))^2 k_3 \gamma}.
\end{equation}
\textbf{Solution Set 8: } For the choice of parameters such as \(c_1 = \frac{\eta^2+1}{2}\), \(c_2 = -\frac{1}{2}\), and \(c_3 = -\frac{(1-\eta^2)^2}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)= \operatorname{cn}(\epsilon,\eta) \pm \operatorname{dn}(\epsilon,\eta)$ yielding:
\[
\Phi(\epsilon) = m_1 \operatorname{cn}(\epsilon,\eta) \pm \operatorname{dn}(\epsilon,\eta).
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{cn}(\epsilon, \eta) \to \operatorname{sech}{\epsilon}$ and $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{(\epsilon)}$
\begin{equation}\label{E45}
    \Phi(\epsilon) = m_1 (\operatorname{sech}(\epsilon) \pm \operatorname{sech}(\epsilon)),
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in multidimensional wave transforms:
\begin{equation}\label{E46}
  \phi_8(\mu, \nu, \tau) =m_1 (\text{tanh}(\alpha \mu + \beta \nu - \gamma \tau) + \text{sech}(\alpha \mu + \beta \nu - \gamma \tau)) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)},
\end{equation}
\begin{equation}\label{E47}
   \psi_8(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 (\text{tanh}(\alpha \mu + \beta \nu - \gamma \tau) + \text{sech}(\alpha \mu + \beta \nu - \gamma \tau))^2}{k_3 \gamma}.
\end{equation} 
\textbf{Solution Set 9: } For the choice of parameters such as \(c_1 = \frac{\eta^2+1}{2}\), \(c_2 = \frac{\eta^2-1}{2}\), and \(c_3 = \frac{\eta^2-1}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)=\frac{\operatorname{dn}(\epsilon,\eta)}{1 \pm \operatorname{dn}(\epsilon,\eta)} $, yielding:
\[
\Phi(\epsilon)=m_1\frac{\operatorname{dn}(\epsilon,\eta)}{1 \pm \operatorname{dn}(\epsilon,\eta)}.
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{dn}(\epsilon, \eta) \to \operatorname{sech}{\epsilon}$ and $\operatorname{sn}(\epsilon, \eta) \to \operatorname{tanh}{(\epsilon)}$
\begin{equation}\label{E48}
    \Phi(\epsilon) = m_1 \frac{\operatorname{sech}(\epsilon)}{1 \pm \operatorname{tanh}(\epsilon)},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in the form of multidimensional wave transforms:
\begin{equation}\label{E49}
  \phi_9(\mu, \nu, \tau) =\frac{m_1 \text{sech}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{1 + \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau)},
\end{equation}
\begin{equation}\label{E50}
   \psi_9(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 \text{sech}^2(\alpha \mu + \beta \nu - \gamma \tau)}{(1 + \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau))^2 k_3 \gamma}.
   \end{equation}
\textbf{Solution Set 10: } For the choice of parameters such as \(c_1 = \frac{\eta^2+1}{2}\), \(c_2 = \frac{1-\eta^2}{2}\), and \(c_3 = \frac{1-\eta^2}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)= \frac{\operatorname{cn}(\epsilon, \eta)}{1 \pm \operatorname{sn}(\epsilon, \eta)}$, yielding:
\[
\Phi(\epsilon) = m_1 \frac{\operatorname{cn}(\epsilon, \eta)}{1 \pm \operatorname{sn}(\epsilon, \eta)}.
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{cn}(\epsilon, \eta) \to \operatorname{sech}{(\epsilon)}$ and $\operatorname{sn}(\epsilon, \eta) \to \operatorname{tanh}{(\epsilon)}$
\begin{equation}\label{E51}
    \Phi(\epsilon) = m_1 \frac{\operatorname{sech}(\epsilon)}{1 \pm \operatorname{tanh}(\epsilon)},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in the form of multidimensional wave transforms:
\begin{equation}\label{E52}
  \phi_{10}(\mu, \nu, \tau) =\frac{m_1 \text{sech}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{1 + \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau)},
\end{equation}
\begin{equation}
    \psi_{10}(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 \text{sech}^2(\alpha \mu + \beta \nu - \gamma \tau)}{(1 + \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau))^2 k_3 \gamma}.
\end{equation}
\textbf{Solution Set 11: } For the choice of parameters such as \(c_1 = \frac{\eta^2+1}{2}\), \(c_2 = \frac{1-\eta^2}{2}\), and \(c_3 = \frac{1-\eta^2}{4}\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)= \frac{\operatorname{sn}(\epsilon, \eta)}{\operatorname{dn}(\epsilon, \eta) \pm \operatorname{cn}(\epsilon, \eta)}$, yielding:
\[
\Phi(\epsilon) = m_1 \frac{\operatorname{sn}(\epsilon, \eta)}{\operatorname{dn}(\epsilon, \eta) \pm \operatorname{cn}(\epsilon, \eta)}.
\]

As the elliptic modulus parameter $\eta \to 1$ the Jacobi elliptic function $\operatorname{sn}(\epsilon, \eta) \to \operatorname{tanh}{(\epsilon)}$ and $\operatorname{dn},\operatorname{cn} (\epsilon, \eta) \to \operatorname{tanh}{(\epsilon)}$
\begin{equation}\label{E54}
    \Phi(\epsilon) = m_1 \frac{\operatorname{tanh}(\epsilon)}{\operatorname{sech}(\epsilon) \pm \operatorname{sech}(\epsilon)},
\end{equation}
substituting back this solution into Eq. \eqref{E2}, we obtain the solutions in the form of multidimensional wave transforms:
\begin{equation}\label{E55}
  \phi_{11}(\mu, \nu, \tau) =\frac{1}{2} \frac{m_1 \text{tanh}(\alpha \mu + \beta \nu - \gamma \tau) e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{\text{sech}(\alpha \mu + \beta \nu - \gamma \tau)},
\end{equation}
\begin{equation}\label{E56}
   \psi_{11}(\mu, \nu, \tau)=-\frac{1}{4} \frac{k_4 \alpha m_1^2 \text{tanh}^2(\alpha \mu + \beta \nu - \gamma \tau)}{\text{sech}^2(\alpha \mu + \beta \nu - \gamma \tau) k_3 \gamma}.
   \end{equation}
\textbf{Solution Set 12: }
For the choice of parameters such as \(c_1 = 0\), \(c_2 = 2\), and \(c_3 = 0\), the wave profile is derived using the Jacobi elliptic function $R(\epsilon)= \frac{D}{\epsilon}$, yielding:
\[
\Phi(\epsilon) = m_1 \frac{D}{\epsilon},
\]
substituting back this solution into Eq. \eqref{E2} to obtain the solutions in the form of multidimensional wave transforms:
\begin{equation}\label{E57}
  \phi_{12}(\mu, \nu, \tau) =\frac{m_1 D e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)}}{\alpha \mu + \beta \nu - \gamma \tau},
\end{equation}
\begin{equation}\label{E58}
   \psi_{12}(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 D^2}{(\alpha \mu + \beta \nu - \gamma \tau)^2 k_3 \gamma}.
   \end{equation}
\textbf{Solution Set 13: } 
For the choice of parameters such as \(c_1 = 0\), \(c_2 = 1\), and \(c_3 = 0\), the wave profile is described using \(R(\epsilon) =  D e^\epsilon\)
\[
\Phi(\epsilon) = m_1 D e^\epsilon,
\]
substituting back this solution into Eq. \eqref{E2} to obtain the solutions in the form of multidimensional wave transforms:
\begin{equation}\label{E59}
  \phi_{13}(\mu, \nu, \tau) =m_1 D e^{\alpha \mu + \beta \nu - \gamma \tau} e^{i (\mu \lambda_1 + \nu \lambda_2 + \tau \lambda_3)},
\end{equation}
\begin{equation}\label{E60}
   \psi_{13}(\mu, \nu, \tau)=-\frac{k_4 \alpha m_1^2 D^2 e^{2(\alpha \mu + \beta \nu - \gamma \tau)}}{k_3 \gamma}.
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
\begin{figure}[H]
    \centering
    % Disable subfigure numbering
    \captionsetup[subfigure]{labelformat=empty}

    % Subfigures for sigma = 0.962
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p1a.png}
        \caption{(a) 3D Portait}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p1b.png}
        \caption{(b)Contour Plot}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p1c.png}
        \caption{(c) 2D-Portrait for $\sigma=0.962$}
    \end{subfigure}
\vspace{-0.3cm}
    \caption{Graphical structure showing anti-Kink soliton solution $\phi_1(\mu,\nu,\tau)$ Eq. \eqref{E25}}
    \label{fig:p_1}
\end{figure}
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[H]
    \centering
    % Disable subfigure numbering
    \captionsetup[subfigure]{labelformat=empty}

    % Subfigures for sigma = 0.962
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p6a.png}
        \caption{(a) 3D Portait}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p6b.png}
        \caption{(b)Contour Plot}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p6c.png}
        \caption{(c) 2D-Portrait for $\sigma=0.962$}
    \end{subfigure}
    \vspace{-0.3cm}
    \caption{Graphical results showing Dark Soliton solution $\phi_6(\mu,\nu,\tau)$ Eq. \eqref{E40}}
    \label{fig:p_1}
\end{figure}
%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[H]
    \centering
    % Disable subfigure numbering
    \captionsetup[subfigure]{labelformat=empty}

    % Subfigures for sigma = 0.962
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p8a.png}
        \caption{(a) 3D Portait}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p8b.png}
        \caption{(b)Contour Plot}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{p8c.png}
        \caption{(c) 2D-Portrait for $\sigma=0.962$}
    \end{subfigure}
\vspace{-0.3cm}
    \caption{Graphical structure showing Bright-Dark soliton solution for $\phi_8(\mu,\nu,\tau)$ Eq. \eqref{E46}}
    \label{fig:p_1}
\end{figure}
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%
\begin{figure}[H]
    \centering
    % Disable subfigure numbering
    \captionsetup[subfigure]{labelformat=empty}

    % Subfigures for sigma = 0.962
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{psi8a.png}
        \caption{(a) 3D Portait}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{psi8b.png}
        \caption{(b)Contour Plot}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.5\textwidth}
        \centering
        \includegraphics[width=\textwidth,height=0.2\textheight,keepaspectratio]{psi8c.png}
        \caption{(c) 2D-Portrait for $\sigma=0.962$}
    \end{subfigure}
\vspace{-0.3cm}
    \caption{Graphical structure showing results for the solution $\psi_8(\mu,\nu,\tau)$ Eq. \eqref{E47}}
    \label{fig:psi_1}
\end{figure}
\section{Conclusion}\label{S5}
We conducted a thorough analysis of the system's dynamics, which improved our understanding of the qualitative characteristics of the (2+1)-dimension coupled Davey-Stewartson Fokas system. We created graphic representations in two and three dimensions to better showcase these discoveries and highlight the intricate architecture of the found information. We used the Jacobi elliptic function method to look into the soliton solutions of the coupled Davey Stewartson Fokas system. Our goals were to find new soliton solutions and understand how the system works dynamically. By using the Jacobi elliptic function method along with dynamical system approaches, we are able to come up with new wave solutions that are consistent with the original PDE and provide us with new information about how the model works. Future research on comparable systems will benefit greatly from the approaches used in this work, which have proven to be successful in resolving challenging nonlinear issues. The findings have potential implications in a number of scientific and engineering fields in addition to advancing theory. In the future, it will be easier to study nonlinear evolution equations with fractional-order derivatives and systems with more dimensions using the methods described here. This study paves the way for a more profound investigation of nonlinear wave phenomena and their practical applications.
\section*{Data Availability} {{For the current study, no data sharing is applicable because no data were investigated or developed for this article.}}
\section*{Conflicts of Interest}{The authors declare that they have no conflict of interest.}
\section*{Author's contribution}
{Muhammad Hammad: Contributed to methodology development, software implementation, and manuscript writing, including reviewing and editing. Amjad Hussain: Responsible for conceptualization, original draft preparation, editing, and conducting the investigation.}

\printbibliography
\end{document}
