% AB3_withRK_initialization.m
% MST 383/683 
% Fall 2020
%
% AB3 with Runge-Kutta initialization
% Approximates solution to y' = -y + 2cos(t), y(0)= 1, for t in [0,10]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear,clc

% Defining preliminary variables
h = 0.1; 
t = 0:h:10;
yRK3 = zeros(1,length(t));
yRK3(1) = 1;
yRK2 = yRK3;

f = @(t,y) -y + 2*cos(t);

% Startup with RK3
for n = 1:2
    k1 = f(t(n),yRK3(n));
    k2 = f(t(n) + h/2 , yRK3(n) + (h/2)*k1);
    k3 = f(t(n) + h , yRK3(n) - k1*h + 2*k2*h);
    F = (h/6)*(k1 + 4*k2 + k3);
    yRK3(n+1) = yRK3(n) + F;
end

% Iterating AB3
for n = 3:( length(t) -1 )
    yRK3(n+1) = yRK3(n) + (h/12)*(23*f(t(n),yRK3(n)) - 16*f(t(n-1),yRK3(n-1)) + 5*f(t(n-2),yRK3(n-2)));
end

% Startup with RK2
for n = 1:2
    k1 = f(t(n),yRK2(n));
    k2 = f(t(n) + h , yRK2(n) + h*k1);
    F = (h/2)*(k1 + k2);
    yRK2(n+1) = yRK2(n) + F;
end

% Iterating AB3
for n = 3:( length(t) -1 )
    yRK2(n+1) = yRK2(n) + (h/12)*(23*f(t(n),yRK2(n)) - 16*f(t(n-1),yRK2(n-1)) ...
        + 5*f(t(n-2),yRK2(n-2)));
end

exact = sin(t)+cos(t);
error_RK3 = abs(exact - yRK3);
error_RK2 = abs(exact - yRK2);

% Solutions
figure, hold on
set(gca,'FontSize',12)
plot(t,exact,'k')
plot(t,yRK3,'--r')
plot(t,yRK2,'-.b')
xlabel('t');
ylabel('y');

% Error
figure, hold on
set(gca,'FontSize',12)
plot(t,error_RK3,'--r')
plot(t,error_RK2,'-.b')
xlabel('t');
ylabel('Error');