disp('Robertson example');

F = @(Y,t) [ -0.04*Y(1) + 1e4*Y(2)*Y(3); ...
							0.04*Y(1) - 1e4*Y(2)*Y(3) - 3e7*Y(2)^2; ...
							3e7*Y(2)^2];
T = 1; N = 100;  
x0 = [1 0 0]';

g = @(X)F(X,t); opts.TolX = opts.TolFun = 1e-12;
Ystationary = fsolve(g,[1 0.0000365 0]',opts) 

% jesteśmy więc b. blisko punktu stacjonarnego, 
% gdzie wartości własne to 0 (skąd degeneracja Jakobianu i problem w fsolve) 
% oraz dwie ujemne: -0.405 oraz -2189.6

t = linspace(0,T,N+1);
Y = lsode(F,x0,t);
plot(t,Y);

xe = zeros(length(x0),size(t)); xe(:,1) = x0; 	% ... dla Eulera
xie = xe; 							% ... dla iEulera

h = 0.001; 

n = (T/h)/N; FF = @(t,X) F(X,t);
for i = 1:N
	% wykonujemy JEDEN krok schematu...
	xe(:,i+1) 	= euler(FF, xe(:,i), t(i), t(i+1), n );
	xie(:,i+1) 	= impliciteuler(FF, xie(:,i), t(i), t(i+1), n );
end

plot(t,Y(:,1),'g');
hold on;
plot(t,xe(1,:),'r');
plot(t,xie(1,:),'b');
hold off;
legend('lsode','Euler','iEuler');

