[X,Y] = meshgrid(-3*pi-1:0.005:3*pi+1,-4:0.005:4);
Z = 0.5.*Y.*Y-cos(X)+1;

hold on;
[C,h] = contour(X,Y,Z,[2 2]);

P = zeros(7,2);
P(:,1)=linspace(-3*pi,3*pi,7);
plot(P(:,1),P(:,2),'.','Color','black');

Max_Time= 100;
b=0.1;

Xinit= -10;
Yinit= 3.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--red')

Xinit= -10;
Yinit= 2.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--magenta')

Xinit= -10;
Yinit= 2;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--blue')

Xinit= -10;
Yinit= 1.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--cyan')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xinit= 10;
Yinit= -3.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--red')

Xinit= 10;
Yinit= -2.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--magenta')

Xinit= 10;
Yinit= -2;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--blue')

Xinit= 10;
Yinit= -1.5;
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinit]);
plot(W(:,1),W(:,2),'--cyan')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bolzano shooting
b = 0.1;
eps = 1.0e-7;
Max_Time = 100;

Xinit = -10;

Yinitu = 2.5; % Upper
Yinitl = 2.0; % Lower

Yinitm = 0.0; % Ezt szeretnenk megtalalni
% Keressuk azt a kezdeti erteket ami a (pi,0) fixponthoz 
%eps kozelsegben megy el!

% Ide irjuk az eljarast

disp(strcat('Nyertes Yinit:',num2str(Yinitm)));

% Rajzoljuk is ki!
[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinitm]);
plot(W(:,1),W(:,2),'--black')


hold off;
