MATLAB: The code is entering infinite loop when i am trying to solve differential equation using ode45 command , can some one help me fix this

infinite loop coming while solving using ode45 command.

 function dE_dz =asim(z,E1,~,~)
 c=3*1e8;
L=800*10^-9
pi=3.1415926535;
Dk=20.94;
w1=(2*pi*c)/L;
%w1=fs*[-N/2:N/2-1];
deff=2.0657e-12;
Lc=2*pi/Dk;
Lnl=50*Lc;
n2=5*10^-16;
I0=10^-9;
 dE_dz = (-1i*conj(E1).*E1.^2.*w1*deff)/(c*n2*Dk)+ 1i*2*pi*(n2*I0)*(Lnl/L)*((abs(E1)).^2).*E1 +...
    + 1i*4*pi*(n2*I0)*(Lnl/L).*((-w1*deff.*E1.^2*exp(-1i*Dk*z))/(c*n2*Dk)).^2.*E1; 
end
 clc; clear all; close all; clf;
 Pi=1.98;                               %input('Enter the value of input power in mW  ')
t=120*10^-15                           %input('Enter the value of input pulse width in seconds  ')
tau=-300*10^-15:0.1*10^-15:300*10^-15;                         %input('Enter time period with upper(U), lower(L) and interval between upper and  lower interval(I) in this format L:I:U')
dt=1e-15/t
fs=1/dt;
pi=3.1415926535;
Ao=sqrt(Pi)                            % amplitude of an input gaussian pulse
h=2000                                 % small step as defined by  split step algorithm
                                       %for ii=0.1:0.1:(s/10)              %1*10^-7:1*10^-7:s %different fiber lengths
E1=Ao*exp(-(tau/t).^2); % generation of an gaussian input pulse
figure(1)
plot(tau,abs(E1));                         % graph of input pulse
title('Input Pulse'); xlabel('Time in ps'); ylabel('Amplitude');
grid on;
hold on;
 [z,E1]=ode45(@asim,[0 1.0],0.6325)
%q = ode45(@(z,E1)asim(z,E1,~,~),[0 1],E1_0);
% time=
% % % tau1=L*(-300*10^-15:300*10^-15);
plot(z,E1);

Best Answer

Hi,
You make your code inefficient by calculating constant values again and again in every iteration. Try this:
syms E1 z
c=3*1e8;
L=800*10^-9;
pi=3.1415926535;
Dk=20.94;
w1=(2*pi*c)/L;
deff=2.0657e-12;
Lc=2*pi/Dk;
Lnl=50*Lc;
n2=5*10^-16;
I0=10^-9;
dE_dz = (-1i*conj(E1).*E1.^2.*w1*deff)/(c*n2*Dk)+ 1i*2*pi*(n2*I0)*(Lnl/L)*((abs(E1)).^2).*E1 +...
    + 1i*4*pi*(n2*I0)*(Lnl/L).*((-w1*deff.*E1.^2*exp(-1i*Dk*z))/(c*n2*Dk)).^2.*E1; 
matlabFunction(dE_dz,'File','asim_fast')
This code creates a file asim_fast.m containing this code - which does the same calculation much faster:
function dE_dz = asim_fast(E1,z)
%ASIM_FAST
%    DE_DZ = ASIM_FAST(E1,Z)
%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    10-Oct-2018 20:02:27
t2 = E1.^2;
t3 = t2.^2;
t4 = abs(E1);
dE_dz = E1.*t4.^2.*5.891597660294395e-17i-t2.*conj(E1).*1.549567321951994e9i+E1.*t3.*exp(z.*-4.188e1i).*2.829332414080319e2i;
Change your ode call to:
Pi=1.98;
t=120*10^-15;
tau=-300*10^-15:0.1*10^-15:300*10^-15;
dt=1e-15/t;
fs=1/dt; 
Ao=sqrt(Pi);
h=2000;
E1=Ao*exp(-(tau/t).^2);
figure(1)
subplot(2,1,1)
plot(tau,abs(E1));
title('Input Pulse');
xlabel('Time in ps');
ylabel('Amplitude'); 
grid on; 
[z,E1]=ode45(@asim_fast,[0 1.0],0.6325);
subplot(2,1,2)
plot(z,E1)
This will significant increase speed.
I used subplot to better see the both diagrams. Since the result is changing mainly in the imaginary part perhaps another kind of plot (e.g. z over abs(E1) would make more sense - but i have no insight to your problem.
Best regards
Stephan
Related Question