[ODE45 ISSUE] Problem when plotting

[Marcus Sandri]

Hello,

          I'm new in Octave and I'm trying to plot some ODEs. I did read some tutorials and some others topics in this list. So, I wrote a code to plot following ODE:

#  ODE'S to plot...
#
#  dx                      ____
#  ---- = U_0 * N * C \           C_n-1    *  C_n  + (delta - Kappa) * (alpha - alpha_0)
#  dt                       /___      
#               
#
#
#
#  dCn
#  ---- = -4i (bf * t)^2  C_n    + U_0 * (alpha * C_n+1   - alpha * C_n+1)
#   dt                                          
#
#
#5 states (C_n+1,C_n,C_n-1,C_n-2,C_n-3, alpha)
# 
#Initial condictions:
#
# C_n(0) = 1, others_states must be equals: 0
# alpha(0) = 0
# delta = 0
# kappa = 160
# bf = 0.035
# N = 20000
# U_0 = 0.04
# alpha_0 = 20

My code is:

function testing5modes()
global N delta kappa Nmo bf U_0 alpha_0 C
 
#clear all
 
delta=0;
kappa=160;
Nmo=5; %if 5, n+1,n,n-1,n-2,n-3
bf=0.035;
N=20000;
tf=28; % Final time
C=0; %C=0 Bloch, C=1 CARL
U_0=0.04;
alpha_0=20;
 
Ndiag=1000; % Nb of points for the time sampling
tspan=(0:Ndiag)/Ndiag.*tf; % Interval of integration
 
h=1e-6; % Numerical tolerance
options = odeset('RelTol',1e-4,'AbsTol',[1e-4]);
 
y0=zeros(Nmo+1,1);
y0(Nmo+1)=0;
y0(Nmo-3)=1;

 
[t,y]=ode45(@eqn,tspan,y0,options);
 
plot(t,y)
 
endfunction 
 
%*************************************************************************
function dy = eqn(t,y)
global N delta kappa Nmo bf U_0 alpha_0 C
 
dy=zeros(Nmo+1,1);
alpha=y(Nmo+1);
dy(Nmo+1)=U_0.*N.*sum((Nmo.* Nmo-1) +(1i.* delta-kappa).*(alpha-alpha_0));
dy(1:Nmo)=-1i.*4.*(bf.*t).^2.*y(1:Nmo)+U_0.*(conj(alpha).*[y(2:Nmo); 0]-alpha.*[0; y(1:Nmo-1)]);
 
endfunction

Problem:

                         When I plot (t, y), it just plot a "vertical" line. I'd checked too the ' y ' vector and it output just the firsts 3 numbers in the matrix. After that it just show NaN (Not a Number).  Where Am I wrong in this code?

PS: Full code attached

Cheers,

Marcs.

testing.m
Description: Text Data