% The following script will produce a fourth order Adams-Bashforth method.
% It will begin by finding the first four points by means of a 5th order
% Runge-Kutta scheme. Once those values are found it will carry on with A-B
% with step-size h.

function y = adamsbash4(ode_fxns,Tspan,IV,h)
[t_RK,y_RK] = get_points(ode_fxns,Tspan,IV,h);
t_temp = t_RK;
y_temp = y_RK;
y_check = zeros(4,1);
i = 1;j=0;
    while (t_temp(4) < Tspan(2))
    y(:,i+1) = y_temp(:,5)+ h/24*(55*ode_fxns(t_temp(4),y_temp(:,4))-59*...
        ode_fxns(t_temp(3),y_temp(:,3))+37*ode_fxns(t_temp(2),y_temp(:,2))...
        -9*ode_fxns(t_temp(1),y_temp(:,1)));
    
%     if (mod(t_temp(4),11.124/2) == 0)   
%         j=j+1;
%         if (j == 1)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'red');
%         elseif(j ==2)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'blue');
%         elseif(j==3)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'Color',[0 0.5 0]);
%         elseif(j==4)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'magenta');
%         elseif(j==5)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'cyan');
%         elseif(j==6)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'Color',[.5 0 0]);
%         elseif(j==7)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'black');
%         elseif(j==8)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'yellow');    
%         else
%             legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','j=8')
%         end
%         t_start = i+1;
%     end

%     y_check(:) = h/24*(55*ode_fxns(t_temp(1),y_temp(:,1))-59*ode_fxns(t_temp(2),y_temp(:,2))...
%         +37*ode_fxns(t_temp(3),y_temp(:,3))-9*ode_fxns(t_temp(4),y_temp(:,4)));
%     
%     abs(y(3:4,i+1)-y_check(3:4));
    y_temp(:,1:4) = y_temp(:,2:5); t_temp(1:3)= t_temp(2:4);
    y_temp(:,5) = y(:,i+1); t_temp(4) = t_temp(3) + h;
    i = i+1;
    end
    
    
end

function [t,y] = get_points(ode_fxns,Tspan,IV,h)
num_points = 4;
a = Tspan(1);
t = linspace(a,a+3*h,4);
y = zeros(4,1);
Phi = zeros(4,1);
Phi_star = zeros(4,1);
y(:,1) = IV;
for i = 1:num_points
    %----------------------------------------------------------------- 
    k_1 = ode_fxns(t(i),y(:,i));
    k_2 = ode_fxns(t(i)+2/7*h, y(:,i)+2/7*h.*k_1);
    k_3 = ode_fxns(t(i)+7/15*h,y(:,i)+77/900*h.*k_1+343/900*h.*k_2);
    k_4 = ode_fxns(t(i)+35/38*h,y(:,i)+805/1444*h.*k_1...
                   -77175/54872*h.*k_2+97125/54872*h.*k_3);
    Phi(:,i) = (79/490.*k_1+2175/3626.*k_3+2166/9065.*k_4);
    %-----------------------------------------------------------------
    k_5 = ode_fxns(t(i)+h,y(:,i)+h*Phi(:,i));
    Phi_star(:,i) = 229/1470*k_1+1125/1813*k_3+13718/81585*k_4+1/18*k_5;
    y(:,i+1) = y(:,i) + h*Phi_star(:,i);
   
end
    
end