%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         Practice 4          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear Dynamical Systems %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2019.10.09.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Van der Pol oszcillator:
%ddx-mu*(1-x^2)*dx+x=0 ->
%dx=y -> ddx=dy -> dy=mu*(1-x^2)*y-x
close all;

tmax=100;%1000;% 

for mu = 0:0.1:5% mu = 0:50:500%  %calculate the equation for different mu parameters, ode45 gets slow for great mu values
    
    VanDerPol_egyenlet=(@(t, xy)[xy(2); mu*(1-xy(1)^2)*xy(2)-xy(1)]); 

init_xy1=[0.5 0.5];
[t, fi] = ode45(VanDerPol_egyenlet, [0,tmax], init_xy1); %general method
%[t, fi] = ode15s(VanDerPol_egyenlet, [0,tmax], init_xy1); %ode15s can handle very different fast and slow times at the same system (stiff problem), multistep methods, also care about the past of the system (implicit like backward differentiation)    
%[t, fi] = ode113(VanDerPol_egyenlet, [0,tmax], init_xy1); %spec method for non stiff problems

%calculate the system with different initial conditions:
init_xy2=[0.1 0.1];
[t2, fi2] = ode45(VanDerPol_egyenlet, [0,tmax], init_xy2);
%[t2, fi2] = ode15s(VanDerPol_egyenlet, [0,tmax], init_xy2]); 
init_xy3=[2.5 2.5];
[t3, fi3] = ode45(VanDerPol_egyenlet, [0, tmax], init_xy3);
%[t3, fi3] = ode15s(VanDerPol_egyenlet, [0,tmax], init_xy3); 

subplot(2,2,1)
plot(fi(:,1), fi(:,2), 'b') %plot (x,y)
hold on %draw more curves in the same plot
plot(fi2(:,1), fi2(:,2), 'r')
plot(fi3(:,1), fi3(:,2), 'g')

%xlim([-4, 4]);
%ylim([-8, 8]);
xlabel('x');
ylabel('y=dx');
title(['mu=',num2str(mu)]); %write the actual mu value in the title of the plot 
hold off

subplot(2,2,2)
plot(t, fi(:,2), 'b'); %plot (t,y)
hold on % separate 'hold on' is needed for each subplot
plot(t2, fi2(:,2), 'r');
plot(t3, fi3(:,2), 'g');
%xlim([0 100]);
%ylim([-8, 8]);
xlabel('t');
ylabel('y=dx');
hold off

subplot(2,2,3)
plot(fi(:,1), t, 'b'); %plot (x,t)
hold on
plot(fi2(:,1), t2, 'r');
plot(fi3(:,1), t3, 'g');
%xlim([-4, 4]);
%ylim([0 100]);
xlabel('x');
ylabel('t');
hold off

   pause(0.01); %time delay before drawing -> animation
   waitforbuttonpress; %next iteration for pressing some key or clicking
   
end