function set6()
clear all
%close all
% this program uses a predictor-corrector method (6th order Milne's method)
% to solve the system y'=f(x,y). Initial points and points needed in halving 
% step size are supplied by a 4th order Runge-Kutta scheme. Enough function
% values are retained at each step so that the step size can be doubled
% without generating any new points. The Runge-Kutta scheme is in function
% RK4 and the function dy evaluates f(x,y)
% array size declarations
A  = zeros(4,4); PHI = zeros(1,4);
F  = zeros(1,4); Y = F; YY = Y; Y1 = Y; Y2 = Y; Y3 = Y; ZZ = Y; U = Y; US = U;
Z  = zeros(5,4); 
F1 = zeros(7,4); U1 = F1; UB = U1; FB = F1; 
EX = zeros(3000); EY = EX; 
DOC = zeros(3);
% Milne parameters
AL  = 33/10; BE  = -21/5 ; GA  = 39/5 ;
AL1 = 14/45; BE1 =  64/45; GA1 = 24/45;
N5 = 0; %index of progress
DOC(1) = 0; DOC(3) = 1;
NNN = 0; %index for doubler array: need 11 saved values
NNNN = 0; %index for printing
H = .001
fprintf('\n     X(T)          Y(T)           T    ')
N  = 4; % problem size
NN = 5; % scheme order P-1
EP = .00000005;
% set initial parameters
mu = 0.012277471;
% Initialization. Make sure initial time step allows desired 
% accuracy. If not, keep halving step until accuracy is achieved.
    % Initial values
    Y(1) = .994; Y(2:3) = 0; Y(4) = -2.031732629557;
    %Y(1) = .994; Y(2) = 0; Y(3) = 0; Y(4) = -2.113898796645;
    % period
    PERIOD  = 11.12402;
    %PERIOD  =  5.437;
    T0 = 0; %-------------------------------------------------
    ZZ = DP(T0,Y); 
    TEND = 4*PERIOD;
    N5 = 1;  
    S = 10*EP; H = 2*H; % kludge to correct initialization of RK process
    DX(1) = Y(1); DY(1) = Y(3);
while T0 < TEND
    while S > EP  % take 5+1 steps from present and keep halving until done
        H = .5*H;
        T1 = T0; N51 = N5; NNNN = 0;%-------------------------
        U1(1,:) =  Y(:); F1(1,:) = ZZ(:); 
        % get 5 points by Runge-Kutta
        %NN = floor(PERIOD/H)+1; % setup to integrate with fixed-step RK5
        for k = 1:NN
            k1 = k + 1; N51 = N51 + 1;     
            PHI(1,:) = RK5(T1,U1(k,:),H);
            U1(k1,:) = U1(k,:) + H*PHI(1,:);
            T1 = T1 + H;%------------------------------------- 
            F1(k1,:) = DP(T1,U1(k1,:));
            DX(N51) = U1(k1,1);                                                %
            DY(N51) = U1(k1,3);                                                %
        end
        % Predictor                                                          %
        US=U1(1,:)+H*(AL*(F1(2,:)+F1(6,:))+BE*(F1(3,:)+F1(5,:))+GA*F1(4,:));  %
        Z1 = DP(T0,US);                                                       %
        % Corrector 
        U=U1(3,:)+H*(BE1*(F1(6,:)+F1(4,:))+GA1*F1(5,:)+AL1*(F1(3,:)+Z1));     %
        % check accuracy (U - U*)
        S = 0;
        for i = 1:N                                                          %
            S = S + abs(U(i) - US(i));                                       %
        end
        k1 = k1+1; N51 = N51 + 1; T1 = T1 + H; %---------------
        UB(6, :) = U1(1, :); FB(6, :) = F1(1, :);                            %
        for kk = 1:NN                                                        %
            kk1 = kk + 1;                                                    %
            U1(kk, :) = U1(kk1, :); F1(kk, :) = F1(kk1, :);                  %
        end                                                                  %                                                    %
        U1(6, :) =  U; F1(6,:) = DP(T1,U);       
        DX(N51) = U(1);                                                      %
        DY(N51) = U(3); 
        NNN = NN;
    end 
    T0 = T1;%--------------------------------------------------
    N5 = N51;
    %------------------------------------------------------------------- 45
    % Predictor                                                           %
    US=U1(1,:)+H*(AL*(F1(2,:)+F1(6,:))+BE*(F1(3,:)+F1(5,:))+GA*F1(4,:));  %
    Z1 = DP(T0+H,US);                                                       %
    % Corrector 
    U=U1(3,:)+H*(BE1*(F1(6,:)+F1(4,:))+GA1*F1(5,:)+AL1*(F1(3,:)+Z1));     %
    % check accuracy (U - U*)
    S = 0;                                                                %
    for i = 1:N                                                           %
        S = S + abs(U(i) - US(i));                                        %
    end
    if (S > EP) 
        % prepare to halve the stepsize
        Y = U1(3, :); ZZ = F1(3, :);
        T0 = T0 - 3*H; N5 = N5 - 3;%----------------------------
    else
        N5 = N5 + 1;        
        T0 = T0 + H; %------------------------------------------
        ZZ = DP(T0, U);
        DX(N5) = U(1);                                                       %
        DY(N5) = U(3); 
        NNNN = NNNN + 1; 
        T1 = T0;  N51 = N5;%------------------------------------
        if (NNNN == 5)                                                       %
            fprintf('\n %1.15f    %1.15f     %1.15f    %6.6i',U(1), U(2), T0 , N5)    %
            NNNN = 0;                                                        %
        end            
        if (S < .1*EP & NNN >= 11)
            H = 2*H;
            U1(6, :) = U;        F1(6, :) = ZZ;
            U1(4, :) = U1(3, :); F1(4, :) = F1(3, :);
            U1(3, :) = U1(1, :); F1(3, :) = F1(1, :);
            U1(2, :) = UB(5, :); F1(2, :) = FB(5, :);
            U1(1, :) = UB(3, :); F1(1, :) = FB(3, :);
            UB(6, :) = UB(1, :); FB(6, :) = FB(1, :);
            NNN = 6;
        else%no change in stepsize
            NNN = NNN + 1;                                                       %
            for kk = 1:NN                                                        %
                kk1 = kk + 1;                                                    %
                FB(kk, :) = FB(kk1, :);                                          %
                UB(kk, :) = UB(kk1, :);                                          %
            end                                                                  %
            FB(6, :) = F1(1, :);                                                 %
            UB(6, :) = U1(1, :);                                                 %
            for kk = 1:NN                                                        %
                kk1 = kk + 1;                                                    %
                F1(kk, :) = F1(kk1, :);                                          %
                U1(kk, :) = U1(kk1, :);                                          %
            end                                                                  %                                                    %
            U1(6, :) =  U;                                                       %
            F1(6, :) = ZZ;                                                       %                         
        end                                                 
    end
end
% -------------goto 45---------------------------------------------------%
        % 2000
% this section produces the graph
XMAX = 0;XMIN = 0; YMAX = 0; YMIN = 0;
for k = 1:N5
    if ( DX(k) > XMAX)
        XMAX = DX(k);
    end
    if (DX(k) < XMIN)
        XMIN = DX(k);
    end
    if ( DY(k) > YMAX)
        YMAX = DY(k);
    end
    if (DY(k) < YMIN)
        YMIN = DY(k);
    end
end
%N5 = N5*2;
figure
hold on
plot(DX(1:N5), DY(1:N5), 'r-')
for i = 1:4
    plot(DX(floor(i*N5/4)),DY(floor(i*N5/4)), 'b*')
end
hold off
end
%==========================================================================
function YP = DP(T, Y)
mu = 0.012277471; m1 = 0.987722529;
xm = Y(1) + mu; xm1 = xm-1;y2 = Y(3)^2;
R1 = (xm^2 + y2)^(1.5); R2 = ((xm1)^2+y2)^(1.5);
YP(1) = Y(2);
YP(2) = Y(1) + 2*Y(4) - m1*xm  /R1 - mu*xm1 /R2;
YP(3) = Y(4);
YP(4) = Y(3) - 2*Y(2) - (m1/R1 + mu/R2)*Y(3);
end
%==========================================================================
function phi = RK5(x,u,h)
k1 = DP(x,u); h1 = h/4;
k2 = DP(x+h1, u+h1*k1); h2 = 3*h/8;
k3 = DP(x+h2, u+(3*h/32)*k1+(9*h/32)*k2); h3 = 12*h/13;
k4 = DP(x+h3, u+(1932*h/2197)*k1-(7200*h/2197)*k2+(7296*h/2197)*k3);
k5 = DP(x+h , u+(439*h/216)*k1-8*h*k2 + (3680*h/513)*k3 - (845*h/4104)*k4);
%phi(1:4,1) = (25/216)*k1 + (1408/2565)*k3 + (2197/4104)*k4 - (1/5)*k5;
k6=DP(x+h/2,u-(8*h/27)*k1+2*h*k2-(3544*h/2565)*k3+(1859*h/4104)*k4-(11*h/40)*k5);
phi(1:4)=(16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6;
end