#
#  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

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);
y
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