clear all; close all; clc
digits(1000);
format longg

%%%%%%%%%%%%
%parameters%
%%%%%%%%%%%%
%lower and upper determine what points in the orbit will be plotted
%the initial condition starts at n=2
N=vpa(input('Input the value for the parameter N: '));
lower=input('Input the lower for n>1: ');
upper=input('Input the upper for n: ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%zeta%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%z=vpa(input('Input the value for the parameter zeta: '));
z=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%g%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=vpa(input('Input the value for the parameter g>0: '));
%g=.01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%color coding%%%%%%%%%%%%%%%%%%%%%%%%%
%c=input('Input 1 for color coded orbit (colored n mod 3): ');
c=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Level Set%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%level=input('Input 1 for Autonomous dP1 level set overlay: ');
level=0;
%level=1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Asymptotic Expansion%%%%%%%%%%%%%%%%%
%level set is auto. alpha=1 level set through initial condition
%expan=input('Input 1 for asymptotic expansion plot: ');
%expan=0;
expan=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Energy Readout%%%%%%%%%%%%%%%%%%%%%%%
%energ=input('Input 1 for energy readout: ');
energ=0;
%energ=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1/N Expansion%%%%%%%%%%%%%%%%%%%%%%%%
% enn=input('Input 1 for 1/N expansion: ');
%enn=0;
%enn=1;





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%computes the initial condition%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms y
fun=exp(-N*((z/2 * y.^2)+ (g/4*y.^4)));
func=(y.^2).*exp(-N*(z/2 * y.^2+ g/4*y.^4));
funct=(y.^4).*exp(-N*(z/2 * y.^2+ g/4*y.^4));
mu0=vpaintegral(fun, 0,Inf, 'RelTol', 1e-307, 'AbsTol', 0, 'MaxFunctionCalls', Inf);
mu2=vpaintegral(func, 0,Inf, 'RelTol', 1e-307, 'AbsTol', 0, 'MaxFunctionCalls', Inf);
mu4=vpaintegral(funct, 0,Inf, 'RelTol', 1e-307, 'AbsTol', 0, 'MaxFunctionCalls', Inf);
inity=vpa(mu2/mu0);
initx=vpa((mu0*mu4-mu2^2)/(mu0*mu2));


%%%%%%%%%%%%%%%%
%computes orbit%
%%%%%%%%%%%%%%%%
currentx=vpa(initx);
currenty=vpa(inity);

X=vpa(zeros(1,upper-lower+1));
Y=vpa(zeros(1,upper-lower+1));
if isequal(lower,2)
    X(1)=currentx;
    Y(1)=currenty;
end

a=vpa(-z/g);
b=vpa(1/(g*N));
n=upper; 


for k=2:n-1
    tempx=vpa(currentx);
    tempy=vpa(currenty);
    currentx=vpa(a+(b*k)/(tempx)-tempx-tempy);
    currenty=vpa(tempx);
    if k>(lower-2)
        X(k-lower+2)=currentx;
        Y(k-lower+2)=currenty;
        alpha=(vpa(k+1)/N);
        ene=vpa(currentx*currenty*(currentx+currenty+1/g)-alpha*(currentx+currenty)/g-1/g^2);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plots orbit, color coded or not%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on
if isequal(c,1)
   for k=1:upper-lower+1
        if isequal(mod(k,3),1)
            plot(X(k),Y(k), '.', 'MarkerSize',7+13*(k)/(upper-lower), 'Color',[0 0.9 0])            
        elseif isequal(mod(k,3),2)
            plot(X(k),Y(k), '.', 'MarkerSize',7+13*(k)/(upper-lower), 'Color',[1 0.9 0])
        else 
            plot(X(k),Y(k), '.', 'MarkerSize',7+13*(k)/(upper-lower), 'Color',[1 0 0])
        end
   end 
else
    for k=1:upper-lower+1
        plot(X(k),Y(k), '.', 'MarkerSize',7+13*(k)/(upper-lower), 'Color',[1 0 0])
        %change          ^  to '*' 
        % if you have points hiding behind the asymptotic expansion points
    end
end
    
%%%%%%%%%%%%%%%%%%%%%%
%plots AdP1 level set%
%%%%%%%%%%%%%%%%%%%%%%
if isequal(level,1)
    xmin=min(min(X));
    xmax=max(max(X));
    ymin=min(min(Y));
    ymax=max(max(Y));
    ex=initx;
    ey=inity;
    energy=ex*ey*(ex+ey+1/g)-(ex+ey)/g-1/(g^2);
    syms x y
    f=@(x,y)  x.*y.*(x.^1+y.^1+1/g)-(x.^1+y.^1)/g-1/(g^2)-energy;
    fimplicit(f, [xmin xmax ymin ymax], '-b');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%computes and plots truncated asymptotic expansion%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(expan,1)
   asymp=vpa(zeros(1,upper-lower+1));
   hold on
   for j=lower:upper
        s=vpa(-(j/N)*(g/4));
        z0=vpa((1-sqrt(1-48*s))/(24*s));
        z1=vpa((2*z0*(z0-1)^2)/(3*(2-z0)^4));
        z2=vpa(z0*(z0-1)*(56-294*z0+5146*z0^2-434*z0^3+126*z0^4)/(9*(2-z0)^9));
        z3=((2*z0*(z0-1))*(51744000*(z0-1)^5+257456640*(z0-1)^6+217647360*(z0-1)^7))/(362880*(2-z0)^14);         
        asymp(j-lower+1)=vpa((j/N)*(z0+(z1/(j^2))+(z2/(j^4))+0*(z3/(j^6))));
        if  (j-lower)>0
            plot(asymp(j-lower+1),asymp(j-lower), '+', 'MarkerSize',20, 'Color',[0 1 0])
        end 
      
   end
end
 

 if isequal(energ,1)  
   energy=zeros(upper-lower+1);
   for k=1:upper-lower+1
       %energy(k)=(X(k)^2)*Y(k)+(Y(k)^2)*X(k)+(1/g)*(X(k)*Y(k))-((k+lower-1)/(N*g))*(X(k)+Y(k));
       energy(k)=(X(k)^2)*Y(k)+(Y(k)^2)*X(k)+(1/g)*(X(k)*Y(k))-(1/g)*(X(k)+Y(k));
       disp('energy');
       disp(energy(k));
       disp('x');
       disp(X(k));
       disp('y');
       disp(Y(k));
   end
   
 end

