N = 38;
figure;
plot([0 N^2 N^2 0 0],[N^2 N^2 0 0 N^2]);hold on;
IntegersXs = [];
IntegersYs = [];
for i=1:2:N,
    
    if isprime(i) && isprime(N - i);
    
    plot([0 i i 0 0],[i i 0 0 i],'r-');hold on;
    plot([i i+(N-i)^2 i+(N-i)^2 i i],[i+(N-i)^2 i+(N-i)^2 i i i+(N-i)^2],'r-','LineWidth',3);hold on;
    plot(N^2 - [0 i i 0 0],N^2 - [i i 0 0 i],'r-');hold on;
    plot(N^2 - [i i+(N-i)^2 i+(N-i)^2 i i],N^2 - [i+(N-i)^2 i+(N-i)^2 i i i+(N-i)^2],'r-','LineWidth',3);hold on;
    else
        plot([0 i i 0 0],[i i 0 0 i],'b-');hold on;
        plot([i i+(N-i)^2 i+(N-i)^2 i i],[i+(N-i)^2 i+(N-i)^2 i i i+(N-i)^2],'b-','LineWidth',3);hold on;
        plot(N^2 - [0 i i 0 0],N^2 - [i i 0 0 i],'b-');hold on;
        plot(N^2 - [i i+(N-i)^2 i+(N-i)^2 i i],N^2 - [i+(N-i)^2 i+(N-i)^2 i i i+(N-i)^2],'b-','LineWidth',3);hold on;
    end

    text(N + 3, i+(N-i)^2 + 8,strcat(num2str(i),'+',num2str(N-i),'^2'));
    Xs = 1:(N-i)^2;
    Ys = (N-i)^2./Xs;
    XsRot = i + (N-i)^2 - Ys;
    YsRot = i + (N-i)^2 - Xs;
    plot(XsRot,YsRot,'LineWidth',1);hold on;
    plot(N^2 - XsRot,N^2 - YsRot,'LineWidth',1);hold on;
    IntegersXs(end+1:end+length(XsRot((mod(Ys,1) < 0.001)))) = XsRot((mod(Ys,1) < 0.001));
    IntegersYs(end+1:end+length(YsRot((mod(Ys,1) < 0.001)))) = YsRot((mod(Ys,1) < 0.001));
%plot([(N/2 - i)^2 (N/2 - i)^2+(N/2 +i)^2 (N/2 - i)^2+(N/2 +i)^2 (N/2 - i)^2 (N/2 - i)^2 ], [(N/2 - i)^2 + (N/2 +i)^2 (N/2 - i)^2 + (N/2 +i)^2 0 0 (N/2 - i) + (N/2 +i)^2 ]);hold on;
end


    plot(IntegersXs,IntegersYs,'Marker','o','MarkerSize',7,'LineStyle','none','MarkerFaceColor',[0 1 0]);
    plot(N^2 - IntegersXs,N^2 - IntegersYs,'Marker','o','MarkerSize',7,'LineStyle','none','MarkerFaceColor',[0 1 0]);


set(gca,'XLim',[0 N^2],'YLim',[0 N^2]);