in order to solve coupled partial differential equation using FFT , my code is showing error as -

Error using horzcat Requested 240x1094476 (3.9GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

Error in ode45 (line 484) yout = [yout, zeros(neq,chunk,dataType)];

Error in odefun2 (line 54) [z, E_omega] = ode45(@odefun,zspan,E_omega);

>> for the code-

%first electric field envelop in fourier space tic; clc clear all close all Po=9;

Ao=sqrt(Po); C=0; %pi=3.1415926535; x=-60*10^-15:1*10^-15:59*10^-15;

x = 10*x; i=sqrt(-1); Wo=120*10^-15;%initial pulse width in second u=Ao*exp(-(x./Wo).^2) figure(1) plot(abs(u),'b'); title('Input Pulse'); xlabel('Time'); ylabel('Amplitude'); hold on grid on E1_omega=fft(fftshift(u)); disp(E1_omega)

%second electric field envelop in fourier space P1=9; A1=sqrt(P1); C=0; %pi=3.1415926535; % x =-60*10^-15:1*10^-15:59*10^-15 i=sqrt(-1); Wo=120*10^-15;%initial pulse width in second u=A1*exp(-(x./Wo).^2) figure(2) plot(abs(u),'b'); title('Input Pulse'); xlabel('Time'); ylabel('Amplitude'); hold on grid on E2_omega=fft(fftshift(u)); disp(E2_omega)

% patching of the E1 and E2 electric fields with: E_omega=[E1_omega,E2_omega] disp(E_omega) figure(3) plot(abs(E_omega),'r');

zspan=[0 30*10^-3] E_omega_0=[0 0.0625] [z, E_omega] = ode45(@odefun,zspan,E_omega);

E1_omega = E_omega(:, 1:end/2); E2_omega = E_omega(:, end/2+1:end); E1_t = ifft(E1_omega, [], 2); E2_t = ifft(E2_omega, [], 2);

figure(101) subplot(2, 2, 1) [X, Y] = meshgrid(1:size(E1_omega, 2), z); mesh(X, Y, abs(fftshift(E1_omega, 2)).^2); subplot(2, 2, 2) [X, Y] = meshgrid(1:size(E1_omega, 2), z); mesh(X, Y,abs(fftshift(E2_omega, 2)).^2); subplot(2, 2, 3) [X, Y] = meshgrid(x, z); mesh(X, Y,abs(fftshift(E1_t, 2)).^2); subplot(2, 2, 4) [X, Y] = meshgrid(x, z); mesh(X, Y,abs(fftshift(E2_t, 2)).^2);

% g=angle(E_omega) % K=unwrap(g) figure(4) plot(z,E_omega,'linewidth',4) xlabel('z') ylabel('phase (radian)') toc

IN WHICH ODE FUNCTION IS DEFINED AS-

function dE_omega_dz = odefun(z, E_omega,~,~) dE_omega_dz=zeros(length(E_omega),1);

Dk=20.94; l=800; % lambda 800nm c=3*10^8; pi=3.1415926535; omega=2*pi*c/l; n_2=5*10^-16; %5*10^-16cm^2/W L_NL=3.7245*10^-3; LGVM=0.6*10^-3;

I0=50*10^9;

to=125e-12; % initial pulse widthin second dt=1e-12/to;

dw=1/l/dt*2*pi; w=(-1*l/2:1:l/2-1)*dw;

% you would have to split the fields back in two: E1_omega = E_omega(1:end/2); E2_omega = E_omega(end/2+1:end);

% go back to time space to calculate the nonlinear part: E1_t = ifft(E1_omega); E2_t = ifft(E2_omega);

% and calculate the derivatives: dE_omega_dz(1:length(E_omega)/2) = fft(1i*conj(E1_t).*E2_t.*exp(1i*Dk*z) ... + 1i*2*pi*n_2*I0*L_NL/l*(abs(E1_t.^2 + ... 2*abs(E2_t.^2)).*E1_t)); dE_omega_dz(length(E_omega)/2+1:length(E_omega)) = 1i*omega*L_NL/LGVM * E2_t + fft(1i*E1_t.*E1_t.*exp(-1i*Dk*z).... + 1i*4*pi*n_2*I0*L_NL/l*(2*abs(E1_t.^2 + ... abs(E2_t.^2)).*E1_t));

end