%
% 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(-3,3,30);
x2 = linspace(-3,3,30);
[x,y] = meshgrid(x1,x2);

% definition du system
k=1; g = 9.8; l=1; 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

figure1 = figure; axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); grid(axes1,'on'); set(axes1,'FontSize',14);
ylabel('$x_2$','FontSize',20,'Interpreter','latex');
xlabel('$x_1$','FontSize',20,'Interpreter','latex');
quiver(x,y,u,v); 
axis tight equal; 
hold on;


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

% seconde trajectoire
X0 = [-0.3;0.0];
[ts,ys] = ode45(f,[0 10],X0);
plot(ys(:,1),ys(:,2),'b','LineWidth',2);
plot(ys(1,1),ys(1,2),'*b','LineWidth',2,'MarkerSize',6);