clear;

%Define what part of the plane to solve the problem in
mini=0;
maxi=1;
side=50;
scale=(maxi-mini)/side;
[X,Y]=meshgrid(linspace(mini,maxi,side));

%Define the vector valued function uprim according to pdf
alpha1=0.8;
alpha2=3;
alpha3=1;
epsilon=0;
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]);

%Stepsize and time.
h=0.02;
t=0:h:2;

%Allocate answer arrays
Y1=zeros(side,side,numel(t));
Y2=zeros(side,side,numel(t));

%initialize arrays wiht boundary values
Y1(:,:,1)=ones(size(Y));
Y2(:,:,1)=1-((X-0.5).^2+(Y-0.5).^2);
Y1([1,end],:,:)=0;
Y1(:,[1,end],:)=0;
Y2([1,end],:,:)=0;
Y2(:,[1,end],:)=0;

for j=2:numel(t)+1      %Timestepping
    [gx1,gy1]=gradient(Y1(:,:,j-1),scale,scale); %Calculate gradient of u1 field
    [gx2,gy2]=gradient(Y2(:,:,j-1),scale,scale); %and u2 field
    du1=del2(Y1(:,:,j-1),scale,scale); %Calculate laplace of u1
    du2=del2(Y2(:,:,j-1),scale,scale);
    for i=0:numel(X)-1   %Loop over all squares
        i1=mod(i,side)+1;       %for indexing vectors
        i2=floor(i/side)+1;

        gu1=[gx1(i1,i2);gy1(i1,i2)];  %Define gradient of u1 field in the point
        gu2=[gx2(i1,i2);gy2(i1,i2)];  %Same for u2 field
        u=[Y1(i1,i2,j-1);Y2(i1,i2,j-1)];    %Fetch values in the point at the specified time

        k1=uprim(u,gu1,gu2,du1(i1,i2),du2(i1,i2));      %Calculate next step
        k2=uprim(u+h*k1,gu1,gu2,du1(i1,i2),du2(i1,i2)); %Using simple fixed point iteration

        slope=(k1+k2)/2;    %Average them

        Y1(i1,i2,j)=Y1(i1,i2,j-1)+slope(1)*h;   %Add values to the answers
        Y2(i1,i2,j)=Y2(i1,i2,j-1)+slope(2)*h;
    end
    Y1([1,end],:,j)=1;  %Set boundaries to 0
    Y1(:,[1,end],j)=0;
    Y2([1,end],:,j)=0;
    Y2(:,[1,end],j)=0;
end
%% Record da movie

min1=min(Y1(:));
max1=max(Y1(:));
min2=min(Y1(:));
max2=max(Y2(:));

for j=1:numel(t)
    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,10)