%3. gyakorlat: Van der Pol oszcillátor:
%ddx-mu*(1-x^2)*dx+x=0 ->
%dx=y -> ddx=dy -> dy=mu*(1-x^2)*y-x
close all;

tmax=100;%100;%1000;% 
%kezdeti értékek
init_xy1=[0.5 0.5];
init_xy2=[0.1 0.1];
init_xy3=[2 2];

for mu = 0:0.1:10%5%mu = 0:50:500% %különbözõ mu paraméterekre számolom ki az egyenletet, az ode45 belassul nagy mu értékeknél
%
%Explicit Euler módszerrel:
% mu=8-nél elszáll, dt=0.05-re 
dt=0.01; %0.05 %idõbeli lépésköz
t=0:dt:tmax;
t2=t;
t3=t;
fi=zeros(length(t),2);
fi(1,:)=init_xy1;
fi2=zeros(length(t),2);
fi2(1,:)=init_xy2;
fi3=zeros(length(t),2);
fi3(1,:)=init_xy3;

for i=2:1:length(t)
    fi(i,:)= fi(i-1,:)+dt*[fi(i-1,2), mu*(1-fi(i-1,1)^2)*fi(i-1,2)-fi(i-1,1)];
    fi2(i,:)= fi2(i-1,:)+dt*[fi2(i-1,2), mu*(1-fi2(i-1,1)^2)*fi2(i-1,2)-fi2(i-1,1)];
    fi3(i,:)= fi3(i-1,:)+dt*[fi3(i-1,2), mu*(1-fi3(i-1,1)^2)*fi3(i-1,2)-fi3(i-1,1)];
end
%}
    
%{
%MatLab differenciál egyenlet megoldókkal (ode) :
    VanDerPol_egyenlet=(@(t, xy)[xy(2); mu*(1-xy(1)^2)*xy(2)-xy(1)]); 

[t, fi] = ode45(VanDerPol_egyenlet, [0,tmax], init_xy1); %általános módszer
%[t, fi] = ode15s(VanDerPol_egyenlet, [0,tmax], init_xy1); %az ode15s tudja kezelni a nagyon ha a rendszerben nagyon gyors és nagyon lassú változások is vannak ("stiff problem")  
% ezek a módszerek a rendszer múltját is figyelik, változtatják a lépésközt ennek alapján 

%a rendszer kiszámolása különbözõ kezdeti értékekkel
[t2, fi2] = ode45(VanDerPol_egyenlet, [0,tmax], init_xy2);
%[t2, fi2] = ode15s(VanDerPol_egyenlet, [0,tmax], init_xy2); 

[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), 'bx-') %plot (x,y)
hold on %a görbék egymásra rajzolása
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)]); %az aktuális mu érték a plot címében
hold off

subplot(2,2,2)
plot(t, fi(:,2), 'bx-'); %plot (t,y)
hold on % külön 'hold on' kell minden subplotra
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, 'bx-'); %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); %késleltetés kirajzolás elött -> animáció
   waitforbuttonpress; %következõ iteráció kattintásra
   
end