num

[soumeyahakima]

%% 4-Component Coupled Nonlinear Schrödinger System (p=4)
clear; clc; close all;

% Spatial domain
Nx = 256; % number of spatial points
x = linspace(-10,10,Nx);
dx = x(2)-x(1);

% Temporal domain
dt = 0.001; % time step
Nt = 500; % number of time steps
t = (0:Nt-1)*dt;

% Coupling coefficients (symmetric positive)
A = [1 0.5 0.3 0.2;
0.5 1 0.4 0.3;
0.3 0.4 1 0.5;
0.2 0.3 0.5 1];

p = 4; % exponent

% Initial conditions (Gaussian-type)
u1 = exp(-x.^2);
u2 = exp(-x.^2/2);
u3 = exp(-x.^2/4);
u4 = exp(-x.^2/8);

% FFT wave numbers
k = (2pi/(Nxdx)) * [0:Nx/2-1 -Nx/2:-1];
K2 = (1i*k).^2; % second derivative in Fourier space

% Precompute linear operator
L = exp(-1iK2dt);

% Allocate space to save snapshots for plotting
U1 = zeros(Nx, Nt); U2 = U1; U3 = U1; U4 = U1;
U1(:,1) = u1.'; U2(:,1) = u2.'; U3(:,1) = u3.'; U4(:,1) = u4.';

% Time evolution (Split-Step Fourier Method)
for n = 2:Nt
% Nonlinear step (half)
NL1 = sum(A(1,:)'.abs([u1;u2;u3;u4]).^(p) . abs(u1).^2).*u1;
NL2 = sum(A(2,:)'.abs([u1;u2;u3;u4]).^(p) . abs(u2).^2).*u2;
NL3 = sum(A(3,:)'.abs([u1;u2;u3;u4]).^(p) . abs(u3).^2).*u3;
NL4 = sum(A(4,:)'.abs([u1;u2;u3;u4]).^(p) . abs(u4).^2).*u4;

u1 = u1.*exp(-1i*NL1*dt/2);
u2 = u2.*exp(-1i*NL2*dt/2);
u3 = u3.*exp(-1i*NL3*dt/2);
u4 = u4.*exp(-1i*NL4*dt/2);

% Linear step (FFT)
u1 = ifft(L.*fft(u1));
u2 = ifft(L.*fft(u2));
u3 = ifft(L.*fft(u3));
u4 = ifft(L.*fft(u4));

% Nonlinear step (half)
NL1 = sum(A(1,:)'.*abs([u1;u2;u3;u4]).^(p) .* abs(u1).^2).*u1;
NL2 = sum(A(2,:)'.*abs([u1;u2;u3;u4]).^(p) .* abs(u2).^2).*u2;
NL3 = sum(A(3,:)'.*abs([u1;u2;u3;u4]).^(p) .* abs(u3).^2).*u3;
NL4 = sum(A(4,:)'.*abs([u1;u2;u3;u4]).^(p) .* abs(u4).^2).*u4;

u1 = u1.*exp(-1i*NL1*dt/2);
u2 = u2.*exp(-1i*NL2*dt/2);
u3 = u3.*exp(-1i*NL3*dt/2);
u4 = u4.*exp(-1i*NL4*dt/2);

% Save snapshots
U1(:,n) = u1.'; U2(:,n) = u2.'; U3(:,n) = u3.'; U4(:,n) = u4.';

end

%% Plot 3D figure for u1
figure;
mesh(t, x, abs(U1));
xlabel('t'); ylabel('x'); zlabel('|u_1(x,t)|');
title('3D Evolution of |u_1(x,t)|');
shading interp;

% Save as PNG
saveas(gcf,'u1_3D.png');