%
% Représentation d'un système d'ordre 2 dans le plan de phase
% - tracé du champ de vecteurs
% - tracé de trajectoires
%

% plages valeurs en x1 et x2
x1 = linspace(-pi/2,pi/2,30);
x2 = linspace(-3,3,30);
[x,y] = meshgrid(x1,x2);

% definition du system
k=1; g = 9.8; l=2; m=1;
f = @(t,X) [X(2); -k/m*X(2)-g/l*sin(X(1))]; 

% tracé du champ de vecteurs
u = zeros(size(x));
v = zeros(size(x));
for i = 1:numel(x)
    Xpoint = f(0,[x(i); y(i)]);
    u(i) = Xpoint(1);
    v(i) = Xpoint(2);
end

quiver(x,y,u,v,1); 
ylabel('$x_2$','FontSize',20,'Interpreter','latex');
xlabel('$x_1$','FontSize',20,'Interpreter','latex');
axis([min(x1) max(x1) min(x2) max(x2)]);
grid on; 
hold on;


% tracé d'une solution dans le plan
X0 = [1;1];
[t,X] = ode45(f,[0 10],X0);
plot(X(:,1),X(:,2),'r','LineWidth',2);
plot(X(1,1),X(1,2),'*r','LineWidth',2,'MarkerSize',6) 

for i=1:5
X0 = ginput(1);
[t,X] = ode45(f,[0 10],X0);
plot(X(:,1),X(:,2),'r','LineWidth',2);
plot(X(1,1),X(1,2),'*r','LineWidth',2,'MarkerSize',6) 
end


