% Solution of a system of 2 ODes
% Numerical method: modified Euler (Runge-Kutta, order 2)

clearvars;
close all;

% x' = a x + b y
% y' = c x + d y

% Parameters

a = 0;
b = 1;
c = -1;
d = 0;

% Functions

F=@(x,y) a*x+b*y;
G=@(x,y) c*x+d*y;

% Time definitions

Tmax = 100;
dt = 0.01;
t = 0:dt:Tmax;

x = zeros(1,length(t));
y = zeros(1,length(t));

x(1) = 1;
y(1) = 0;

for j=1:length(t)-1
    k1x = F(x(j),y(j));
    k1y = G(x(j),y(j));
    ax = x(j)+k1x*dt;
    ay = y(j)+k1y*dt;
    k2x = F(ax,ay);
    k2y = G(ax,ay);
    x(j+1)=x(j)+(k1x+k2x)*dt/2;
    y(j+1)=y(j)+(k1y+k2y)*dt/2;    
end

figure
hold on
plot(t,x,'-b','linewidth',2);
plot(t,y,'-r','linewidth',2);
axis([-0.2 Tmax -3.2 3.2]);
set(gca,'fontsize',24);
xlabel('t');
ylabel('');
legend('x','y');

figure
hold on
plot(x,y,'-b','linewidth',2);
plot([0 0],[-100 100],'--','Color',[.7 .7 .7],'linewidth',1);
plot([-100 100],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
axis([-2 2 -2 2]);
set(gca,'fontsize',24);
xlabel('x');
ylabel('y');

figure
for j=1:length(t)   
    plot([0 0],[-100 100],'--','Color',[.7 .7 .7],'linewidth',1);
    hold on
    plot([-100 100],[0 0],'--','Color',[.7 .7 .7],'linewidth',1);
    plot(x(1:j),y(1:j),'-b','linewidth',2);
    plot(x(j),y(j),'or','linewidth',3);
    axis([-2 2 -2 2]);
    set(gca,'fontsize',24);
    xlabel('x');
    ylabel('y');        
    pause(0.00001);
    hold off
end






