[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

done = false;

while ~done
Yinitm = (Yinitu + Yinitl)/2;

[t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinitm]);

if W(end,1) > pi
   Yinitu=Yinitm;
else
    Yinitl=Yinitm;
end

% \exists (x_i,y_i) \in K_eps(pi) => done:=true
for i=1:size(W)
   if abs(W(i,1) - pi) < eps && abs(W(i,2)) < eps
        done=true;        
    end
end

end

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')
%}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% EZ BE FOG DOGLENI! 100 idore! MEG 1 eps re is! Az ok a linearizált egyenletben
% levo exponencialis novekedes a nyeregpont kornyeken.

% % Bolzano shooting "Naiv matematikus" verzio
% b = 0.1;
% eps = 0.01;
% Max_Time = 12;
% 
% Xinit = -10;
% 
% Yinitu = 2.5; % Upper
% Yinitl = 2.0; % Lower
% 
% done = false;
% 
% while ~done
% Yinitm = (Yinitu + Yinitl)/2;
%     
% [t,W] = ode45(@(t,y) damped_pendulum_equat(t,y,b),[0 Max_Time],[Xinit Yinitm]);
% 
% if W(end,1) > pi
%    Yinitu=Yinitm;
% else
%     Yinitl=Yinitm;
% end
% 
% % Megallasi feltetel naivan megvalasztva
% if abs(W(end,1) - pi) < eps
%     done=true;        
% end
% 
% end
% 
% 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;