%Runge-Kutta version

clear;

side=50;
[X,Y]=meshgrid(linspace(0,1,side));

x=ones(size(Y));
y=1-((X-0.5).^2+(Y-0.5).^2);

alpha1=0.8;
alpha2=1.2;
alpha3=0.5;
epsilon=0.01;
Beta1=-[0.1;0.1];
Beta2=-[0.1;0.1];
uprim=@(u,gu1,gu2,du1,du2)([
    -alpha1*u(1)*u(2)-dot(Beta1,gu1)+epsilon*du1;
    alpha2*u(1)*u(2)-alpha3*u(2)-dot(Beta2,gu2)+epsilon*du2]);
h=0.05;
t=0:h:2;

Y1=zeros(side,side,numel(t));
Y2=zeros(side,side,numel(t));

Y1(:,:,1)=x;
Y2(:,:,1)=y;
for j=2:numel(t)+1
    [gx1,gy1]=gradient(Y1(:,:,j-1),h,h);
    [gx2,gy2]=gradient(Y2(:,:,j-1),h,h);
    du1=del2(Y1(:,:,j-1),h,h);
    du2=del2(Y2(:,:,j-1),h,h);
    k1=zeros(side,side,2);
    k2=zeros(side,side,2);
    k3=zeros(side,side,2);
    for i=0:numel(X)-1
        i1=mod(i,side)+1;
        i2=floor(i/side)+1;

        gu1=[gx1(i1,i2);gy1(i1,i2)];
        gu2=[gx2(i1,i2);gy2(i1,i2)];
        u=[Y1(i1,i2,j-1);Y2(i1,i2,j-1)];

        k1(i1,i2,:)=uprim(u,gu1,gu2,du1(i1,i2),du2(i1,i2));
    end
    Yk1=cat(3,Y1(:,:,j-1),Y2(:,:,j-1))+h*k1/2;
    [gx1,gy1]=gradient(Yk1(:,:,1),h,h);
    [gx2,gy2]=gradient(Yk1(:,:,2),h,h);
    du1=del2(Yk1(:,:,1),h,h);
    du2=del2(Yk1(:,:,2),h,h);
    for i=0:numel(X)-1
        i1=mod(i,side)+1;
        i2=floor(i/side)+1;

        gu1=[gx1(i1,i2);gy1(i1,i2)];
        gu2=[gx2(i1,i2);gy2(i1,i2)];
        u=[Yk1(i1,i2,1);Yk1(i1,i2,2)];

        k2(i1,i2,:)=uprim(Yk1,gu1,gu2,du1(i1,i2),du2(i1,i2));
    end
    Yk2=cat(3,Y1(:,:,j-1),Y2(:,:,j-1))+h*k2/2;
    [gx1,gy1]=gradient(Yk2(:,:,1),h,h);
    [gx2,gy2]=gradient(Yk2(:,:,2),h,h);
    du1=del2(Yk2(:,:,1),h,h);
    du2=del2(Yk2(:,:,2),h,h);
    for i=0:numel(X)-1
        i1=mod(i,side)+1;
        i2=floor(i/side)+1;

        gu1=[gx1(i1,i2);gy1(i1,i2)];
        gu2=[gx2(i1,i2);gy2(i1,i2)];
        u=[Yk2(i1,i2,1);Yk2(i1,i2,2)];
        k3(i1,i2,:)=uprim(u,gu1,gu2,du1(i1,i2),du2(i1,i2));
    end
    Yk3=cat(3,Y1(:,:,j-1),Y2(:,:,j-1))+h*k2;
    [gx1,gy1]=gradient(Yk3(:,:,1),h,h);
    [gx2,gy2]=gradient(Yk3(:,:,2),h,h);
    du1=del2(Yk3(:,:,1),h,h);
    du2=del2(Yk3(:,:,2),h,h);
    for i=0:numel(X)-1
        i1=mod(i,side)+1;
        i2=floor(i/side)+1;

        gu1=[gx1(i1,i2);gy1(i1,i2)];
        gu2=[gx2(i1,i2);gy2(i1,i2)];
        u=[Yk3(i1,i2,1);Yk3(i1,i2,2)];

        k4(i1,i2,:)=uprim(u,gu1,gu2,du1(i1,i2),du2(i1,i2));
    end
    slope=(k1+2*k2+2*k3+k4)/6;

    Y1(:,:,j)=Y1(:,:,j-1)+slope(:,:,1)*h;
    Y2(:,:,j)=Y2(:,:,j-1)+slope(:,:,2)*h;
end
%% Record da movie
min1=min(Y1(:));
min2=min(Y2(:));
max1=max(Y1(:));
max2=max(Y2(:));

for j=1:numel(t)
    clf
    subplot(2,1,1)
    imagesc(Y1(:,:,j))
    colorbar
    axis off
    subplot(2,1,2)
    imagesc(Y2(:,:,j))
    colorbar
    axis off
    M(j)=getframe(figure(1));
end

%% PLAY!!!
movie(M,3,5)