% Equazioni di Eulero
tstar = 120;
ts = 300;
k = tstar/ts;
t = 0;
I = [2;1;2/3];
f = @(y) [(I(2)-I(3))/(I(2)*I(3))*y(2)*y(3);...
	  (I(3)-I(1))/(I(3)*I(1))*y(3)*y(1);...
	  (I(1)-I(2))/(I(1)*I(2))*y(1)*y(2)];
j = @(y) [0,(I(2)-I(3))/(I(2)*I(3))*y(3),(I(2)-I(3))/(I(2)*I(3))*y(2);...
	  (I(3)-I(1))/(I(3)*I(1))*y(3),0,(I(3)-I(1))/(I(3)*I(1))*y(1);...
	  (I(1)-I(2))/(I(1)*I(2))*y(2),(I(1)-I(2))/(I(1)*I(2))*y(1),0];
y0 = [1;1;1];
% punto medio implicito
F = @(y,yn) y-yn-k*f((yn+y)/2);
J = @(y,yn) eye(length(y))-k*j((yn+y)/2)/2;
% trapezi
%F = @(y,yn) y-yn-k*(f(yn)+f(y))/2;
%J = @(y,yn) eye(length(y))-k*j(y)/2;
y = zeros(length(y0),ts+1);
y(:,1) = y0;
tol = k^2/100;
for n = 1:ts
  y(:,n+1) = y(:,n);
  delta = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  while (norm(delta,inf) > tol)
    y(:,n+1) = y(:,n+1)+delta;
    delta = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  end
end
plot(linspace(0,tstar,ts+1),y(1,:).^2+y(2,:).^2+y(3,:).^2,'bx')
xlabel('t')
ylabel('y_1(t)^2+y_2(t)^2+y_3(t)^2')
