% Declare the relevant parameters
h=0.1;
tmin=0.0;
tmax=1.0;

alpha = 0.0;

nmax=(tmax-tmin)/h;

t(1) = tmin;
w_rk4(1) = alpha;
w_ab(1) = alpha;
w_am(1) = alpha;
y(1) = alpha;

% Beginning of Loop 1

for i=2:4

% Advance time
    t(i) = t(i-1) + h;

% Runge Kutta 4
    k1 = h*deriv(t(i-1),w_rk4(i-1));
    k2 = h*deriv(t(i-1)+h/2,w_rk4(i-1)+k1/2.0);
    k3 = h*deriv(t(i-1)+h/2,w_rk4(i-1)+k2/2.0);
    k4 = h*deriv(t(i-1)+h,w_rk4(i-1)+k3);
    w_rk4(i)=w_rk4(i-1)+(k1 + 2*k2 + 2*k3 + k4)/6.0;

% Get first 3 values for Adam Bashforth Moulton
    w_ab(i) =  w_rk4(i);
    w_am(i) =  w_rk4(i);

% Exact Solution
    y(i) = exact(t(i));
    d_rk4(i) = abs(y(i) - w_rk4(i));
    d_am(i) = abs(y(i) - w_am(i));

end


% Beginning of Main Loop

for i=5:(nmax+1)

% Advance time
    t(i) = t(i-1) + h;

% Adam Bashforth Step (Predictor)
  w_ab(i) = w_am(i-1) + h/24 * (55*deriv(t(i-1),w_am(i-1)) - 59*deriv(t(i-2),w_am(i-2)));
  w_ab(i) = w_ab(i) + h/24 * (37*deriv(t(i-3),w_am(i-3)) - 9*deriv(t(i-4),w_am(i-4)));

% Adam Moulton Step (Corrector)
  w_am(i) = w_am(i-1) + h/24 * (9*deriv(t(i),w_ab(i)) +19*deriv(t(i-1),w_am(i-1)));
  w_am(i) = w_am(i) + h/24 * (-5*deriv(t(i-2),w_am(i-2))+deriv(t(i-3),w_am(i-3)));

% Runge Kutta 4
    k1 = h*deriv(t(i-1),w_rk4(i-1));
    k2 = h*deriv(t(i-1)+h/2,w_rk4(i-1)+k1/2.0);
    k3 = h*deriv(t(i-1)+h/2,w_rk4(i-1)+k2/2.0);
    k4 = h*deriv(t(i-1)+h,w_rk4(i-1)+k3);
    w_rk4(i)=w_rk4(i-1)+(k1 + 2*k2 + 2*k3 + k4)/6.0;

% Exact Solution
    y(i) = exact(t(i));
    d_rk4(i) = abs(y(i) - w_rk4(i));
    d_am(i) = abs(y(i) - w_am(i));

end 

% Exact Solution, fine grid

  tlast=tmin;
  y_fine(1) = alpha;
  h=h/10

% Beginning of Loop

  for i=2:(10*nmax+1)
    t_fine(i) = tlast + h;
    y_fine(i) = exact(t_fine(i));
    tlast=t_fine(i);
  end



% Plot Results for w, y

plot(t,w_rk4,'xm');
hold on;
plot(t,w_am,'og');
hold on;
plot(t_fine,y_fine,'r');

% Label the plot 

title('Solution to dy/dt=t*exp(3*t)-2y');
xlabel('Time');
ylabel('y');
legend('RK4','adam bashforth moulton','exact')

% Save the results (for each timestep) in a file

file = zeros(size(t,2),6);

file(:,1)=t';
file(:,2)=y';
file(:,3)=w_rk4';
file(:,4)=d_rk4';
file(:,5)=w_am';
file(:,6)=d_am';

save output_pc file -ASCII;
