clear
clc
global h

%parametri da modificare
T=7;      %tempo finale
u0=[0.5; 0.6]; %condizione iniziale Lotka-Volterra
%u0=[1; 0]; %condizione iniziale per oscillatore armonico
%u0=[2; 0; 0; -0.2]; %condizione iniziale Terra-Sole
N=100;   %numero nodi

%altri parametri
h=T/(N-1);  %passo
t(1)=0;     %inizializzazione vettore dei tempi
u(:,1)=u0;  %inizializzazione condizione iniziale

%schema
for k=1:N-1
  t(k+1)=t(k)+h;
  %u(:,k+1)=EULERO_ESPL(u(:,k));
  %u(:,k+1)=HEUN(u(:,k));
  %u(:,k+1)=RK_IV(u(:,k)); 
  %u(:,k+1)=EULERO_IMPL(u(:,k));
  %u(:,k+1)=CRANK_NIC(u(:,k));
end

%%%%%%%%%%%% PLOT %%%%%%%%%%%%%%%%%%%%%%
%le prime due componenti separate u1(t), u2(t)
subplot(2,1,1)
plot(t,u(1,:))
subplot(2,1,2)
plot(t,u(2,:))

%la curva soluzione t->(u1,u2)
figure
plot(u(1,1),u(2,1),'ro')  %pallino rosso su condizione iniziale 
axis equal
hold on
plot(u(1,:),u(2,:))

%curva soluzione (filmino)
break
figure
plot(0,0,'ko','markersize',10,'MarkerFaceColor','black')
hold on
h=plot(u(1,1),0,'ko','markersize',10,'MarkerFaceColor','red','erasemode','normal') %oscillat armon
%h=plot(u(1,1),u(1,2),'ko','markersize',10,'MarkerFaceColor','red','erasemode','normal')  %terra-sole
axis([-1 2.1 -1 1])
for k=2:N-1
  set(h,'XData',u(1,k)) %oscillat armonico
  %set(h,'XData',u(1,k),'YData',u(2,k)) %terra sole
  drawnow
end