function [dz,dLambda,du,da_oa] = ForcedElements1(P,Lambda,mu,z,u,OrdersAll,MaxPower)
% ForcedElements1: calculate the forced eccentricities and mean longitudes
% by a disturbing function expansion.

%Difference from ForcedElements:
% rather than writing down a series of terms, use a general expression for
% the disturbing function terms contributions and sum them up.

%Inputs: P - orbital periods, Lambda - mean longitudes, mu - planets/star
%mass ratio, z - free eccentricities (if they are not given, zero values
%are assumed), u - free inclinations (=I*exp(i*Omega) (if they are not given, zero values
%are assumed). Labmda, z and u can be given as columns of multiple values.
%OrdersAll: the number of j values taken in the expansion. If given an
%empty value [], the code will decide by itself what is the required
%number. MaxPower: the maximal (joint) power of eccentricities and
%inclinations taken in the expansion. Default is 4.

% Yair Judkovsky, 6.12.2020. Coded with lots of blood, swear and ink :)

if ~exist('z','var')
    z = zeros(size(P));
end

if ~exist('u','var')
    u = zeros(size(P));
end

if ~exist('OrdersAll','var')
    OrdersAll = [];
end

if ~exist('MaxPower','var')
    MaxPower = 4;
end

%mean motions
n = 2*pi./P;

%absolute values and angles of eccentricities
e = abs(z);
pom = angle(z);

%absolute values and angles of inclinations
I = abs(u);
Om = angle(u);

%allocate space for 3d matrices that will hold the values of dLambda and dz
%according to these dimensios: orders, data points, number of planets
dLambdaMatrix = zeros(size(z,1),length(P));
dzMatrix = zeros(size(z,1),length(P));
duMatrix = zeros(size(z,1),length(P));
da_oaMatrix = zeros(size(z,1),length(P));

%for each pair of planets, calculate the 2nd order forced
%eccentricity. secular terms are ignored.
for j1 = 1:length(P)
    for j2 = (j1+1):length(P)
        
        %ratio of semi-major axes
        alph = (n(j2)/n(j1))^(2/3)*((1+mu(j1))/(1+mu(j2)))^(1/3);
        
        %write the orbital elements for this specific pair
        Lambda1 = Lambda(:,j1);
        Lambda2 = Lambda(:,j2);
        
        %write the complex eccentricities, magnitude and angle of the
        %currently analyzed pair        
        e1 = e(:,j1);
        e2 = e(:,j2);
        pom1 = pom(:,j1);
        pom2 = pom(:,j2);
        
        %write the complex inclinations, magnitude and angle of the
        %currently analyzed pair
        I1 = I(:,j1);
        I2 = I(:,j2);
        Om1 = Om(:,j1);
        Om2 = Om(:,j2);
        
        %translate from I to s - as defined in the disturbing function
        %expansion (Solar System Dynamics)
        s1 = sin(I1/2);
        s2 = sin(I2/2);
        
        %write the relative masses and the mean motions of the currently analyzed pair
        mu1 = mu(j1);
        mu2 = mu(j2);
        n1 = n(j1);
        n2 = n(j2);      
                
        %set number of orders, if not given, by the proximity to closest
        %MMRs
        if isempty(OrdersAll)
            J = round((1:MaxPower)*P(j2)/(P(j2)-P(j1)));
            Orders = max([J 6])+1;
        else
            Orders = OrdersAll;
        end
        
        %prepare the structure Terms, which will include the table elements
        Terms = [];
        clear AllTerms
        
        %here, instead of the loop on j, make a loop on the different
        %terms, to be continued
        for j = -Orders:Orders
            
            
            %% CalculateSMAFunctions;
            
            %Calculate Laplace coefficients and their semi-major axis functions - Solar
            %System Dynamics (Murray & Dermott 1999), Appendix B
            
            for DummyIdx = 1:1
                
                
            
