% Modelo "Strogatz" (solución)
clc
close all
clear all

disp('Sistema de 2 Ecuaciones Diferenciales (Coordenadas Cartesianas)')
disp('Modelo "Strogatz"')
 
disp('...................................................................')
 
disp('Ecuaciones: ')
 
disp('dx/dt = x*(1-x^2-y^2)-w*y') 
disp('dy/dt = y*(1-x^2-y^2)+w*x')
disp('r = (x^2 + y^2)^0.5') 
disp('Theta = Arctg(y/x)')

disp('...................................................................')
 
% 1) Introducir: 

%A) Una matriz de ceros llamada "z" de dos columnas, donde los elementos de las mismas representen los sucesivos valores 
%de las variables "x" e "y" (1º y 2º columna, respectivamente). 

z=zeros(1,2);

%B) Los comandos necesarios para indicar los valores de las condiciones iniciales (x0, y0, tiempo final (tf), veloicidad angular (w)) 
%desde la Command Window.

z(1)=input('Valor inicial de x (x0): ');
z(2)=input('Valor inicial de y (y0): ');
tf=input('Tiempo final, tf: ');
w=input('Velocidad angular, w: '); %Ej. 2*pi

% 2) 
%A) Definir el sistema de ecuaciones diferenciales "fg" para "dx/dt" y "dy/dt", donde z(1) sea "x" y z(2) sea "y".

fg=@(t,z) [z(1)*(1-z(1)^2-z(2)^2)-w*z(2);z(2)*(1-z(1)^2-z(2)^2)+w*z(1)]; % z(1) es "x", z(2) es "y"

%B) Resolver numéricamente el sistema de ecuaciones diferenciales "fg" utilizando el solver ode45. 

tspan=[0 tf]; 
[t,z]=ode23s(fg,tspan,z);

% 3) Graficar los cursos temporales de x e y vs. tiempo (individualmente en un gráfico y superpuestos en otro).

% "x" vs. "t"
figure
plot(t,z(:,1),'r','LineWidth',2)
grid
xlabel('t')
ylabel('x');
title('x vs. t')

% "y" vs. "t"
figure
plot(t,z(:,2),'b','LineWidth',2)
grid
xlabel('t')
ylabel('y');
title('y vs. t')

% "x" e "y" vs. "t"
figure
plot(t,z(:,1),'r',t,z(:,2),'b','LineWidth',2)
grid
xlabel('t')
ylabel('x,y');
legend('x vs t','y vs t','Location','northeast')
title('x e y vs. t')

% 4) 
%A) Calcular el valor del radio (r) para cada instante de tiempo.

r=(z(:,1).^2+z(:,2).^2).^0.5; % r = (x^2 + y^2)^(1/2)

%B) Graficar r vs. tiempo. Titular el gráfico y nombrar adecuadamente los ejes. 

r=(z(:,1).^2+z(:,2).^2).^0.5; % r = (x^2 + y^2)^(1/2)
figure
plot(t,r,'m','LineWidth',2)
grid
xlabel('t')
ylabel('r');
title('r vs. t')

% 5)
%A) Calcular el valor del ánglo Theta para cada instante de tiempo.

Theta=atan(z(:,2)./z(:,1)); % Theta = Arctg(y/x)

%B) Graficar Theta vs. tiempo. Titular el gráfico y nombrar adecuadamente los ejes. 

% "Theta" vs. "t"
figure
plot(t,Theta,'k','LineWidth',2)
grid
xlabel('t')
ylabel('Theta');
title('Theta = Arctg (y/x) vs. t')


% 6) 
%A) Graficar el plano de fases y vs X.

% y vs. x  (Plano de Fases)
figure
grid
hold on
plot(z(:,1),z(:,2),'g','LineWidth',2)

%B) Marcar la posición inicial y final de las coordenadas x e y, distinguéndolas con diferentes colores/símbolos. 

plot(z(1,1),z(1,2),'o','MarkerEdgeColor','k',...
             'MarkerFaceColor','r',...
             'MarkerSize',7)
plot(z(end,1),z(end,2),'o','MarkerEdgeColor','r',...
             'MarkerFaceColor','b',...
             'MarkerSize',7)

%C) Limitar, cuando sea necesario, los ejes del gráfico para obtener una buena visulaización del mismo 
%(sobre todo necesario cuando se parte desde "fuera" del ciclo límite, sino se puede comentar y no tener en cuenta). 

%Ejemplo para X0=2 y0=-2
axis([-3 3 -3 3])

%D) Nombrar adecuadamente los ejes del gráfico y titular el mismo. 
xlabel('x');
ylabel('y');
title('y vs. x (Plano de Fases)')

% 7) Simular para varios conjuntos de parámentros, teniendo en cuenta diferentes valores de x0 e y0 (tanto de "adentro" 
%como de "afuera" del ciclo límite), así como para diferentes valores de tiempo y de w.