% Equazioni di Eulero
tstar = 120;
ts = 300;
k = tstar/ts;
t = 0;
y0 = [1;1;1];
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];
% RK2
y = zeros(3,ts+1);
y(:,1) = y0;
for n = 1:ts
  f1 = f(y(:,n));
  f2 = f(y(:,n)+k*f1);
  y(:,n+1) = y(:,n)+k/2*(f1+f2);
end
plot(linspace(0,tstar,ts+1),y(1,:).^2+y(2,:).^2+y(3,:).^2,'g*')
axis([0,tstar,2.8,3.3])
% trapezi
F = @(y,yn) y-yn-k/2*f(yn)-k/2*f(y);
J = @(y,yn) eye(length(y))-k/2*j(y);
y = zeros(3,ts+1);
y(:,1) = y0;
tol = 1e-14;
for n = 1:ts
  y(:,n+1) = y(:,n);
  errest = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  while (norm(errest,inf) > tol)
    y(:,n+1) = y(:,n+1)+errest;
    errest = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  end
  t = t+k;
end
hold on
plot(linspace(0,tstar,ts+1),y(1,:).^2+y(2,:).^2+y(3,:).^2,'ro')
% punto medio implicito
y0 = [1;1;1];
F = @(y,yn) y-yn-k*f((yn+y)/2);
J = @(y,yn) eye(length(y))-k*j((yn+y)/2)/2;
y = zeros(3,ts+1);
y(:,1) = y0;
tol = 1e-14;
for n = 1:ts
  y(:,n+1) = y(:,n);
  errest = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  while (norm(errest,inf) > tol)
    y(:,n+1) = y(:,n+1)+errest;
    errest = -J(y(:,n+1),y(:,n))\F(y(:,n+1),y(:,n));
  end
  t = t+k;
end
plot(linspace(0,tstar,ts+1),y(1,:).^2+y(2,:).^2+y(3,:).^2,'bx')
legend('Heun','trapezi','punto medio implicito')
xlabel('t')
ylabel('y_1(t)^2+y_2(t)^2+y_3(t)^2')