%                 Aj = LaplaceCoeffs(alph,0.5,j,0);
%                 DAj = LaplaceCoeffs(alph,0.5,j,1);
%                 D2Aj = LaplaceCoeffs(alph,0.5,j,2);
%                 D3Aj = LaplaceCoeffs(alph,0.5,j,3);
%                 D4Aj = LaplaceCoeffs(alph,0.5,j,4);
%                 D5Aj = LaplaceCoeffs(alph,0.5,j,5);
%                 
%                 Ajp1 = LaplaceCoeffs(alph,0.5,j+1,0);
%                 DAjp1 = LaplaceCoeffs(alph,0.5,j+1,1);
%                 D2Ajp1 = LaplaceCoeffs(alph,0.5,j+1,2);
%                 D3Ajp1 = LaplaceCoeffs(alph,0.5,j+1,3);
%                 D4Ajp1 = LaplaceCoeffs(alph,0.5,j+1,4);
%                 D5Ajp1 = LaplaceCoeffs(alph,0.5,j+1,5);
%                 
%                 Ajp2 = LaplaceCoeffs(alph,0.5,j+2,0);
%                 DAjp2 = LaplaceCoeffs(alph,0.5,j+2,1);
%                 D2Ajp2 = LaplaceCoeffs(alph,0.5,j+2,2);
%                 D3Ajp2 = LaplaceCoeffs(alph,0.5,j+2,3);
%                 D4Ajp2 = LaplaceCoeffs(alph,0.5,j+2,4);
%                 D5Ajp2 = LaplaceCoeffs(alph,0.5,j+2,5);
%                 
%                 Ajm1 = LaplaceCoeffs(alph,0.5,j-1,0);
%                 DAjm1 = LaplaceCoeffs(alph,0.5,j-1,1);
%                 D2Ajm1 = LaplaceCoeffs(alph,0.5,j-1,2);
%                 D3Ajm1 = LaplaceCoeffs(alph,0.5,j-1,3);
%                 D4Ajm1 = LaplaceCoeffs(alph,0.5,j-1,4);
%                 D5Ajm1 = LaplaceCoeffs(alph,0.5,j-1,5);
%                 
%                 Ajm2 = LaplaceCoeffs(alph,0.5,j-2,0);
%                 DAjm2 = LaplaceCoeffs(alph,0.5,j-2,1);
%                 D2Ajm2 = LaplaceCoeffs(alph,0.5,j-2,2);
%                 D3Ajm2 = LaplaceCoeffs(alph,0.5,j-2,3);
%                 D4Ajm2 = LaplaceCoeffs(alph,0.5,j-2,4);
%                 D5Ajm2 = LaplaceCoeffs(alph,0.5,j-2,5);
%                 
%                 Ajm3 = LaplaceCoeffs(alph,0.5,j-3,0);
%                 DAjm3 = LaplaceCoeffs(alph,0.5,j-3,1);
%                 D2Ajm3 = LaplaceCoeffs(alph,0.5,j-3,2);
%                 D3Ajm3 = LaplaceCoeffs(alph,0.5,j-3,3);
%                 D4Ajm3 = LaplaceCoeffs(alph,0.5,j-3,4);
%                 D5Ajm3 = LaplaceCoeffs(alph,0.5,j-3,5);
%                 
%                 Ajm4 = LaplaceCoeffs(alph,0.5,j-4,0);
%                 DAjm4 = LaplaceCoeffs(alph,0.5,j-4,1);
%                 D2Ajm4 = LaplaceCoeffs(alph,0.5,j-4,2);
%                 D3Ajm4 = LaplaceCoeffs(alph,0.5,j-4,3);
%                 D4Ajm4 = LaplaceCoeffs(alph,0.5,j-4,4);
%                 D5Ajm4 = LaplaceCoeffs(alph,0.5,j-4,5);
%                 
%                 
%                 Bj = LaplaceCoeffs(alph,1.5,j,0);
%                 DBj = LaplaceCoeffs(alph,1.5,j,1);
%                 D2Bj = LaplaceCoeffs(alph,1.5,j,2);
%                 D3Bj = LaplaceCoeffs(alph,1.5,j,3);
%                 
%                 Bjm1 = LaplaceCoeffs(alph,1.5,j-1,0);
%                 DBjm1 = LaplaceCoeffs(alph,1.5,j-1,1);
%                 D2Bjm1 = LaplaceCoeffs(alph,1.5,j-1,2);
%                 D3Bjm1 = LaplaceCoeffs(alph,1.5,j-1,3);
%                 
%                 Bjm2 = LaplaceCoeffs(alph,1.5,j-2,0);
%                 DBjm2 = LaplaceCoeffs(alph,1.5,j-2,1);
%                 D2Bjm2 = LaplaceCoeffs(alph,1.5,j-2,2);
%                 D3Bjm2 = LaplaceCoeffs(alph,1.5,j-2,3);
%                 
%                 Bjm3 = LaplaceCoeffs(alph,1.5,j-3,0);
%                 DBjm3 = LaplaceCoeffs(alph,1.5,j-3,1);
%                 D2Bjm3 = LaplaceCoeffs(alph,1.5,j-3,2);
%                 D3Bjm3 = LaplaceCoeffs(alph,1.5,j-3,3);
%                 
%                 
%                 
%                 Bjp1 = LaplaceCoeffs(alph,1.5,j+1,0);
%                 DBjp1 = LaplaceCoeffs(alph,1.5,j+1,1);
%                 D2Bjp1 = LaplaceCoeffs(alph,1.5,j+1,2);
%                 D3Bjp1 = LaplaceCoeffs(alph,1.5,j+1,3);
%                 
%                 Bjp2 = LaplaceCoeffs(alph,1.5,j+2,0);
%                 DBjp2 = LaplaceCoeffs(alph,1.5,j+2,1);
%                 D2Bjp2 = LaplaceCoeffs(alph,1.5,j+2,2);
%                 D3Bjp2 = LaplaceCoeffs(alph,1.5,j+2,3);
%                 
%                 Cj = LaplaceCoeffs(alph,2.5,j,0);
%                 DCj = LaplaceCoeffs(alph,2.5,j,1);
%                 
%                 Cjm2 = LaplaceCoeffs(alph,2.5,j-2,0);
%                 DCjm2 = LaplaceCoeffs(alph,2.5,j-2,1);
%                 
%                 Cjp2 = LaplaceCoeffs(alph,2.5,j+2,0);
%                 DCjp2 = LaplaceCoeffs(alph,2.5,j+2,1);
                
                
            
                sVec = [0.5*ones(1,42) 1.5*ones(1,24) 2.5*ones(1,6)];
                jVec = [j*ones(1,6) (j+1)*ones(1,6) (j+2)*ones(1,6) (j-1)*ones(1,6)  (j-2)*ones(1,6) (j-3)*ones(1,6) (j-4)*ones(1,6) (j)*ones(1,4) (j-1)*ones(1,4) (j-2)*ones(1,4) (j-3)*ones(1,4) (j+1)*ones(1,4) (j+2)*ones(1,4) j j j-2 j-2 j+2 j+2];
                DVec = [0:5 0:5 0:5 0:5 0:5 0:5 0:5 0:3 0:3 0:3 0:3 0:3 0:3 0:1 0:1 0:1];
                CoeffsVec = LaplaceCoeffs(alph,sVec,jVec,DVec);
                Aj = CoeffsVec(1);
                DAj = CoeffsVec(2);
                D2Aj = CoeffsVec(3);
                D3Aj = CoeffsVec(4);
                D4Aj  = CoeffsVec(5);
                D5Aj  = CoeffsVec(6);
                
                Ajp1 = CoeffsVec(7);
                DAjp1  = CoeffsVec(8);
                D2Ajp1  = CoeffsVec(9);
                D3Ajp1  = CoeffsVec(10);
                D4Ajp1  = CoeffsVec(11);
                D5Ajp1  = CoeffsVec(12);
                
                Ajp2  = CoeffsVec(13);
                DAjp2 = CoeffsVec(14);
                D2Ajp2  = CoeffsVec(15);
                D3Ajp2  = CoeffsVec(16);
                D4Ajp2  = CoeffsVec(17);
                D5Ajp2  = CoeffsVec(18);
                
                Ajm1  = CoeffsVec(19);
                DAjm1  = CoeffsVec(20);
                D2Ajm1 = CoeffsVec(21);
                D3Ajm1  = CoeffsVec(22);
                D4Ajm1  = CoeffsVec(23);
                D5Ajm1  = CoeffsVec(24);
                
                Ajm2  = CoeffsVec(25);
                DAjm2  = CoeffsVec(26);
                D2Ajm2  = CoeffsVec(27);
                D3Ajm2  = CoeffsVec(28);
                D4Ajm2  = CoeffsVec(29);
                D5Ajm2  = CoeffsVec(30);
                
                Ajm3  = CoeffsVec(31);
                DAjm3  = CoeffsVec(32);
                D2Ajm3  = CoeffsVec(33);
                D3Ajm3  = CoeffsVec(34);
                D4Ajm3  = CoeffsVec(35);
                D5Ajm3  = CoeffsVec(36);
                
                Ajm4  = CoeffsVec(37);
                DAjm4  = CoeffsVec(38);
                D2Ajm4  = CoeffsVec(39);
                D3Ajm4  = CoeffsVec(40);
                D4Ajm4  = CoeffsVec(41);
                D5Ajm4  = CoeffsVec(42);
                
                
                Bj  = CoeffsVec(43);
                DBj  = CoeffsVec(44);
                D2Bj  = CoeffsVec(45);
                D3Bj = CoeffsVec(46);
                
                Bjm1  = CoeffsVec(47);
                DBjm1  = CoeffsVec(48);
                D2Bjm1  = CoeffsVec(49);
                D3Bjm1  = CoeffsVec(50);
                
                Bjm2  = CoeffsVec(51);
                DBjm2  = CoeffsVec(52);
                D2Bjm2  = CoeffsVec(53);
                D3Bjm2  = CoeffsVec(54);
                
                Bjm3  = CoeffsVec(55);
                DBjm3  = CoeffsVec(56);
                D2Bjm3  = CoeffsVec(57);
                D3Bjm3  = CoeffsVec(58);
                
                
                
                Bjp1  = CoeffsVec(59);
                DBjp1  = CoeffsVec(60);
                D2Bjp1  = CoeffsVec(61);
                D3Bjp1  = CoeffsVec(62);
                
                Bjp2 = CoeffsVec(63);
                DBjp2  = CoeffsVec(64);
                D2Bjp2  = CoeffsVec(65);
                D3Bjp2  = CoeffsVec(66);
                
                Cj  = CoeffsVec(67);
                DCj  = CoeffsVec(68);
                
                Cjm2 = CoeffsVec(69);
                DCjm2 = CoeffsVec(70);
                
                Cjp2 = CoeffsVec(71);
                DCjp2 = CoeffsVec(72);
                
            
            
            
            
            %0th order semi-major axis functions
            f1 = 0.5*Aj;
            Df1 = 0.5*DAj;
            
            f2 =   1/8*(-4*j^2*Aj  +2*alph*DAj             +alph^2*D2Aj);
            Df2 = 1/8*(-4*j^2*DAj+2*DAj+2*alph*D2Aj+2*alph*D2Aj+alph^2*D3Aj);
            
            f3 =   -0.25*alph*Bjm1                   -0.25*alph*Bjp1;
            Df3 = -0.25*Bjm1-0.25*alph*DBjm1-0.25*Bjp1-0.25*alph*DBjp1;
            
            f4 =   1/128*((-9*j^2+16*j^4)*Aj  +(-8*j^2)*alph*DAj            +(-8*j^2)*alph^2*D2Aj                      +4*alph^3*D3Aj                          +alph^4*D4Aj);
            Df4 = 1/128*((-9*j^2+16*j^4)*DAj+(-8*j^2)*(DAj+alph*D2Aj)+(-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+4*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4+D5Aj));
            
            f5 =   1/32*((16*j^4)*Aj   +(4-16*j^2)*alph*DAj             +(14-8*j^2)*alph^2*D2Aj                       +(8)*alph^3*D3Aj                           +(1)*alph^4*D4Aj);
            Df5 = 1/32*((16*j^4)*DAj +(4-16*j^2)*(DAj+alph*D2Aj) +(14-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj) +(8)*(3*alph^2*D3Aj+alph^3*D4Aj) +(1)*(4*alph^3*D4Aj+alph^4*D5Aj));
            
            f6 =   1/128*((-17*j^2+16*j^4)*Aj  +(24-24*j^2)*alph*DAj            +(36-8*j^2)*alph^2*D2Aj                      +(12)*alph^3*D3Aj                          +(1)*alph^4*D4Aj);
            Df6 = 1/128*((-17*j^2+16*j^4)*DAj+(24-24*j^2)*(DAj+alph*D2Aj)+(36-8*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+(12)*(3*alph^2*D3Aj+alph^3*D4Aj)+(1)*(4*alph^3*D4Aj+alph^4*D5Aj));
            
            f7 = 1/16*((-2+4*j^2)*alph*Bjm1                +(-4)*alph^2*DBjm1                           +(-1)*alph^3*D2Bjm1)...
                +1/16*((-2+4*j^2)*alph*Bjp1                 +(-4)*alph^2*DBjp1                            +(-1)*alph^3*D2Bjp1);
            Df7 = 1/16*((-2+4*j^2)*(Bjm1+alph*DBjm1)+(-4)*(2*alph*DBjm1+alph^2*D2Bjm1)+(-1)*(3*alph^2*D2Bjm1+alph^3*D3Bjm1))...
                +1/16*((-2+4*j^2)*(Bjp1+alph*DBjp1)  +(-4)*(2*alph*DBjp1+alph^2*D2Bjp1)  +(-1)*(3*alph^2*D2Bjp1+alph^3*D3Bjp1));
            
            f8 =   alph^2*(3/16*Cjm2+0.75*Cj+3/16*Cjp2);
            Df8 = 2*alph*(3/16*Cjm2+0.75*Cj+3/16*Cjp2)+alph^2*(3/16*DCjm2+0.75*DCj+3/16*DCjp2);
            
            f9 =   0.25*alph*(Bjm1+Bjp1)                           +3/8*alph^2*Cjm2                         +15/4*alph^2*Cj                    +3/8*alph^2*Cjp2;
            Df9 = 0.25*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))+3/8*(2*alph*Cjm2+alph^2*DCjm2)+15/4*(2*alph*Cj+alph^2*DCj)+3/8*(2*alph*Cjp2+alph^2*DCjp2);
            
            f10 =   0.25*((2+6*j+4*j^2)*Ajp1   -2*alph*DAjp1              -alph^2*D2Ajp1);
            Df10 = 0.25*((2+6*j+4*j^2)*DAjp1-2*(DAjp1+alph*D2Ajp1)-(2*alph*D2Ajp1+alph^2*D3Ajp1));
            
            f11 =   1/32*((-6*j-26*j^2-36*j^3-16*j^4)*Ajp1  +(6*j+12*j^2)*alph*DAjp1               +(-4+7*j)*alph^2*D2Ajp1                          -6*alph^3*D3Ajp1                             -alph^4*D4Ajp1);
            Df11 = 1/32*((-6*j-26*j^2-36*j^3-16*j^4)*DAjp1+(6*j+12*j^2)*(DAjp1+alph*D2Ajp1)+(-4+7*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)-6*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1));
            
            f12 =   1/32*((4+2*j-22*j^2-36*j^3-16*j^4)*Ajp1  +(-4+22*j+20*j^2)*alph*DAjp1               +(-22+7*j+8*j^2)*alph^2*D2Ajp1                          -10*alph^3*D3Ajp1                              -alph^4*D4Ajp1);
            Df12 = 1/32*((4+2*j-22*j^2-36*j^3-16*j^4)*DAjp1+(-4+22*j+20*j^2)*(DAjp1+alph*D2Ajp1)+(-22+7*j+8*j^2)*(2*alph*D2Ajp1+alph^2*D3Ajp1)-10*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1));
            
            f13   = 1/8*((-6*j-4*j^2)*alph*(Bj+Bjp2)                      +4*alph^2*(DBj+DBjp2)                                    +alph^3*(D2Bj+D2Bjp2));
            Df13 = 1/8*((-6*j-4*j^2)*((Bj+Bjp2)+alph*(DBj+DBjp2))+4*(2*alph*(DBj+DBjp2)+alph^2*(D2Bj+D2Bjp2))+(3*alph^2*(D2Bj+D2Bjp2)+alph^3*(D3Bj+D3Bjp2)));
            
            f14 = alph*Bjp1;
            Df14 = Bjp1+alph*DBjp1;
            
            f15 =   0.25*((2-4*j^2)*alph*Bjp1             +4*alph^2*DBjp1                          +alph^3*D2Bjp1);
            Df15 = 0.25*((2-4*j^2)*(Bjp1+alph*DBjp1)+4*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1));
            
            f16 =   0.5*(-alph)*Bjp1                    +3*(-alph^2)*Cj                        +1.5*(-alph^2)*Cjp2;
            Df16 = 0.5*((-1)*Bjp1+(-alph)*DBjp1)+3*((-2*alph)*Cj+(-alph^2)*DCj)+1.5*((-2*alph)*Cjp2+(-alph^2)*DCjp2);
            
            f17 =   1/64*((12+64*j+109*j^2+72*j^3+16*j^4)*Ajp2  +(-12-28*j-16*j^2)*alph*DAjp2               +(6-14*j-8*j^2)*alph^2*D2Ajp2                          +8*alph^3*D3Ajp2                              +alph^4*D4Ajp2);
            Df17 = 1/64*((12+64*j+109*j^2+72*j^3+16*j^4)*DAjp2+(-12-28*j-16*j^2)*(DAjp2+alph*D2Ajp2)+(6-14*j-8*j^2)*(2*alph*D2Ajp2+alph^2*D3Ajp2)+8*(3*alph^2*D3Ajp2+alph^3*D4Ajp2)+(4*alph^3*D4Ajp2+alph^4*D5Ajp2));
            
            f18 =   1/16*((12-15*j+4*j^2)*alph*Bjm1              +(8-4*j)*alph^2*DBjm1                           +alph^3*D2Bjm1);
            Df18 = 1/16*((12-15*j+4*j^2)*(Bjm1+alph*DBjm1)+(8-4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f19 =   1/8*((6*j-4*j^2)*alph*Bj         +(-4+4*j)*alph^2*DBj                      -alph^3*D2Bj);
            Df19 = 1/8*((6*j-4*j^2)*(Bj+alph*DBj)+(-4+4*j)*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f20 =   1/16*((3*j+4*j^2)*alph*Bjp1             -4*j*alph^2*DBjp1                         +alph^3*D2Bjp1);
            Df20 = 1/16*((3*j+4*j^2)*(Bjp1+alph*DBjp1)-4*j*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1));
            
            f21 =   1/8*((-12+15*j-4*j^2)*alph*Bjm1              +(-8+4*j)*alph^2*DBjm1                          -alph^3*D2Bjm1);
            Df21 = 1/8*((-12+15*j-4*j^2)*(Bjm1+alph*DBjm1)+(-8+4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f22 =   0.25*((6*j+4*j^2)*alph*Bj         -4*alph^2*DBj                      -alph^3*D2Bj);
            Df22 = 0.25*((6*j+4*j^2)*(Bj+alph*DBj)-4*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f23 =   0.25*((6*j+4*j^2)*alph*Bjp2             -4*alph^2*DBjp2                         -alph^3*D2Bjp2);
            Df23 = 0.25*((6*j+4*j^2)*(Bjp2+alph*DBjp2)-4*(2*alph*DBjp2+alph^2*D2Bjp2)-(3*alph^2*D2Bjp2+alph^3*D3Bjp2));
            
            f24 =   0.25*((-6*j+4*j^2)*alph*Bj         +(4-4*j)*alph^2*DBj                      +alph^3*D2Bj);
            Df24 = 0.25*((-6*j+4*j^2)*(Bj+alph*DBj)+(4-4*j)*(2*alph*DBj+alph^2*D2Bj)+(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f25 =   1/8*((-3*j-4*j^2)*alph*Bjp1             +4*j*alph^2*DBjp1                         -alph^3*D2Bjp1);
            Df25 = 1/8*((-3*j-4*j^2)*(Bjp1+alph*DBjp1)+4*j*(2*alph*DBjp1+alph^2*D2Bjp1)-(3*alph^2*D2Bjp1+alph^3*D3Bjp1));
            
            f26 =   0.5*alph*Bjp1             +0.75*alph^2*Cj                    +1.5*alph^2*Cjp2;
            Df26 = 0.5*(Bjp1+alph*DBjp1)+0.75*(2*alph*Cj+alph^2*DCj)+1.5*(2*alph*Cjp2+alph^2*DCjp2);
            
            
            %1st order semi-major axis functions
            f27 =   0.5*(-2*j*Aj-alph*DAj);
            Df27 = 0.5*(-2*j*DAj-DAj-alph*D2Aj);
            
            f28 =   1/16*((2*j-10*j^2+8*j^3)*Aj  +(3-7*j+4*j^2)*alph*DAj            +(-2-2*j)*alph^2*D2Aj                      -alph^3*D3Aj);
            Df28 = 1/16*((2*j-10*j^2+8*j^3)*DAj+(3-7*j+4*j^2)*(DAj+alph*D2Aj)+(-2-2*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj));
            
            f29 =   1/8*((8*j^3)*Aj  +(-2-4*j+4*j^2)*alph*DAj            +(-4-2*j)*alph^2*D2Aj                      -alph^3*D3Aj);
            Df29 = 1/8*((8*j^3)*DAj+(-2-4*j+4*j^2)*(DAj+alph*D2Aj)+(-4-2*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj));
            
            f30 =   0.25*((1+2*j)*alph*(Bjm1+Bjp1)                          +alph^2*(DBjm1+DBjp1));
            Df30 = 0.25*((1+2*j)*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))+(2*alph*(DBjm1+DBjp1)+alph^2*(D2Bjm1+D2Bjp1)));
            
            f31 =   0.5*((-1+2*j)*Ajm1+alph*DAjm1);
            Df31 = 0.5*((-1+2*j)*DAjm1+DAjm1+alph*D2Ajm1);
            
            f32 =   1/8*((4-16*j+20*j^2-8*j^3)*Ajm1  +(-4+12*j-4*j^2)*alph*DAjm1                +(3+2*j)*alph^2*D2Ajm1                           +alph^3*D3Ajm1);
            Df32 = 1/8*((4-16*j+20*j^2-8*j^3)*DAjm1+(-4+12*j-4*j^2)*(DAjm1+alph*D2Ajm1)+(3+2*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1));
            
            f33 =   1/16*((-2-j+10*j^2-8*j^3)*Ajm1   +(2+9*j-4*j^2)*alph*DAjm1               +(5+2*j)*alph^2*D2Ajm1                           +alph^3*D3Ajm1);
            Df33 = 1/16*((-2-j+10*j^2-8*j^3)*DAjm1+(2+9*j-4*j^2)*(DAjm1+alph*D2Ajm1)+(5+2*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1));
            
            f34 =   0.25*(-2*j*alph*(Bjm2+Bj)                       -alph^2*(DBjm2+DBj));
            Df34 = 0.25*(-2*j*((Bjm2+Bj)+alph*(DBjm2+DBj))-(2*alph*(DBjm2+DBj)+alph^2*(D2Bjm2+D2Bj)));
            
            f35 =   1/16*((1-j-10*j^2-8*j^3)*Ajp1  +(-1-j-4*j^2)*alph*DAjp1               +(3+2*j)*alph^2*D2Ajp1                          +alph^3*D3Ajp1);
            Df35 = 1/16*((1-j-10*j^2-8*j^3)*DAjp1+(-1-j-4*j^2)*(DAjp1+alph*D2Ajp1)+(3+2*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)+(3*alph^2*D3Ajp1+alph^3*D4Ajp1));
            
            f36 =   1/16*((-8+32*j-30*j^2+8*j^3)*Ajm2  +(8-17*j+4*j^2)*alph*DAjm2                +(-4-2*j)*alph^2*D2Ajm2                           -alph^3*D3Ajm2);
            Df36 = 1/16*((-8+32*j-30*j^2+8*j^3)*DAjm2+(8-17*j+4*j^2)*(DAjm2+alph*D2Ajm2)+(-4-2*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)-(3*alph^2*D3Ajm2+alph^3*D4Ajm2));
            
            f37 =   0.25*((-5+2*j)*alph*Bjm1              -alph^2*DBjm1);
            Df37 = 0.25*((-5+2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1));
            
            f38 =   0.25*(-2*j*alph*Bj         +alph^2*DBj);
            Df38 = 0.25*(-2*j*(Bj+alph*DBj)+(2*alph*DBj+alph^2*D2Bj));
            
            f39 =   0.5*((-1-2*j)*alph*Bjm1              -alph^2*DBjm1);
            Df39 = 0.5*((-1-2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1));
            
            f40 =   0.5*((-1-2*j)*alph*Bjp1             -alph^2*DBjp1);
            Df40 = 0.5*((-1-2*j)*(Bjp1+alph*DBjp1)-(2*alph*DBjp1+alph^2*D2Bjp1));
            
            f41 =   0.5*((5-2*j)*alph*Bjm1              +alph^2*DBjm1);
            Df41 = 0.5*((5-2*j)*(Bjm1+alph*DBjm1)+(2*alph*DBjm1+alph^2*D2Bjm1));
            
            f42 =   0.5*(2*j*alph*Bjm2              +alph^2*DBjm2);
            Df42 = 0.5*(2*j*(Bjm2+alph*DBjm2)+(2*alph*DBjm2+alph^2*D2Bjm2));
            
            f43 =   0.5*(2*j*alph*Bj         +alph^2*DBj);
            Df43 = 0.5*(2*j*(Bj+alph*DBj)+(2*alph*DBj+alph^2*D2Bj));
            
            f44 =   0.5*(2*j*alph*Bj         -alph^2*DBj);
            Df44 = 0.5*(2*j*(Bj+alph*DBj)-(2*alph*DBj+alph^2*D2Bj));
            
            %2nd order semi-major axis functions
            f45 =   1/8*((-5*j+4*j^2)*Aj  +(-2+4*j)*alph*DAj            +alph^2*D2Aj);
            Df45 = 1/8*((-5*j+4*j^2)*DAj+(-2+4*j)*(DAj+alph*D2Aj)+(2*alph*D2Aj+alph^2*D3Aj));
            
            f46 =   1/96*((22*j-64*j^2+60*j^3-16*j^4)*Aj  +(16-46*j+48*j^2-16*j^3)*alph*DAj            +(-12+9*j)*alph^2*D2Aj                      +4*j*alph^3*D3Aj                          +alph^4*D4Aj);
            Df46 = 1/96*((22*j-64*j^2+60*j^3-16*j^4)*DAj+(16-46*j+48*j^2-16*j^3)*(DAj+alph*D2Aj)+(-12+9*j)*(2*alph*D2Aj+alph^2*D3Aj)+4*j*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj));
            
            f47 =   1/32*((20*j^3-16*j^4)*Aj  +(-4-2*j+16*j^2-16*j^3)*alph*DAj            +(-2+11*j)*alph^2*D2Aj                      +(4+4*j)*alph^3*D3Aj                          +alph^4*D4Aj);
            Df47 = 1/32*((20*j^3-16*j^4)*DAj+(-4-2*j+16*j^2-16*j^3)*(DAj+alph*D2Aj)+(-2+11*j)*(2*alph*D2Aj+alph^2*D3Aj)+(4+4*j)*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj));
            
            f48 =   1/16*((2+j-4*j^2)*alph*(Bjm1+Bjp1)                          -4*j*alph^2*(DBjm1+DBjp1)                                         -alph^3*(D2Bjm1+D2Bjp1));
            Df48 = 1/16*((2+j-4*j^2)*((Bjm1+Bjp1)+alph*(DBjm1+DBjp1))-4*j*(2*alph*(DBjm1+DBjp1)+alph^2*(D2Bjm1+D2Bjp1))-(3*alph^2*(D2Bjm1+D2Bjp1)+alph^3*(D3Bjm1+D3Bjp1)));
            
            f49 =   0.25*((-2+6*j-4*j^2)*Ajm1  +(2-4*j)*alph*DAjm1                -alph^2*D2Ajm1);
            Df49 = 0.25*((-2+6*j-4*j^2)*DAjm1+(2-4*j)*(DAjm1+alph*D2Ajm1)-(2*alph*D2Ajm1+alph^2*D3Ajm1));
            
            f50 =   1/32*((20-86*j+126*j^2-76*j^3+16*j^4)*Ajm1  +(-20+74*j-64*j^2+16*j^3)*alph*DAjm1                +(14-17*j)*alph^2*D2Ajm1                           +(-2-4*j)*alph^3*D3Ajm1                              -alph^4*D4Ajm1); %in this term there appears Aj- instead of Aj-1; this is explained in the erratum http://ssdbook.maths.qmul.ac.uk/errors.pdf
            Df50 = 1/32*((20-86*j+126*j^2-76*j^3+16*j^4)*DAjm1+(-20+74*j-64*j^2+16*j^3)*(DAjm1+alph*D2Ajm1)+(14-17*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(-2-4*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1));
            
            f51 =   1/32*((-4+2*j+22*j^2-36*j^3+16*j^4)*Ajm1  +(4+6*j-32*j^2+16*j^3)*alph*DAjm1                +(-2-19*j)*alph^2*D2Ajm1                           +(-6-4*j)*alph^3*D3Ajm1                               -alph^4*D4Ajm1);
            Df51 = 1/32*((-4+2*j+22*j^2-36*j^3+16*j^4)*DAjm1+(4+6*j-32*j^2+16*j^3)*(DAjm1+alph*D2Ajm1)+(-2-19*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(-6-4*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1));
            
            f52 =   1/8*((-2*j+4*j^2)*alph*(Bjm2+Bj)                       +4*j*alph^2*(DBjm2+DBj)                                     +alph^3*(D2Bjm2+D2Bj));
            Df52 = 1/8*((-2*j+4*j^2)*((Bjm2+Bj)+alph*(DBjm2+DBj))+4*j*(2*alph*(DBjm2+DBj)+alph^2*(D2Bjm2+D2Bj))+(3*alph^2*(D2Bjm2+D2Bj)+alph^3*(D3Bjm2+D3Bj)));
            
            f53 =   1/8*((2-7*j+4*j^2)*Ajm2   +(-2+4*j)*alph*DAjm2               +alph^2*D2Ajm2);
            Df53 = 1/8*((2-7*j+4*j^2)*DAjm2+(-2+4*j)*(DAjm2+alph*D2Ajm2)+(2*alph*D2Ajm2+alph^2*D3Ajm2));
            
            f54 =   1/32*((-32+144*j-184*j^2+92*j^3-16*j^4)*Ajm2  +(32-102*j+80*j^2-16*j^3)*alph*DAjm2                +(-16+25*j)*alph^2*D2Ajm2                           +(4+4*j)*alph^3*D3Ajm2                              +alph^4*D4Ajm2);
            Df54 = 1/32*((-32+144*j-184*j^2+92*j^3-16*j^4)*DAjm2+(32-102*j+80*j^2-16*j^3)*(DAjm2+alph*D2Ajm2)+(-16+25*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(4+4*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2));
            
            f55 =   1/96*((12-14*j-40*j^2+52*j^3-16*j^4)*Ajm2  +(-12-10*j+48*j^2-16*j^3)*alph*DAjm2                +(6+27*j)*alph^2*D2Ajm2                           +(8+4*j)*alph^3*D3Ajm2                               +alph^4*D4Ajm2);
            Df55 = 1/96*((12-14*j-40*j^2+52*j^3-16*j^4)*DAjm2+(-12-10*j+48*j^2-16*j^3)*(DAjm2+alph*D2Ajm2)+(6+27*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(8+4*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2));
            
            f56 =   1/16*((3*j-4*j^2)*alph*(Bjm3+Bjm1)                            -4*j*alph^2*(DBjm3+DBjm1)                                         -alph^3*(D2Bjm3+D2Bjm1));
            Df56 = 1/16*((3*j-4*j^2)*((Bjm3+Bjm1)+alph*(DBjm3+DBjm1))-4*j*(2*alph*(DBjm3+DBjm1)+alph^2*(D2Bjm3+D2Bjm1))-(3*alph^2*(D2Bjm3+D2Bjm1)+alph^3*(D3Bjm3+D3Bjm1)));
            
            f57 =   0.5*alph*Bjm1;
            Df57 = 0.5*(Bjm1+alph*DBjm1);
            
            f58 =   1/8*((-14+16*j-4*j^2)*alph*Bjm1              +4*alph^2*DBjm1                           +alph^3*D2Bjm1);
            Df58 = 1/8*((-14+16*j-4*j^2)*(Bjm1+alph*DBjm1)+4*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f59 =   1/8*((2-4*j^2)*alph*Bjm1              +4*alph^2*DBjm1                          +alph^3*D2Bjm1);
            Df59 = 1/8*((2-4*j^2)*(Bjm1+alph*DBjm1)+4*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f60 =   0.75*(-alph^2)*(Cjm2+Cj);
            Df60 = 0.75*((-2*alph)*(Cjm2+Cj)+(-alph^2)*(DCjm2+DCj));
            
            f61 =   0.5*(-alph)*Bjm1                     +0.75*(-alph^2)*Cjm2                             +15/4*(-alph^2)*Cj;
            Df61 = 0.5*((-1)*Bjm1+(-alph)*DBjm1)+0.75*((-2*alph)*Cjm2+(-alph^2)*DCjm2)+15/4*((-2*alph)*Cj+(-alph^2)*DCj);
            
            f62 =   -alph*Bjm1;
            Df62 = -Bjm1-alph*DBjm1;
            
            f63 =   0.25*((14-16*j+4*j^2)*alph*Bjm1              -4*alph^2*DBjm1                          -alph^3*D2Bjm1);
            Df63 = 0.25*((14-16*j+4*j^2)*(Bjm1+alph*DBjm1)-4*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f64 =   0.25*((-2+4*j^2)*alph*Bjm1              -4*alph^2*DBjm1                          -alph^3*D2Bjm1);
            Df64 = 0.25*((-2+4*j^2)*(Bjm1+alph*DBjm1)-4*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f65 =   0.5*alph*Bjm1              +3*alph^2*Cjm2                         +1.5*alph^2*Cj;
            Df65 = 0.5*(Bjm1+alph*DBjm1)+3*(2*alph*Cjm2+alph^2*DCjm2)+1.5*(2*alph*Cj+alph^2*DCj);
            
            f66 =   0.5*alph*Bjm1              +1.5*alph^2*Cjm2                         +3*alph^2*Cj;
            Df66 = 0.5*(Bjm1+alph*DBjm1)+1.5*(2*alph*Cjm2+alph^2*DCjm2)+3*(2*alph*Cj+alph^2*DCj);
            
            f67 =   0.5*(-alph)*Bjm1                     +15/4*(-alph^2)*Cjm2                             +0.75*(-alph^2)*Cj;
            Df67 = 0.5*((-1)*Bjm1+(-alph)*DBjm1)+15/4*((-2*alph)*Cjm2+(-alph^2)*DCjm2)+0.75*((-2*alph)*Cj+(-alph^2)*DCj);
            
            f68 =   1/96*((4-2*j-26*j^2-4*j^3+16*j^4)*Ajp1  +(-4-2*j+16*j^3)*alph*DAjp1               +(6-3*j)*alph^2*D2Ajp1                          +(-2-4*j)*alph^3*D3Ajp1                             -alph^4*D4Ajp1);
            Df68 = 1/96*((4-2*j-26*j^2-4*j^3+16*j^4)*DAjp1+(-4-2*j+16*j^3)*(DAjp1+alph*D2Ajp1)+(6-3*j)*(2*alph*D2Ajp1+alph^2*D3Ajp1)+(-2-4*j)*(3*alph^2*D3Ajp1+alph^3*D4Ajp1)-(4*alph^3*D4Ajp1+alph^4*D5Ajp1));
            
            f69 =   1/96*((36-186*j+238*j^2-108*j^3+16*j^4)*Ajm3  +(-36+130*j-96*j^2+16*j^3)*alph*DAjm3                +(18-33*j)*alph^2*D2Ajm3                           +(-6-4*j)*alph^3*D3Ajm3                               -alph^4*D4Ajm3);
            Df69 = 1/96*((36-186*j+238*j^2-108*j^3+16*j^4)*DAjm3+(-36+130*j-96*j^2+16*j^3)*(DAjm3+alph*D2Ajm3)+(18-33*j)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(-6-4*j)*(3*alph^2*D3Ajm3+alph^3*D4Ajm3)-(4*alph^3*D4Ajm3+alph^4*D5Ajm3));
            
            f70 =   1/8*((-14*j+4*j^2)*alph*Bjm2              -8*alph^2*DBjm2                          -alph^3*D2Bjm2);
            Df70 = 1/8*((-14*j+4*j^2)*(Bjm2+alph*DBjm2)-8*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2));
            
            f71 =   1/8*((-2*j+4*j^2)*alph*Bj         -alph^3*D2Bj);
            Df71 = 1/8*((-2*j+4*j^2)*(Bj+alph*DBj)-(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f72 =   1/8*((-2-j+4*j^2)*alph*Bjm1              +4*j*alph^2*DBjm1                          +alph^3*D2Bjm1);
            Df72 = 1/8*((-2-j+4*j^2)*(Bjm1+alph*DBjm1)+4*j*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f73 =   1/8*((-2-j+4*j^2)*alph*Bjp1             +4*j*alph^2*DBjp1                         +alph^3*D2Bjp1);
            Df73 = 1/8*((-2-j+4*j^2)*(Bjp1+alph*DBjp1)+4*j*(2*alph*DBjp1+alph^2*D2Bjp1)+(3*alph^2*D2Bjp1+alph^3*D3Bjp1));
            
            f74 =   0.25*((2*j-4*j^2)*alph*Bjm2              -4*j*alph^2*DBjm2                          -alph^3*D2Bjm2);
            Df74 = 0.25*((2*j-4*j^2)*(Bjm2+alph*DBjm2)-4*j*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2));
            
            f75 =   0.25*((2*j-4*j^2)*alph*Bj         -4*j*alph^2*DBj                      -alph^3*D2Bj);
            Df75 = 0.25*((2*j-4*j^2)*(Bj+alph*DBj)-4*j*(2*alph*DBj+alph^2*D2Bj)-(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f76 =   0.25*((14*j-4*j^2)*alph*Bjm2              +8*alph^2*DBjm2                          +alph^3*D2Bjm2);
            Df76 = 0.25*((14*j-4*j^2)*(Bjm2+alph*DBjm2)+8*(2*alph*DBjm2+alph^2*D2Bjm2)+(3*alph^2*D2Bjm2+alph^3*D3Bjm2));
            
            f77 =   0.25*((2*j-4*j^2)*alph*Bj         +alph^3*D2Bj);
            Df77 = 0.25*((2*j-4*j^2)*(Bj+alph*DBj)+(3*alph^2*D2Bj+alph^3*D3Bj));
            
            f78 =   1/8*((-3*j+4*j^2)*alph*Bjm3              +4*j*alph^2*DBjm3                           +alph^3*D2Bjm3);
            Df78 = 1/8*((-3*j+4*j^2)*(Bjm3+alph*DBjm3)+4*j*(2*alph*DBjm3+alph^2*D2Bjm3)+(3*alph^2*D2Bjm3+alph^3*D3Bjm3));
            
            f79 =   1/8*((-3*j+4*j^2)*alph*Bjm1              +4*j*alph^2*DBjm1                          +alph^3*D2Bjm1);
            Df79 = 1/8*((-3*j+4*j^2)*(Bjm1+alph*DBjm1)+4*j*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f80 =   1.5*alph^2*Cj;
            Df80 = 1.5*(2*alph*Cj+alph^2*DCj);
            
            f81 =   1.5*alph^2*Cjm2;
            Df81 = 1.5*(2*alph*Cjm2+alph^2*DCjm2);
            
            %3rd order semi-major axis functions
            f82 =   1/48*((-26*j+30*j^2-8*j^3)*Aj  +(-9+27*j-12*j^2)*alph*DAj            +(6-6*j)*alph^2*D2Aj                     -alph^3*D3Aj);
            Df82 = 1/48*((-26*j+30*j^2-8*j^3)*DAj+(-9+27*j-12*j^2)*(DAj+alph*D2Aj)+(6-6*j)*(2*alph*D2Aj+alph^2*D3Aj)-(3*alph^2*D3Aj+alph^3*D4Aj));
            
            f83 =   1/16*((-9+31*j-30*j^2+8*j^3)*Ajm1  +(9-25*j+12*j^2)*alph*DAjm1                +(-5+6*j)*alph^2*D2Ajm1                           +alph^3*D3Ajm1);
            Df83 = 1/16*((-9+31*j-30*j^2+8*j^3)*DAjm1+(9-25*j+12*j^2)*(DAjm1+alph*D2Ajm1)+(-5+6*j)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(3*alph^2*D3Ajm1+alph^3*D4Ajm1));
            
            f84 =   1/16*((8-32*j+30*j^2-8*j^3)*Ajm2  +(-8+23*j-12*j^2)*alph*DAjm2                +(4-6*j)*alph^2*D2Ajm2                           -alph^3*D3Ajm2);
            Df84 = 1/16*((8-32*j+30*j^2-8*j^3)*DAjm2+(-8+23*j-12*j^2)*(DAjm2+alph*D2Ajm2)+(4-6*j)*(2*alph*D2Ajm2+alph^2*D3Ajm2)-(3*alph^2*D3Ajm2+alph^3*D4Ajm2));
            
            f85 =   1/48*((-6+29*j-30*j^2+8*j^3)*Ajm3  +(6-21*j+12*j^2)*alph*DAjm3                +(-3+6*j)*alph^2*D2Ajm3                           +alph^3*D3Ajm3);
            Df85 = 1/48*((-6+29*j-30*j^2+8*j^3)*DAjm3+(6-21*j+12*j^2)*(DAjm3+alph*D2Ajm3)+(-3+6*j)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(3*alph^2*D3Ajm3+alph^3*D4Ajm3));
            
            f86 =   0.25*((3-2*j)*alph*Bjm1              -alph^2*DBjm1);
            Df86 = 0.25*((3-2*j)*(Bjm1+alph*DBjm1)-(2*alph*DBjm1+alph^2*D2Bjm1));
            
            f87 =   0.25*(2*j*alph*Bjm2              +alph^2*DBjm2);
            Df87 = 0.25*(2*j*(Bjm2+alph*DBjm2)+(2*alph*DBjm2+alph^2*D2Bjm2));
            
            f88 =   0.5*((-3+2*j)*alph*Bjm1              +alph^2*DBjm1);
            Df88 = 0.5*((-3+2*j)*(Bjm1+alph*DBjm1)+(2*alph*DBjm1+alph^2*D2Bjm1));
            
            f89 =   0.5*(-2*j*alph*Bjm2              -alph^2*DBjm2);
            Df89 = 0.5*(-2*j*(Bjm2+alph*DBjm2)-(2*alph*DBjm2+alph^2*D2Bjm2));
            
            %4th order semi-major axis functions
            f90 =   1/384*((-206*j+283*j^2-120*j^3+16*j^4)*Aj   +(-64+236*j-168*j^2+32*j^3)*alph*DAj           +(48-78*j+24*j^2)*alph^2*D2Aj                     +(-12+8*j)*alph^3*D3Aj                           +alph^4*D4Aj);
            Df90 = 1/384*((-206*j+283*j^2-120*j^3+16*j^4)*DAj+(-64+236*j-168*j^2+32*j^3)*(DAj+alph*D2Aj)+(48-78*j+24*j^2)*(2*alph*D2Aj+alph^2*D3Aj)+(-12+8*j)*(3*alph^2*D3Aj+alph^3*D4Aj)+(4*alph^3*D4Aj+alph^4*D5Aj));
            
            f91 =   1/96*((-64+238*j-274*j^2+116*j^3-16*j^4)*Ajm1 +(64-206*j+156*j^2-32*j^3)*alph*DAjm1                 +(-36+69*j-24*j^2)*alph^2*D2Ajm1                           +(10-8*j)*alph^3*D3Ajm1                              -alph^4*D4Ajm1);
            Df91 = 1/96*((-64+238*j-274*j^2+116*j^3-16*j^4)*DAjm1+(64-206*j+156*j^2-32*j^3)*(DAjm1+alph*D2Ajm1)+(-36+69*j-24*j^2)*(2*alph*D2Ajm1+alph^2*D3Ajm1)+(10-8*j)*(3*alph^2*D3Ajm1+alph^3*D4Ajm1)-(4*alph^3*D4Ajm1+alph^4*D5Ajm1));
            
            f92 =   1/64*((52-224*j+259*j^2-112*j^3+16*j^4)*Ajm2  +(-52+176*j-144*j^2+32*j^3)*alph*DAjm2                +(26-60*j+24*j^2)*alph^2*D2Ajm2                           +(-8+8*j)*alph^3*D3Ajm2                               +alph^4*D4Ajm2);
            Df92 = 1/64*((52-224*j+259*j^2-112*j^3+16*j^4)*DAjm2+(-52+176*j-144*j^2+32*j^3)*(DAjm2+alph*D2Ajm2)+(26-60*j+24*j^2)*(2*alph*D2Ajm2+alph^2*D3Ajm2)+(-8+8*j)*(3*alph^2*D3Ajm2+alph^3*D4Ajm2)+(4*alph^3*D4Ajm2+alph^4*D5Ajm2));
            
            f93 =   1/96*((-36+186*j-238*j^2+108*j^3-16*j^4)*Ajm3   +(36-146*j+132*j^2-32*j^3)*alph*DAjm3               +(-18+51*j-24*j^2)*alph^2*D2Ajm3                           +(6-8*j)*alph^3*D3Ajm3                              -alph^4*D4Ajm3);
            Df93 = 1/96*((-36+186*j-238*j^2+108*j^3-16*j^4)*DAjm3+(36-146*j+132*j^2-32*j^3)*(DAjm3+alph*D2Ajm3)+(-18+51*j-24*j^2)*(2*alph*D2Ajm3+alph^2*D3Ajm3)+(6-8*j)*(3*alph^2*D3Ajm3+alph^3*D4Ajm3)-(4*alph^3*D4Ajm3+alph^4*D5Ajm3));
            
            f94 =   1/384*((24-146*j+211*j^2-104*j^3+16*j^4)*Ajm4  +(-24+116*j-120*j^2+32*j^3)*alph*DAjm4                +(12-42*j+24*j^2)*alph^2*D2Ajm4                           +(-4+8*j)*alph^3*D3Ajm4                               +alph^4*D4Ajm4);
            Df94 = 1/384*((24-146*j+211*j^2-104*j^3+16*j^4)*DAjm4+(-24+116*j-120*j^2+32*j^3)*(DAjm4+alph*D2Ajm4)+(12-42*j+24*j^2)*(2*alph*D2Ajm4+alph^2*D3Ajm4)+(-4+8*j)*(3*alph^2*D3Ajm4+alph^3*D4Ajm4)+(4*alph^3*D4Ajm4+alph^4*D5Ajm4));
            
            f95 =   1/16*((16-17*j+4*j^2)*alph*Bjm1              +(-8+4*j)*alph^2*DBjm1                           +alph^3*D2Bjm1);
            Df95 = 1/16*((16-17*j+4*j^2)*(Bjm1+alph*DBjm1)+(-8+4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)+(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f96 =   1/8*((10*j-4*j^2)*alph*Bjm2              +(4-4*j)*alph^2*DBjm2                           -alph^3*D2Bjm2);
            Df96 = 1/8*((10*j-4*j^2)*(Bjm2+alph*DBjm2)+(4-4*j)*(2*alph*DBjm2+alph^2*D2Bjm2)-(3*alph^2*D2Bjm2+alph^3*D3Bjm2));
            
            f97 =   1/16*((-3*j+4*j^2)*alph*Bjm3              +4*j*alph^2*DBjm3                          +alph^3*D2Bjm3);
            Df97 = 1/16*((-3*j+4*j^2)*(Bjm3+alph*DBjm3)+4*j*(2*alph*DBjm3+alph^2*D2Bjm3)+(3*alph^2*D2Bjm3+alph^3*D3Bjm3));
            
            f98 =   3/8*alph^2*Cjm2;
            Df98 = 3/8*(2*alph*Cjm2+alph^2*DCjm2);
            
            f99 =   1/8*((-16+17*j-4*j^2)*alph*Bjm1              +(8-4*j)*alph^2*DBjm1                          -alph^3*D2Bjm1);
            Df99 = 1/8*((-16+17*j-4*j^2)*(Bjm1+alph*DBjm1)+(8-4*j)*(2*alph*DBjm1+alph^2*D2Bjm1)-(3*alph^2*D2Bjm1+alph^3*D3Bjm1));
            
            f100 =   0.25*((-10*j+4*j^2)*alph*Bjm2              +(-4+4*j)*alph^2*DBjm2                          +alph^3*D2Bjm2);
            Df100 = 0.25*((-10*j+4*j^2)*(Bjm2+alph*DBjm2)+(-4+4*j)*(2*alph*DBjm2+alph^2*D2Bjm2)+(3*alph^2*D2Bjm2+alph^3*D3Bjm2));
            
            f101 =   1/8*((3*j-4*j^2)*alph*Bjm3              -4*j*alph^2*DBjm3                          -alph^3*D2Bjm3);
            Df101 = 1/8*((3*j-4*j^2)*(Bjm3+alph*DBjm3)-4*j*(2*alph*DBjm3+alph^2*D2Bjm3)-(3*alph^2*D2Bjm3+alph^3*D3Bjm3));
            
            f102 =   1.5*(-alph^2)*Cjm2;
            Df102 = 1.5*((-2*alph)*Cjm2+(-alph^2)*DCjm2);
            
            f103 =   9/4*(alph^2*Cjm2);
            Df103 = 9/4*(2*alph*Cjm2+alph^2*DCjm2);
            
            
            end
            
            
            %% PrepareExpansionTerms;
            
            %Write down the disturbing function terms powers and resonant orders - Solar
            %System Dynamics (Murray & Dermott 1999), Appendix B
            
            for DummyIdx = 1:1
            
            if ~isfield(Terms,'k') %constant fields, not depending on the value of alpha or j, should be written once and for all
                %  if isempty(Terms) %constant fields, not depending on the value of alpha or j, should be written once and for all
                
                Terms.k = [zeros(1,38) ... %0th order
                    ones(1,22) ... %1st order
                    2*ones(1,46) ... %2nd order
                    3*ones(1,10) ... %3rd order
                    4*ones(1,19)];       %4th order
                
                Terms.A = [0 2 0 0 0 4 2 0 2 0 2 0 0 0 0 ...
                    1 3 1 1 1 0 2 0 0 0 2 2 1 0 2 1 1 1 0 2 1 0 0 ... %0th order
                    1 3 1 1 1 0 2 0 0 0 2 1 1 0 1 1 1 0 0 0 1 0 ... %1st order
                    2 4 2 2 2 1 3 1 1 1 0 2 0 0 0 0 2 0 0 0 ...
                    0 2 0 0 0 0 2 0 0 0 3 1 1 1 2 2 1 1 1 1 0 0 0 1 1 0 ... %2nd order
                    3 2 1 0 1 0 1 0 1 0 ... %3rd order
                    4 3 2 1 0 2 1 0 0 2 1 0 0 2 1 0 0 0 0 ... %4th order
                    ];
                
                Terms.Ap = [0 0 2 0 0 0 2 4 0 2 0 2 0 0 0 ...
                    1 1 3 1 1 0 0 2 0 0 2 0 1 2 0 1 1 1 2 0 1 2 0 ...  %0th order
                    0 0 2 0 0 1 1 3 1 1 1 2 0 1 0 0 0 1 1 1 0 1 ... %1st order
                    0 0 2 0 0 1 1 3 1 1 2 2 4 2 2 0 0 2 0 0 ...
                    0 0 2 0 0 0 0 2 0 0 1 3 1 1 0 0 1 1 1 1 2 2 0 1 1 0 ... %2nd order
                    0 1 2 3 0 1 0 1 0 1 ... %3rd order
                    0 1 2 3 4 0 1 2 0 0 1 2 0 0 1 2 0 0 0 ... %4th order
                    ];
                
                Terms.B = [0 0 0 2 0 0 0 0 2 2 0 0 4 0 2 ...
                    0 0 0 2 0 1 1 1 3 1 0 2 2 2 1 1 1 1 1 0 0 0 2 ...  %0th order
                    0 0 0 2 0 0 0 0 2 0 0 0 2 2 1 1 1 1 1 1 0 0 ... %1st order
                    0 0 0 2 0 0 0 0 2 0 0 0 0 2 0 2 2 2 4 2 ...
                    1 1 1 3 1 0 0 0 2 0 0 0 2 2 1 1 1 1 1 1 1 1 3 0 0 1 ... %2nd order
                    0 0 0 0 2 2 1 1 0 0 ... %3rd order
                    0 0 0 0 0 2 2 2 4 1 1 1 3 0 0 0 2 1 0 ... %4th order
                    ];
                
                Terms.Bp = [0 0 0 0 2 0 0 0 0 0 2 2 0 4 2 ...
                    0 0 0 0 2 1 1 1 1 3 0 0 0 0 1 1 1 1 1 2 2 2 2 ...  %0th order
                    0 0 0 0 2 0 0 0 0 2 0 0 0 0 1 1 1 1 1 1 2 2 ... %1st order
                    0 0 0 0 2 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2 ...
                    1 1 1 1 3 2 2 2 2 4 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 3 ... %2nd order
                    0 0 0 0 0 0 1 1 2 2 ... %3rd order
                    0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 3 4 ... %4th order
                    ];
                
                Terms.C = [zeros(1,15) ...
                    ones(1,5) zeros(1,5) 2 -2 -1 0 -2 1 1 -1 0 -2 -1 0 0 ...  %0th order
                    1 1 1 1 1 0 0 0 0 0 2 -1 -1 0 1 1 -1 0 0 0 -1 0 ... %1st order
                    2 2 2 2 2 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 ...
                    0 0 0 0 0 0 0 0 0 0 3 -1 -1 1 2 2 1 1 -1 1 0 0 0 -1 1 0 ... %2nd order
                    3 2 1 0 1 0 1 0 1 0 ... %3rd order
                    4 3 2 1 0 2 1 0 0 2 1 0 0 2 1 0 0 0 0 ... %4th order
                    ];
                
                Terms.Cp = [zeros(1,15) ...
                    -ones(1,5) zeros(1,5) -2 0 -1 -2 0 -1 -1 -1 -2 0 -1 -2 0 ...  %0th order
                    0 0 0 0 0 1 1 1 1 1 -1 2 0 -1 0 0 0 1 1 -1 0 -1 ... %1st order
                    0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 0 0 0 0 0 ...
                    0 0 0 0 0 0 0 0 0 0 -1 3 1 -1 0 0 1 1 1 -1 2 2 0 1 -1 0 ... %2nd order
                    0 1 2 3 0 1 0 1 0 1 ... %3rd order
                    0 1 2 3 4 0 1 2 0 0 1 2 0 0 1 2 0 0 0 .... %4th order
                    ];
                
                Terms.D = [zeros(1,15) ...
                    zeros(1,5) ones(1,5) 0 2 2 2 1 -1 1 1 1 0 0 0 2 ...  %0th order
                    zeros(1,10) 0 0 2 2 -1 1 1 -1 1 1 0 0 ... %1st order
                    zeros(1,15) 2 2 2 2 2 1 1 1 1 1 0 0 0 0 0 ...
                    0 0 2 2 -1 1 -1 1 1 1 -1 1 3 0 0 -1 ... %2nd order
                    0 0 0 0 2 2 1 1 0 0 ... %3rd order
                    0 0 0 0 0 2 2 2 4 1 1 1 3 0 0 0 2 1 0 ... %4th order
                    ];
                
                Terms.Dp = [zeros(1,15) ...
                    zeros(1,5) -ones(1,5) 0 0 0 0 1 1 -1 1 1 2 2 2 -2 ...  %0th order
                    zeros(1,10) 0 0 0 0 1 -1 1 1 -1 1 2 2 ... %1st order
                    zeros(1,20) 1 1 1 1 1 2 2 2 2 2 ...
                    0 0 0 0 1 -1 1 -1 1 1 1 -1 -1 2 2 3 ... %2nd order
                    0 0 0 0 0 0 1 1 2 2 ... %3rd order
                    0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 3 4 ... %4th order
                    ];
                
                Terms.N = length(Terms.A);
            end
            
            Terms.f = [f1 f2 f2 f3 f3 f4 f5 f6 f7 f7 f7 f7 f8 f8 f9 ...
                f10 f11 f12 f13 f13 f14 f15 f15 f16 f16 f17 f18 f19 f20 f21 f22 f23 f24 f25 f18 f19 f20 f26 ...  %0th order
                f27 f28 f29 f30 f30 f31 f32 f33 f34 f34 f35 f36 f37 f38 f39 f40 f41 f42 f43 f44 f37 f38 ... %1st order
                f45 f46 f47 f48 f48 f49 f50 f51 f52 f52 f53 f54 f55 f56 f56 f57 f58 f59 f60 f61 f62 f63 f64 f65 f66 f57 f58 f59 f67 f60 ...
                f68 f69 f70 f71 f72 f73 f74 f75 f76 f77 f78 f79 f80 f70 f71 f81 ... %2nd order
                f82 f83 f84 f85 f86 f87 f88 f89 f86 f87 ... %3rd order
                f90 f91 f92 f93 f94 f95 f96 f97 f98 f99 f100 f101 f102 f95 f96 f97 f103 f102 f98 ... %4th order
                
                ];
            
            %indirect terms due to an external perturber. The prefactor of alph
            %comes from the definition of the disturbing function (Solar System
            %Dynamics equation 6.44).
            
            Terms.fE = alph*[(j==1)*[-1 0.5 0.5 1 1 1/64  -0.25 1/64 -0.5 -0.5 -0.5 -0.5 0 0 -1] ...
                (j==-2)*[-1  0.75 0.75 1 1] (j==-1)*[-2 1 1 1 1] (j==-1)*(-1/64)+(j==-3)*(-81/64) ...
                (j==1)*(-1/8) 0  (j==-1)*(-1/8) (j==1)*0.25 0 (j==-2)*(-2) 0 (j==-1)*0.25 ...
                (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==-1)*(-1) ... %up to here - 0th order indirect
                (j==-1)*[-0.5 3/8 0.25 0.5 0.5]+(j==1)*[1.5 0 -0.75 -1.5 -1.5] ...
                (j==2)*[-2 1 1.5 2 2] (j==-2)*(-0.75) (j==1)*3/16+(j==3)*(-27/16)...
                (j==1)*1.5 0 (j==1)*3 (j==-1)*(-1) (j==1)*(-3) (j==2)*(-4) 0 0  (j==1)*1.5 0 ... %up to here - 1st order indirect
                (j==-1)*[-3/8 3/8 3/16 3/8 3/8]+(j==1)*[-1/8 -1/24 1/16 1/8 1/8] (j==2)*[3 0 -9/4 -3 -3] ...
                (j==1)*[-1/8 1/16 -1/24 1/8 1/8]+(j==3)*[-27/8 27/16 27/8 27/8 27/8]...
                (j==1)*[-1 0.5 0.5 0 1] (j==1)*[2 -1 -1 -1 -1] (j==1)*[-1 0.5 0.5 1 0]...
                (j==-2)*(-2/3) (j==2)*(0.25)+(j==4)*(-8/3) (j==2)*3 0 ...
                (j==1)*(-0.25) (j==-1)*(-0.75) (j==2)*6 0 (j==2)*(-6) 0 ...
                (j==3)*(-27/4) (j==1)*(-0.25) 0 (j==2)*3 0 0 ... %up to here - 2nd order indirect
                (j==-1)*(-1/3)+(j==1)*(-1/24) (j==2)*(-0.25) (j==1)*(-1/16)+(j==3)*(81/16) ...
                (j==2)*(-1/6)+(j==4)*(-16/3) (j==1)*(-0.5) (j==2)*(-2) (j==1)*1 (j==2)*4 (j==1)*(-0.5) (j==2)*(-2) ... %up to here - 3rd order indirect
                (j==-1)*(-125/384)+(j==1)*(-3/128) (j==2)*(-1/12) (j==1)*(-3/64)+(j==3)*(-27/64) ...
                (j==2)*(-1/12)+(j==4)*(8) (j==3)*(-27/128)+(j==5)*(-3125/384) (j==1)*(-3/8) (j==2)*(-1) (j==3)*(-27/8) 0 ...
                (j==1)*(0.75) (j==2)*(2) (j==3)*(27/4) 0 (j==1)*(-3/8) (j==2)*(-1) (j==3)*(-27/8) 0 0 0 ... %up to here - 4th order indirect
                ];
            
            %Indirect terms due to an internal perturber. The prefactor of 1/alph^2
            %comes from the definition of the disturbing function (Solar System
            %Dynamics equation 6.45).
            
            
            Terms.fI = 1/alph^2*[(j==1)*[-1 0.5 0.5 1 1 1/64  -0.25 1/64 -0.5 -0.5 -0.5 -0.5 0 0 -1] ...
                +(j==-2)*[-1  0.75 0.75 1 1] (j==-1)*[-2 1 1 1 1] (j==-1)*(-1/64)+(j==-3)*(-81/64) ...
                (j==1)*(-1/8) 0  (j==-1)*(-1/8) (j==1)*0.25 0 (j==-2)*(-2) 0 (j==-1)*0.25 ...
                (j==1)*(-1/8) 0 (j==-1)*(-1/8) (j==-1)*(-1) ... %up to here - 0th order indirect
                (j==-1)*[-2 1.5 1 2 2] ...
                (j==0)*[1.5 -0.75 0 -1.5 -1.5]+(j==2)*[-0.5 0.25 3/8 0.5 0.5] ...
                (j==0)*3/16+(j==-2)*(-27/16)...
                (j==3)*(-0.75) 0 (j==0)*1.5 0 (j==-1)*(-4) 0 (j==2)*(-1) (j==0)*3 (j==0)*(-3) 0 (j==0)*1.5 ... %up to here - 1st order indirect
                (j==-1)*[-27/8 27/8 27/16 27/8 27/8]+(j==1)*[-1/8 -1/24 1/16 1/8 1/8]...
                (j==2)*[3 -9/4 0 -3 -3] (j==1)*[-1/8 1/16 -1/24 1/8 1/8]+(j==3)*[-3/8 3/16 3/8 3/8 3/8] ...
                (j==1)*[-1 0.5 0.5 0 1] (j==1)*[2 -1 -1 -1 -1] (j==1)*[-1 0.5 0.5 1 0] ...
                (j==0)*(0.25)+(j==-2)*(-8/3) (j==4)*(-2/3) 0 (j==0)*3 (j==1)*(-0.25) (j==-1)*(-27/4)...
                0 (j==0)*6 0 (j==0)*(-6) (j==3)*(-0.75) (j==1)*(-0.25) 0 0 (j==0)*3 0 ... %up to here - 2nd order indirect
                (j==-1)*(-16/3)+(j==1)*(-1/6) (j==0)*(81/16)+(j==2)*(-1/16) (j==1)*(-0.25) ...
                (j==2)*(-1/24)+(j==4)*(-1/3) (j==1)*(-2) (j==2)*(-0.5) (j==1)*4 (j==2)*1 (j==1)*(-2) (j==2)*(-0.5) ... %up to here - 3rd order indirect
                (j==-1)*(-3125/384)+(j==1)*(-27/128) (j==0)*(8)+(j==2)*(-1/12) (j==1)*(-27/64)+(j==3)*(-3/64) ...
                (j==2)*(-1/12) (j==3)*(-3/128)+(j==5)*(-125/384) (j==1)*(-27/8) (j==2)*(-1) (j==3)*(-3/8) 0 ...
                (j==1)*(27/4) (j==2)*(2) (j==3)*(0.75) 0 (j==1)*(-27/8) (j==2)*(-1) (j==3)*(-3/8) 0 0 0 ... %up to here - 4th order indirect
                ];
            
            
            Terms.Df = [Df1 Df2 Df2 Df3 Df3 Df4 Df5 Df6 Df7 Df7 Df7 Df7 Df8 Df8 Df9 ...
                Df10 Df11 Df12 Df13 Df13 Df14 Df15 Df15 Df16 Df16 Df17 Df18 Df19 Df20 Df21 Df22 Df23 Df24 Df25 Df18 Df19 Df20 Df26... %0th order
                Df27 Df28 Df29 Df30 Df30 Df31 Df32 Df33 Df34 Df34 Df35 Df36 Df37 Df38 Df39 Df40 Df41 Df42 Df43 Df44 Df37 Df38 ... %1st order
                Df45 Df46 Df47 Df48 Df48 Df49 Df50 Df51 Df52 Df52 Df53 Df54 Df55 Df56 Df56 Df57 Df58 Df59 Df60 Df61 Df62 Df63 Df64 Df65 Df66 Df57 Df58 Df59 Df67 Df60 ...
                Df68 Df69 Df70 Df71 Df72 Df73 Df74 Df75 Df76 Df77 Df78 Df79 Df80 Df70 Df71 Df81 ... %2nd order
                Df82 Df83 Df84 Df85 Df86 Df87 Df88 Df89 Df86 Df87 ... %3rd order
                Df90 Df91 Df92 Df93 Df94 Df95 Df96 Df97 Df98 Df99 Df100 Df101 Df102 Df95 Df96 Df97 Df103 Df102 Df98 ... %4th order
                ];
            
            Terms.DfE = Terms.fE/alph; %deriving alph*constant is the same as dividing by alph
            Terms.DfI = Terms.fI*(-2/alph); %deriving alph^(-2)*constant is the same as multiplying by (-2/alph)
            
            end
            
            
            
      
            
            %%
            Terms.j = ones(1,Terms.N)*j;
            
            if exist('AllTerms','var')
                AllTerms = CatStructFields(AllTerms,Terms,2);
            else
                AllTerms = Terms;
            end
            
        end
        
        %% one vectorized operation for all terms, including vectorizing on j
        TermsInd = find((AllTerms.j~=0 | AllTerms.k~=0) & (AllTerms.A+AllTerms.Ap+AllTerms.B+AllTerms.Bp)<=MaxPower);
        [dLambda1_pair,dLambda2_pair,dz1_pair,dz2_pair,du1_pair,du2_pair,da_oa1_pair,da_oa2_pair] = DisturbingFunctionMultiTermVariations(n1,n2,alph,Lambda1,Lambda2,mu1,mu2,AllTerms.j(TermsInd),AllTerms.k(TermsInd),e1,e2,pom1,pom2,I1,I2,s1,s2,Om1,Om2,AllTerms.A(TermsInd),AllTerms.Ap(TermsInd),AllTerms.B(TermsInd),AllTerms.Bp(TermsInd),AllTerms.C(TermsInd),AllTerms.Cp(TermsInd),AllTerms.D(TermsInd),AllTerms.Dp(TermsInd),AllTerms.f(TermsInd),AllTerms.fE(TermsInd),AllTerms.fI(TermsInd),AllTerms.Df(TermsInd),AllTerms.DfE(TermsInd),AllTerms.DfI(TermsInd));
        
        
        %% add values of dz and dLambda of the currently analyzed pair
        dLambdaMatrix(:,j1) = dLambdaMatrix(:,j1)+dLambda1_pair;
        dzMatrix(:,j1) = dzMatrix(:,j1)+dz1_pair;
        duMatrix(:,j1) = duMatrix(:,j1)+du1_pair;
        da_oaMatrix(:,j1) = da_oaMatrix(:,j1)+da_oa1_pair;
        
        dLambdaMatrix(:,j2) = dLambdaMatrix(:,j2)+dLambda2_pair;
        dzMatrix(:,j2) = dzMatrix(:,j2)+dz2_pair;
        duMatrix(:,j2) = duMatrix(:,j2)+du2_pair;
        da_oaMatrix(:,j2) = da_oaMatrix(:,j2)+da_oa2_pair;
    end
end

%%
dLambda = dLambdaMatrix;
dz = dzMatrix;
du = duMatrix;
da_oa = da_oaMatrix;







function [dLambda1_j,dz1_j,dLambda2_j,dz2_j] = VariationsPair(n1,n2,alph,Lambda1,Lambda2,mu1,mu2,j,z1,z2,e1,e2,pom1,pom2,u1,u2,s1,s2,Om1,Om2)
%This function is the algebraic "heart" of the code, as it calculates the
%integrals of the rhs of Lagrange's planetary equations to yield variations
%in the eccentricities and mean longitudes.

%the following two flags, for adding additional terms, remain here for
%debugging and backward compatibility purposes, though for operational usage they will remain hard-coded as 1.
AddSecondOrderIncTerms = 1;
AddThirdOrder = 1;


%calculate laplace coefficients
Aj = LaplaceCoeffs(alph,0.5,j,0);
DAj = LaplaceCoeffs(alph,0.5,j,1);
D2Aj = LaplaceCoeffs(alph,0.5,j,2);
D3Aj = LaplaceCoeffs(alph,0.5,j,3);

Ajp1 = LaplaceCoeffs(alph,0.5,j+1,0);
DAjp1 = LaplaceCoeffs(alph,0.5,j+1,1);
D2Ajp1 = LaplaceCoeffs(alph,0.5,j+1,2);
D3Ajp1 = LaplaceCoeffs(alph,0.5,j+1,3);

Ajm1 = LaplaceCoeffs(alph,0.5,j-1,0);
DAjm1 = LaplaceCoeffs(alph,0.5,j-1,1);
D2Ajm1 = LaplaceCoeffs(alph,0.5,j-1,2);
D3Ajm1 = LaplaceCoeffs(alph,0.5,j-1,3);

Ajm2 = LaplaceCoeffs(alph,0.5,j-2,0);
DAjm2 = LaplaceCoeffs(alph,0.5,j-2,1);
D2Ajm2 = LaplaceCoeffs(alph,0.5,j-2,2);
D3Ajm2 = LaplaceCoeffs(alph,0.5,j-2,3);

f1 = 0.5*Aj;
f2 = -0.5*j^2*Aj+0.25*alph*DAj+1/8*alph^2*D2Aj;
f10 = 0.25*(2+6*j+4*j^2)*Ajp1-0.5*alph*DAjp1-0.25*alph^2*D2Ajp1;

f27 = -j*Aj-0.5*alph*DAj;
f31 = (j-0.5)*Ajm1+0.5*alph*DAjm1;

f45 = 1/8*(-5*j+4*j^2)*Aj+0.5*j*alph*DAj+1/8*alph^2*D2Aj;
f49 = 0.25*(-2+6*j-4*j^2)*Ajm1+0.25*alph*(2-4*j)*DAjm1-0.25*alph^2*D2Ajm1;
f53 = 1/8*(2-7*j+4*j^2)*Ajm2+1/8*(-2*alph+4*j*alph)*DAjm2+1/8*alph^2*D2Ajm2;

Df1 = 0.5*DAj;
Df2 = -0.5*j^2*DAj+0.25*DAj+0.25*alph*D2Aj+1/4*alph*D2Aj+1/8*alph^2*D3Aj;
Df10 = 0.25*(2+6*j+4*j^2)*DAjp1-0.5*DAjp1-0.5*alph*D2Ajp1-0.5*alph*D2Ajp1-0.25*alph^2*D3Ajp1;

Df27 = -j*DAj-0.5*DAj-0.5*alph*D2Aj;
Df31 = (j-0.5)*DAjm1+0.5*DAjm1+0.5*alph*D2Ajm1;

Df45 = 1/8*(-5*j+4*j^2)*DAj+0.5*j*DAj+0.5*j*alph*D2Aj+1/4*alph*D2Aj+1/8*alph^2*D3Aj;
Df49 = 0.25*(-2+6*j-4*j^2)*DAjm1+0.25*(2-4*j)*DAjm1+0.25*alph*(2-4*j)*D2Ajm1-0.5*alph*D2Ajm1-0.25*alph^2*D3Ajm1;
Df53 = 1/8*(2-7*j+4*j^2)*DAjm2+1/8*(-2+4*j)*DAjm2+1/8*(-2*alph+4*j*alph)*D2Ajm2+1/4*alph*D2Ajm2+1/8*alph^2*D3Ajm2;

%second order inclination-related and third-order inclination related terms
Bjm1 = LaplaceCoeffs(alph,1.5,j-1,0);
DBjm1 = LaplaceCoeffs(alph,1.5,j-1,1);
f57 = 0.5*alph*Bjm1;
Df57 = 0.5*Bjm1+0.5*alph*DBjm1;
f62 = -2*f57;
Df62 = -2*Df57;

%define relative angle (as used at Hadden et al. 2019)
Psi = j*Lambda2-j*Lambda1;

%define resonant angles
Lambda_j1 = j*Lambda2+(1-j)*Lambda1;
Lambda_j2 = j*Lambda2+(2-j)*Lambda1;
Lambda_j3 = j*Lambda2+(3-j)*Lambda1;

%define "super mean motions"
nj1 = j*n2+(1-j)*n1;
nj2 = j*n2+(2-j)*n1;
nj3 = j*n2+(3-j)*n1;

%da/a term - inner planet
if j~=0
    Order0Term = ((f1-alph*(j==1)+(f2+0.5*alph*(j==1))*(e1.^2+e2.^2)).*sin(Psi)...
        +(f10-alph*(j==-2))*e1.*e2.*sin(Psi+pom2-pom1))/(-j)/((n2-n1).^2);
else
    Order0Term = 0;
end

Order1Term = (...
    (f27+3/2*alph*(j==1)-0.5*alph*(j==-1)).*e1.*sin(Lambda_j1-pom1)+...
    (f31-2*alph*(j==2)).*e2.*sin(Lambda_j1-pom2)...
    )...
    *(1-j)/((nj1).^2);

Order2Term = (...
    (f45-3/8*alph*(j==-1)-1/8*alph*(j==1))*e1.^2.*sin(Lambda_j2-2*pom1)+...
    (f49+3*alph*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)...
    +(f53-1/8*alph*(j==1)-27/8*alph*(j==3))*e2.^2.*sin(Lambda_j2-2*pom2)...
    )...
    .*(2-j)/((nj2).^2);

Int_da_over_a_term = -3*alph*n1^2*(mu2/(1+mu1))*(Order0Term+Order1Term+Order2Term);

%dR/dalpha terms - inner planet
if j~=0
    Order0Term = (...
        (Df1-(j==1)+(Df2+0.5*(j==1))*(e1.^2+e2.^2)).*sin(Psi)...
        +(Df10-(j==-2))*e1.*e2.*sin(Psi+pom2-pom1))/j./(n2-n1);
else
    Order0Term = 0;
end

Order1Term = (...
    (Df27+3/2*(j==1)-0.5*(j==-1)).*e1.*sin(Lambda_j1-pom1)...
    +(Df31-2*(j==2)).*e2.*sin(Lambda_j1-pom2)...
    )...
    ./(nj1);

Order2Term = (...
    (Df45-3/8*(j==-1)-1/8*(j==1)).*e1.^2.*sin(Lambda_j2-2*pom1)...
    +(Df49+3*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)...
    +(Df53-1/8*(j==1)-27/8*(j==3)).*e2.^2.*sin(Lambda_j2-2*pom2)...
    )...
    ./(nj2);

Int_dR_dalph_term = -2*alph^2*n1*(mu2/(1+mu1))*(Order0Term+Order1Term+Order2Term);

%e-derivative term - inner planet
if j~=0
    Order0Term = (...
        (2*f2+alph*(j==1))*e1.*sin(Psi)+...
        (f10-alph*(j==-2))*e2.*sin(Psi-pom1+pom2)...
        )...
        ./(j*(n2-n1));
else
    Order0Term = 0;
end

Order1Term =  (f27+3/2*alph*(j==1)-0.5*alph*(j==-1))*sin(Lambda_j1-pom1)./(nj1);

Order2Term = (...
    2*(f45-3/8*alph*(j==-1)-1/8*alph*(j==1))*e1.*sin(Lambda_j2-2*pom1)+...
    (f49+3*alph*(j==2))*e2.*sin(Lambda_j2-pom1-pom2)...
    )...
    ./(nj2);

e_term = alph/2.*e1.*n1*(mu2/(1+mu1)).*(Order0Term+Order1Term+Order2Term);


dLambda1_j = Int_da_over_a_term+Int_dR_dalph_term+e_term;

if AddSecondOrderIncTerms
    %inclination terms of the inner planet
    dl1Order2IncTerm1 = 3*mu2/(1+mu1)*alph*n1^2/(nj2)^2*(j-2)*(f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+f57*s2.^2.*sin(Lambda_j2-2*Om2));
    dl1Order2IncTerm2 = -2*mu2/(1+mu1)*alph^2*n1/(nj2)*(Df57*s1.^2.*sin(Lambda_j2-2*Om1)+Df62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+Df57*s2.^2.*sin(Lambda_j2-2*Om2));
    dl1Order2IncTerm3 = mu2/(1+mu1)*alph*n1/nj2*(2*f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2));
    dLambda1_j = dLambda1_j+dl1Order2IncTerm1+dl1Order2IncTerm2+dl1Order2IncTerm3;
end

%z variations - inner planet
if j~=0
    Order0Term = (...
        ((2*f2+alph*(j==1)).*sin(Psi)+1i*j*(0.5*f1-0.5*alph*(j==1))*(-1)*cos(Psi)).*z1...
        +(f10-alph*(j==-2))*exp(1i*(Psi)).*z2/1i...
        )...
        ./(j*(n2-n1));
else
    Order0Term = 0;
end

Order1Term =  (...
    (f27+3/2*alph*(j==1)-0.5*alph*(j==-1))...
    *(...
    (1-0.5*e1.^2).*exp(1i*(Lambda_j1))./(1i)...
    +0.5i*(j-1).*e1.^2.*exp(1i*pom1).*(-1).*cos(Lambda_j1-pom1)...
    )...
    +...
    0.5i*(f31-2*alph*(j==2))*(j-1)*e1.*e2.*exp(1i*pom1).*(-1).*cos(Lambda_j1-pom2)...
    )...
    ./(nj1);

Order2Term = (...
    2*(f45-3/8*alph*(j==-1)-1/8*alph*(j==1)).*exp(1i*(Lambda_j2))/1i.*conj(z1)+...
    (f49+3*alph*(j==2)).*exp(1i*(Lambda_j2))/1i.*conj(z2)...
    )...
    ./(nj2);


dz1_j = 1i*n1*alph*mu2/(1+mu1).*(Order0Term+Order1Term+Order2Term);


%da'/a' term - outer planet
if j~=0
    Order0Term = (...
        -(f1-1/alph^2*(j==1)+(f2+0.5/alph^2*(j==1))*(e1.^2+e2.^2)).*sin(Psi)...
        -(f10-1/alph^2*(j==-2))*e1.*e2.*sin(Psi+pom2-pom1)...
        )...
        /(-j)/((n2-n1).^2);
else
    Order0Term = 0;
end

Order1Term = (...
    -(f27-2/alph^2*(j==-1))*e1.*sin(Lambda_j1-pom1)...
    -(f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2)).*e2.*sin(Lambda_j1-pom2)...
    )...
    *(-j)/((nj1).^2);

Order2Term = (...
    -(f45-27/8/alph^2*(j==-1)-1/8/alph^2*(j==1))*e1.^2.*sin(Lambda_j2-2*pom1)...
    -(f49+3/alph^2*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)...
    -(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3))*e2.^2.*sin(Lambda_j2-2*pom2)...
    )...
    .*(-j)./((nj2).^2);

Int_da_over_a_term = -3*n2^2*(mu1/(1+mu2)).*(Order0Term+Order1Term+Order2Term);

%dR'/dalpha terms - outer planet - corrected
if j~=0
    Order0Term = (...
        (f1+alph*Df1+(j==1)/alph^2+(f2+alph*Df2-(j==1)/(2*alph^2)).*(e1.^2+e2.^2)).*sin(Psi)...
        +(f10+alph*Df10+(j==-2)/alph^2).*e1.*e2.*sin(Psi+pom2-pom1)...
        )...
        /j./(n2-n1);
else
    Order0Term = 0;
end

Order1Term = (...
    (f27+alph*Df27+2/(alph^2)*(j==-1)).*e1.*sin(Lambda_j1-pom1)...
    +(f31+alph*Df31-3/alph^2*(j==0)+1/(2*alph^2)*(j==2)).*e2.*sin(Lambda_j1-pom2)...
    )...
    ./(nj1);

Order2Term = (...
    (f45+alph*Df45+27/8/alph^2*(j==-1)+1/8/alph^2*(j==1)).*e1.^2.*sin(Lambda_j2-2*pom1)...
    +(f49+alph*Df49-3/alph^2*(j==2)).*e1.*e2.*sin(Lambda_j2-pom1-pom2)...
    +(f53+alph*Df53+1/8/alph^2*(j==1)+3/8/alph^2*(j==3)).*e2.^2.*sin(Lambda_j2-2*pom2)...
    )...
    ./(nj2);


Int_dR_dalph_term = +2*n2*(mu1/(1+mu2))*(Order0Term+Order1Term+Order2Term);

%e'-derivative term - outer planet
if j~=0
    Order0Term = (...
        (2*f2+1/alph^2*(j==1))*e2.*sin(Psi)+...
        (f10-1/alph^2*(j==-2))*e1.*sin(Psi+pom2-pom1)...
        )...
        ./(j*(n2-n1));
else
    Order0Term = 0;
end

Order1Term =  (f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2))*sin(Lambda_j1-pom2)...
    ./(nj1);

Order2Term = (...
    (f49+3/alph^2*(j==2))*e1.*sin(Lambda_j2-pom1-pom2)+...
    2*(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3))*e2.*sin(Lambda_j2-2*pom2)...
    )...
    ./(nj2);

e_term = e2/2*n2*(mu1/(1+mu2)).*(Order0Term+Order1Term+Order2Term);

dLambda2_j = (Int_da_over_a_term+Int_dR_dalph_term+e_term);


if AddSecondOrderIncTerms
    %inclination terms of the outer planet
    dl2Order2IncTerm1 = -3*mu1/(1+mu2)*n2^2/(nj2)^2*(j)*(f57*s1.^2.*sin(Lambda_j2-2*Om1)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2)+f57*s2.^2.*sin(Lambda_j2-2*Om2));
    dl2Order2IncTerm2 = +2*mu1/(1+mu2)*n2/(nj2)*((f57+alph*Df57)*s1.^2.*sin(Lambda_j2-2*Om1)+(f62+alph*Df62)*s1.*s2.*sin(Lambda_j2-Om1-Om2)+(f57+alph*Df57)*s2.^2.*sin(Lambda_j2-2*Om2));
    dl2Order2IncTerm3 = mu1/(1+mu2)*n2/nj2*(2*f57*s2.^2.*sin(Lambda_j2-2*Om2)+f62*s1.*s2.*sin(Lambda_j2-Om1-Om2));
    dLambda2_j = dLambda2_j+dl2Order2IncTerm1+dl2Order2IncTerm2+dl2Order2IncTerm3;
end

%z' variations - outer planet
if j~=0
    Order0Term = (...
        ((2*f2+1/alph^2*(j==1)).*sin(Psi)+1i*j*(-0.5*f1+0.5/alph^2*(j==1))*(-1)*cos(Psi)).*z2...
        +(f10-1/alph^2*(j==-2))*exp(-1i*j*(Lambda2-Lambda1))/(-1i).*z1...
        )...
        ./(j*(n2-n1));
else
    Order0Term = 0;
end

Order1Term =  (...
    (f31+3/2/alph^2*(j==0)-0.5/alph^2*(j==2))*(1-0.5*e2.^2).*exp(1i*(Lambda_j1))./(1i)...
    -0.5i*j*(f27-1/alph^2*(j==-1)).*e1.*e2.*exp(1i*pom2).*(-1).*cos(Lambda_j1-pom1)...
    -0.5i*j*(f31+3/2/alph^2*(j==0)-1/2/alph^2*(j==2))*e2.^2.*exp(1i*pom2).*(-1).*cos(Lambda_j1-pom2)...
    )...
    ./(nj1);

Order2Term = (...
    (f49+3/alph^2*(j==2))*exp(1i*(Lambda_j2))/1i.*conj(z1)+...
    2*(f53-1/8/alph^2*(j==1)-3/8/alph^2*(j==3)).*exp(1i*(Lambda_j2))/1i.*conj(z2)...
    )...
    ./(nj2);


dz2_j = 1i*n2*mu1/(1+mu2)*(Order0Term+Order1Term+Order2Term);




if AddThirdOrder
    
    %calculate Laplace Coefficients required for 3rd order resonant terms
    Ajm3 = LaplaceCoeffs(alph,0.5,j-3,0);
    DAjm3 = LaplaceCoeffs(alph,0.5,j-3,1);
    D2Ajm3 = LaplaceCoeffs(alph,0.5,j-3,2);
    D3Ajm3 = LaplaceCoeffs(alph,0.5,j-3,3);
    
    D4Aj = LaplaceCoeffs(alph,0.5,j,4);
    D4Ajm1 = LaplaceCoeffs(alph,0.5,j-1,4);
    D4Ajm2 = LaplaceCoeffs(alph,0.5,j-2,4);
    D4Ajm3 = LaplaceCoeffs(alph,0.5,j-3,4);
    
    %calculate the coefficients from Murray&Dermott appendix B related to
    %3rd order resonant terms
    f82 = 1/48*(-26*j*Aj+30*j^2*Aj-8*j^3*Aj-9*alph*DAj+27*j*alph*DAj-12*j^2*alph*DAj+6*alph^2*D2Aj-6*j*alph^2*D2Aj-alph^3*D3Aj);
    f83 = 1/16*(-9*Ajm1+31*j*Ajm1-30*j^2*Ajm1+8*j^3*Ajm1+9*alph*DAjm1-25*j*alph*DAjm1+12*j^2*alph*DAjm1-5*alph^2*D2Ajm1+6*j*alph^2*D2Ajm1+alph^3*D3Ajm1);
    f84 = 1/16*(8*Ajm2-32*j*Ajm2+30*j^2*Ajm2-8*j^3*Ajm2-8*alph*DAjm2+23*j*alph*DAjm2-12*j^2*alph*DAjm2+4*alph^2*D2Ajm2-6*j*alph^2*D2Ajm2-alph^3*D3Ajm2);
    f85 = 1/48*(-6*Ajm3+29*j*Ajm3-30*j^2*Ajm3+8*j^3*Ajm3+6*alph*DAjm3-21*j*alph*DAjm3+12*j^2*alph*DAjm3-3*alph^2*D2Ajm3+6*j*alph^2*D2Ajm3+alph^3*D3Ajm3);
    
    Df82 = 1/48*(-26*j*DAj+30*j^2*DAj-8*j^3*DAj-9*DAj-9*alph*D2Aj+27*j*DAj+27*j*alph*D2Aj-12*j^2*DAj-12*j^2*alph*D2Aj+12*alph*D2Aj+6*alph^2*D3Aj-12*j*alph*D2Aj-6*j*alph^2*D3Aj-3*alph^2*D3Aj-alph^3*D4Aj);
    Df83 = 1/16*(-9*DAjm1+31*j*DAjm1-30*j^2*DAjm1+8*j^3*DAjm1+9*DAjm1+9*alph*D2Ajm1-25*j*DAjm1-25*j*alph*D2Ajm1+12*j^2*DAjm1+12*j^2*alph*D2Ajm1-10*alph*D2Ajm1-5*alph^2*D3Ajm1+12*j*alph*D2Ajm1+6*j*alph^2*D3Ajm1+3*alph^2*D3Ajm1+alph^3*D4Ajm1);
    Df84 = 1/16*(8*DAjm2-32*j*DAjm2+30*j^2*DAjm2-8*j^3*DAjm2-8*DAjm2-8*alph*D2Ajm2+23*j*DAjm2+23*j*alph*D2Ajm2-12*j^2*DAjm2-12*j^2*alph*D2Ajm2+8*alph*D2Ajm2+4*alph^2*D3Ajm2-12*j*alph*D2Ajm2-6*j*alph^2*D3Ajm2-3*alph^2*D3Ajm2-alph^3*D4Ajm2);
    Df85 = 1/48*(-6*DAjm3+29*j*DAjm3-30*j^2*DAjm3+8*j^3*DAjm3+6*DAjm3+6*alph*D2Ajm3-21*j*DAjm3-21*j*alph*D2Ajm3+12*j^2*DAjm3+12*j^2*alph*D2Ajm3-6*alph*D2Ajm3-3*alph^2*D3Ajm3+12*j*alph*D2Ajm3+6*j*alph^2*D3Ajm3+3*alph^2*D3Ajm3+alph^3*D4Ajm3);
    
    
    D2Bjm1 = LaplaceCoeffs(alph,1.5,j-1,2);
    Bjm2 = LaplaceCoeffs(alph,1.5,j-2,0);
    DBjm2 = LaplaceCoeffs(alph,1.5,j-2,1);
    D2Bjm2 = LaplaceCoeffs(alph,1.5,j-2,2);
    
    
    f86 = 0.75*alph*Bjm1-0.5*j*alph*Bjm1-0.25*alph^2*DBjm1;
    f87 = 0.5*j*alph*Bjm2+0.25*alph^2*DBjm2;
    f88 = -2*f86;
    f89 = -2*f87;
    Df86 = 0.75*Bjm1+0.75*alph*DBjm1  -0.5*j*Bjm1-0.5*j*alph*DBjm1   -0.5*alph*DBjm1-0.25*alph^2*D2Bjm1;
    Df87 = 0.5*j*Bjm2+0.5*j*alph*DBjm2+   0.5*alph*DBjm2+0.25*alph^2*D2Bjm2;
    Df88 = -2*Df86;
    Df89 = -2*Df87;
    
    %calculate the 3rd order resonant terms contribution to the forced eccentricity
    dz1Order3Term1 = mu2/(1+mu1)*alph*n1/(nj3)*exp(1i*pom1).*(3*f82*e1.^2.*exp(1i*(Lambda_j3-3*pom1))          +2*f83.*e1.*e2.*exp(1i*(Lambda_j3-2*pom1-pom2))+1*f84*e2.^2.*exp(1i*(Lambda_j3-pom1-2*pom2)));
    dz2Order3Term1 = mu1/(1+mu2)        *n2/(nj3)*exp(1i*pom2).*(1*f83*e1.^2.*exp(1i*(Lambda_j3-2*pom1-pom2))+2*f84*e1.*e2.*exp(1i*(Lambda_j3-pom1-2*pom2))+3*(f85-(j==4)/(3*alph^2))*e2.^2.*exp(1i*(Lambda_j3-3*pom2)));
    
    %add the contribution of the terms involving both e and I (or s in
    %Murray&Dermott formulae) - these with the coefficients f86-f89
    dz1Order3Term2 = mu2/(1+mu1)*alph*n1/(nj3).*(f86*s1.^2.*exp(1i*(Lambda_j3-2*Om1))+f88*s1.*s2.*exp(1i*(Lambda_j3-Om1-Om2))+f86*s2.^2.*exp(1i*(Lambda_j3-2*Om2)));
    dz2Order3Term2 = mu1/(1+mu2)        *n2/(nj3).*(f87*s1.^2.*exp(1i*(Lambda_j3-2*Om1))+f89*s1.*s2.*exp(1i*(Lambda_j3-Om1-Om2))+f87*s2.^2.*exp(1i*(Lambda_j3-2*Om2)));
    
    %calculate the 3rd order resonant terms contribution to the forced
    %mean longitude
    
    %dn/n contribution
    dl1Order3Term1 = 3*(j-3)*mu2/(1+mu1)*alph*n1^2/(nj3)^2*(f82*e1.^3.*sin(Lambda_j3-3*pom1)+f83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+f84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85-(j==4)*16/3*alph)*e2.^3.*sin(Lambda_j3-3*pom2));
    dl2Order3Term1 = 3*(-j)*mu1/(1+mu2)*n2^2/(nj3)^2*(f82*e1.^3.*sin(Lambda_j3-3*pom1)+f83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+f84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85-(j==4)/(3*alph^2))*e2.^3.*sin(Lambda_j3-3*pom2));
    
    %dR/da and dR/da' contribution
    dl1Order3Term2 = -2*mu2/(1+mu1)*alph^2*n1/(nj3)*(Df82*e1.^3.*sin(Lambda_j3-3*pom1)+Df83*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+Df84*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(Df85-16/3*(j==4))*e2.^3.*sin(Lambda_j3-3*pom2));
    dl2Order3Term2 = 2*mu1/(1+mu2)*n2/(nj3)*((f82+alph*Df82)*e1.^3.*sin(Lambda_j3-3*pom1)+(f83+alph*Df83)*e1.^2.*e2.*sin(Lambda_j3-2*pom1-pom2)+(f84+alph*Df84)*e1.*e2.^2.*sin(Lambda_j3-pom1-2*pom2)+(f85+alph*Df85+2/(3*alph^3)*(j==4))*e2.^3.*sin(Lambda_j3-3*pom2));
    
    %add the contribution of the terms involving both e and I (or s in
    %Murray&Dermott formulae) - these with the coefficients f86-f89
    
    %dn/n contribution
    dl1Order3Term3 = 3*(j-3)*mu2/(1+mu1)*alph*n1^2/(nj3)^2*(f86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) +f87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +f88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +f89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + f86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) +f87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2));
    dl2Order3Term3 = 3*(-j)  *mu1/(1+mu2)*n2^2      ./(nj3)^2*(f86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) +f87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +f88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +f89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + f86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) +f87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2));
    
    %dR/da and dR/da' contribution
    dl1Order3Term4 = -2*mu2/(1+mu1)*alph^2*n1/(nj3)* (Df86*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) + Df87*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +Df88*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +Df89*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + Df86*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) + Df87*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2));
    dl2Order3Term4 = 2*mu1/(1+mu2)*n2/(nj3)*              ((f86+alph*Df86)*e1.*s1.^2.*sin(Lambda_j3-pom1-2*Om1) + (f87+alph*Df87)*e2.*s1.^2.*sin(Lambda_j3-pom2-2*Om1) +(f88+alph*Df88)*e1.*s1.*s2.*sin(Lambda_j3-pom1-Om1-Om2) +(f89+alph*Df89)*e2.*s1.*s2.*sin(Lambda_j3-pom2-Om1-Om2) + (f86+alph*Df86)*e1.*s2.^2.*sin(Lambda_j3-pom1-2*Om2) + (f87+alph*Df87)*e2.*s2.^2.*sin(Lambda_j3-pom2-2*Om2));
    
    %sum all the contributions to delta-lambda
    dl1Order3 = dl1Order3Term1+dl1Order3Term2+dl1Order3Term3+dl1Order3Term4;
    dl2Order3 = dl2Order3Term1+dl2Order3Term2+dl2Order3Term3+dl2Order3Term4;
    
    %sum all the contributions to delta-z
    dz1Order3 = dz1Order3Term1+dz1Order3Term2;
    dz2Order3 = dz2Order3Term1+dz2Order3Term2;
    
    %add the 3rd order contribution to the total forced eccentricities and
    %mean longitudes
    dLambda1_j = dLambda1_j+dl1Order3;
    dz1_j = dz1_j+dz1Order3;
    dLambda2_j = dLambda2_j+dl2Order3;
    dz2_j = dz2_j+dz2Order3;
    
end