% vdp4.m
% integration de     dx/dt =  y + x - x**3   ,     dy/dt =  -x
%       schema de Runge et Kutta 
% la solution exacte est l attracteur de Van der Pol
% fd, iacs, 19 decembre 2004, 16 janvier 2005, 25 mai 08 
%                          il y avait un bug dans le rk4...  
% fd  02 mars 2025   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% xlabel ('vdp4. Cycle limite non lineaire avec le schema de Runge et Kutta.') ;
%% title ('Oscillateur de Van der Pol') ;
% xtitle ('Cnam Paris. Travaux pratiques Iacs.') ;
n=10 ;  dt=1/n
 lmax = 4 ;
hold off 
for l=1:lmax;
if l<=1 ;   hold off ; else ; hold on ; end
if (l==1) ;  x=0.01 ; y=0 ; imax = 20*n ; end ; 
if (l==2) ;  x=1    ; y=-1; imax = 10*n ; end ; 
if (l==3) ;  x=.5   ; y=.5; end ;  
if (l==4) ;  x=-.5  ; y=0;  end ;


traj=zeros(imax+1,2) ; traj(1,1) =  x ; traj(1,2) = y ;

for i=1:imax ;
if (i==1) ; plot ( x, y , 'ob;depart;') ;
if (l==1) ; hold on ; end ; end ; 
px1 = y + x - x**3             ;  py1 = - x  ; 
x1  = x +  0.5 *dt * px1       ;  y1  = y + 0.5 * dt * py1 ;
px2 = y1 + x1 - x1**3          ;  py2 =   - x1  ; 
x2  = x +  0.5 *dt * px2       ;  y2  = y + 0.5 * dt * py2 ;
px3 = y2 + x2 - x2**3          ;  py3 = - x2  ; 
x3  = x + dt * px3             ;  y3  = y + dt * py3 ;
px4 = y3 + x3 - x3**3          ;  py4 = - x3  ; 
xn  = x + dt * ( px1 + 2*px2 + 2*px3+ px4  ) / 6 ;
yn  = y + dt * ( py1 + 2*py2 + 2*py3+ py4  ) / 6 ;
x=xn ; y=yn ; 
traj(i+1,1) = xn  ; traj(i+1,2) = yn ;
%
end ; xn , yn
xn , yn 
  traj ; 

plot ( traj(:,1) , traj(:,2), '*k;trajectoire;') ; hold on 
axis ([ -.12 , .025 ,  -.05 , .22  ]) 
axis ([ -1.5 , 1.9,  -1.7 , 1.5 ]) 
	  plot ( xn, yn , '*r;arrivee;') ; 
end
hold off 
	  
	 
title ('Oscillateur de Van der Pol : quatre trajectoires avec le schema de  Runge et Kutta') ;
%% print -dpdf 05-vdp4.pdf
%% print -djpg 05-vdp4.jpg
print('05-vdp4.jpg', '-djpg', '-S891,630');





