% vdp2.m
% integration de     dx/dt =  y + x - x**3   ,     dy/dt =  -x
%       schema de Heun 
% la solution exacte est l attracteur de Van der Pol
% fd, iacs, 18 septembre  2004, 22 janvier 2006, 25 fevrier 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel ('vdp2. Cycle limite non lineaire avec le schema de Heun.') ;
%%% title ('Oscillateur de Van der Pol') ;
% title ('Cnam Paris. Travaux pratiques Iacs.') ;
n=10 ;  dt=1/n
x=0.01 ,  y=0 , 

hold off 
	  plot ( x, y , 'ob;depart;')  ; hold on 
	  
niter = 6*n 
  traj=zeros(niter+1,2) ; 
traj(1,1) =  x ; traj(1,2) = y ;

for i=1:niter ;
px1 =    y + x - x**3    ; py1 =   - x  ; 
x1  = x + dt * px1       ; y1  = y + dt * py1 ;
px2 = y1 + x1 - x1**3     ;   py2 =   - x1  ; 
xn  = x + 0.5 * dt * (px1+px2)   ;   yn  = y +  0.5 * dt * (py1+py2) ;
x=xn ; y=yn ; 
traj(i+1,1) = xn  ; traj(i+1,2) = yn ;
end ; 
xn , yn 
  traj ; 

plot ( traj(:,1) , traj(:,2), '*k;trajectoire;') ; hold on 
axis ([ -.12 , .025 ,  -.05 , .22  ]) 
title ('Oscillateur de Van der Pol : schema de  Heun') ;
	  plot ( xn, yn , '*r;arrivee;') ; hold off 
	  
%%% print -dpdf 03-vdp2.pdf
%%% print -djpg 03-vdp2.jpg
print('03-vdp2.jpg', '-djpg', '-S891,630');

  


