// TP 2 Activité 2


function [dY] = vanderpol(t,Y)
dY(1) = c * (1-Y(2)^2) * Y(1) - Y(2)
dY(2) = Y(1)
endfunction

//==============================================================
// METHODE POUR RESOUDRE UN SYSTEME AVEC RUNGE-KUTTA ORDRE 4
//==============================================================
function [T,Y] = runge4sys(f, t0,y0, tmax, N)
 // Renvoie dans T un vecteur de N valeurs régulièrement espacées 
 // entre  t0  et  tmax 
 T = linspace(t0,tmax,N);
 h = (tmax-t0)/N;
  // Nb d equations dans le systeme
 nbeq = length(y0);
 // point initial
 Y = zeros(N,nbeq);
 Y(1,:)=y0'
 // Definit la taille de k
 k=zeros(4, nbeq);
 //--------------------------------------
 // Calcul itératif des N points 14
 //--------------------------------------
 for i=1:(N-1)
  // Calcul de k(j,1)
  k1 = h * f(T(i), Y(i,:));  
  // Calcul de k(j,2)
  tmp_Y = Y(i,:)+k1'/2;
  k2 = h * f(T(i), tmp_Y);
  // Calcul de k(j,3)
  tmp_Y = Y(i,:)+k2'/2;
  k3 = h * f(T(i), tmp_Y);  
  // Calcul de k(j,4)
  tmp_Y = Y(i,:)+k3';
  k4 = h * f(T(i), tmp_Y);
  // Calcul de y(i+1) en fonction de y(i) et k(j,i)
  Y(i+1,:) = Y(i,:)+1/6*(k1+2*k2+2*k3+k4)' ;
 end
endfunction


function gradient()
// Tracé du champs de vecteur issu de l’équation de VAN DER POL : 
n = 30; dx = 6; dy = 4;
x = linspace(-dx,dx,n);
y = linspace(-dy,dy,n);
clf();
fchamp(vanderpol,0,x,y,1,[-dx,-dy,dx,dy],"031")
xselect()
endfunction

function solution()
// 2. Résolution de l’équation différentielle :
m = 500 ; T = 30; t = linspace(0,T,m);
couleurs = [21 2 3 4 5 6 19 28 32 9 13 22 18 21 12 30 27]; // 17 couleurs
num = -1;
while %t
[c_i,c_x,c_y]=xclick();
if c_i == 0 then
plot2d(c_x, c_y, -9, "000")
u0 = [c_x;c_y];
[u] = ode(u0, 0, t, vanderpol);
num = modulo(num+1,length(couleurs));
plot2d(u(1,:),u(2,:),couleurs(num+1),"000")
elseif c_i == 2 then
break
end
end
endfunction