%%
%% Listing 6.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

clear;

delta = 9/7;
TOL = 1e-3; % jaką wartość $||F||^2$ uznajemy za "małą"?

% Definicja funkcji prawej strony przydatna do zabawy z polem wprost w skrypcie
% jest dopuszczalna tylko w Octave: MATLAB wymaga zdefiniowania tej funkcji 
% w osobnym pliku
function [f1, f2, f3] = pole(u1, u2, u3, delta)
f1 = 1 - u1 - delta*(u2.^2+u3.^2);
f2 = -u2 + delta*u1.*u2;
f3 = -u3 + delta*u1.*u3;
end

% szukamy punktów stacjonarnych w prostokącie $(u_2,u_3)\in [-A,A]\times [-A,A]$...
N = 30; A = 0.6;
u2 = linspace(-A,A,N);
u3 = linspace(-A,A,N);
[U2, U3] = meshgrid(u2,u3);

% ...i dla u1 = [0.7..0.8]
for u1 = linspace(0.7,1.2,50)
	[F1 F2 F3] = pole(u1, U2, U3, delta);
	G = double(F1.^2 + F2.^2 + F3.^2 < TOL);
	% czy znaleźliśmy choć jeden punkt, w którym F jest małe?
	if max(G(:)) > 0
		fprintf(stderr,'Mala norma dla u1 = %g',u1); % MATLAB nie zna stderr...
		[i, j] = find(G);
		fprintf(stderr,' i np. (u2,u3) = (%5.3g,%5.3g)\n', ...
					U2(i(1),j(1)), U3(i(1),j(1)));
		mesh(U2,U3,G); 
		title(['u1=', num2str(u1)]);
		pause(1);
	end
end

