%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%       Antes de rodar este arquivo é necessário       %%%%%%%%%%%
%%%%%%%%%   primeiramente rodar o arquivo "principalnegativo"  %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all; clc
global mi
syms x1 x2

% Dinâmica do sistema oscilador de Van der Pol
x1ponto = x2;
x2ponto = -x1 + mi*(1-x1^2)*x2;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nesta parte é utilizado as variáveis x1 e x2 como variáveis simbólicas 
% para poder utilizar o comando ezsurf, que plota a superfície da função.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Funções de energia propostas 
v1 = .5*x1^2+.5*x2^2;
v2 = .5*(x1^2+(mi*((x1^3)/3-x1)+x2)^2);

figure (1)
ezsurf(v1, [-3,3], [-3,3]) % plota a superfície de V1
zlabel('Função energia V1')

figure (2)
ezsurf(v2, [-10,10], [-10,10]) % plota a superfície de V2
zlabel('Função energia V2')

% Cálculo das derivadas das funções energia
vponto1 = x1ponto*diff(v1,x1) +x2ponto*diff(v1,x2);
vponto2 = x1ponto*diff(v2,x1) +x2ponto*diff(v2,x2);

figure (3)
ezsurf(vponto1, [-3,3], [-3,3]) % plota a superfície de V1ponto
zlabel('V1ponto')

figure (4)
ezsurf(vponto2, [-10,10], [-10,10]) % plota a superfície de V2ponto
zlabel('V2ponto')

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Agora são calculadas as funções energia V1 e V2 e também as suas respectivas derivadas
% utilizando os valores de x1 e x2 calculados pela ODE do arquivo "principalnegativo".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Funções de energia propostas 
V1 =.5*x(:,1).^2+.5*x(:,2).^2;
V2 = .5*(x(:,1).^2+(mi*((x(:,1).^3)/3-x(:,1))+x(:,2)).^2);

V1ponto = x(:,1).*x(:,2)-x(:,2).*(x(:,1)+x(:,2).*(x(:,1).^2-1));
V2ponto = x(:,2).*(x(:,1)+(x(:,1).^2-1).*(x(:,1).^3/3-x(:,1)+x(:,2)))-(x(:,1)+x(:,2).*(x(:,1).^2-1)).*(x(:,1).^3/3-x(:,1)+x(:,2));

% Gráficos 

% resposta da função energia V1 no tempo
figure(5)
plot(tempo,V1) 
xlabel('Tempo [s]')
ylabel('V1')
grid

% função energia V1 com as variáveis de estado x1 e x2
figure(6)
plot3(x(:,1),x(:,2),V1) 
hold on
xlabel('x1')
ylabel('x2')
zlabel('V1')
grid

% resposta de V1ponto no tempo
figure(7)
plot(tempo,V1ponto) 
xlabel('Tempo [s]')
ylabel('V1ponto')
grid

% V1ponto com as variáveis de estado x1 e x2
figure(8)
plot3(x(:,1),x(:,2),V1ponto) 
hold on
xlabel('x1')
ylabel('x2')
zlabel('V1ponto')
grid

% resposta da função energia V2 no tempo
figure(9)
plot(tempo,V2) 
xlabel('Tempo [s]')
ylabel('V2')
grid

% função energia V2 com as variáveis de estado x1 e x2
figure(10)
plot3(x(:,1),x(:,2),V2) 
hold on
xlabel('x1')
ylabel('x2')
zlabel('V2')
grid

% resposta de V2ponto no tempo
figure(11)
plot(tempo,V2ponto) 
xlabel('Tempo [s]')
ylabel('V2ponto')
grid

% V2ponto com as variáveis de estado x1 e x2
figure(12)
plot3(x(:,1),x(:,2),V2ponto) 
hold on
xlabel('x1')
ylabel('x2')
zlabel('V2ponto')
grid