%
% Programa de chamada para o ODE solver
% 
clear;
clc;
% chama o metodo numerico
[T,Y] = ode15s('vdp1000',[0 3000],[2 0]);
% graficos
subplot(3,1,1), plot(T,Y(:,1),'-b'), title('Primeira componente'), xlabel('t'), ylabel('Y(1)')
subplot(3,1,2), plot(T,Y(:,2),'-g'), title('Segunda componente'), xlabel('t'), ylabel('Y(2)')
subplot(3,1,3), plot(Y(:,1),Y(:,2),'r'), title('Retrato de fase'), xlabel('Y(1)'), ylabel('Y(2)')
axis([min(Y(:,1)) max(Y(:,1)) min(Y(:,2)) max(Y(:,2))])