% vdp5.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, 18 septembre, 19 decembre 2004, 2 janvier 2005, 18 oct 2006.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel ('vdp5. Schema de Runge et Kutta d ordre 4.') ;
title ('Oscillateur de Van der Pol') ;
% title ('Cnam Paris. Travaux pratiques Iacs.') ;
n=10 ;  dt=1/n
vec=zeros(2,2) ; 
x=0.01 ; y=0;
x=1    ; y=-1;
x=.5   ; y=.5;
x=-.5  ; y=0;
gset terminal postscript ;   gset output "vanderpol01.ps" ;
axis ([ -1.5 , 1.9,  -1.7 , 1.5 ], 'square') % gset xrange[-1.5: 1.9] ;  gset yrange[-1.7:1.5] ; 
for i=1:10*n;
if i<=1 ;   hold off ; plot(x , y, '*g;depart;') ;  hold on ; 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 ;
%
vec(1,1)=x ; vec(2,1)=xn ; 
vec(1,2)=y ; vec(2,2)=yn ;   
if i<=1 ;   
plot ( vec(:,1) , vec(:,2), '*r;trajectoire;') ; else ; 
plot ( vec(:,1) , vec(:,2), '*r;;') ; end ; 
x=xn ; y=yn ; 
end ; xn , yn 
hold on ;  plot(xn , yn , '*b;arrivee;') ; hold off ;  
