% Code for lecture 16 final example
clear, close all;

%setup
h = 0.1; y0 = 1; T = 8;
N = T/h; t = 0:h:T;
%pre-allocate arrays
xx = zeros(size(t));
y_ee = xx; y_ie = xx; y_mp = xx;
ysolRK = xx; ymult = xx; ymultun = xx;

%assign initial value
y_ee(1) = y0; y_ie(1) = y0; y_mp(1) = y0; ysolRK(1) = y0;

%reference solution with ode45
opts = odeset('AbsTol', 1e-10, 'RelTol', 1e-10);
g = @(t) sin(t.^2); gg = @(t,y) sin(t.^2).*y; % change and play
%g = @(t) cos(t.^2); gg = @(t,y) cos(t.^2).*y;
%g = @(t) 1./(1+t.^2); gg = @(t,y) 1./(t.^2+1).*y;

[tref,yref] = ode45( @(t,y) gg(t,y), [0, T], 1, opts);
yode45 = ode45( @(t,y) gg(t,y), [0, T], 1, opts);

yex = deval(yode45,t);

ymult(1:4) = yex(1:4); ymultun(1:3) = yex(1:3); % for multistep

%numerical solutions
for ii = 1:N
    %explicit Euler
    y_ee(ii+1) = (1 + h*g(t(ii)))*y_ee(ii);
    %implicit Euler
    y_ie(ii+1) = y_ie(ii)/(1 - h*g(t(ii+1)));
    %implicit Midpoint Rule
    y_mp(ii+1) = (1 + h*g(t(ii)+h/2)/2)*y_mp(ii)/(1 - h*g(t(ii)+h/2)/2);
    
    %Runge-Kutta 4
    k1 = gg(t(ii),ysolRK(ii));
    k2 = gg(t(ii)+h/2,ysolRK(ii) + h*k1/2);
    k3 = gg(t(ii)+h/2,ysolRK(ii) + h*k2/2);
    k4 = gg(t(ii)+h,ysolRK(ii) + h*k3);
    ysolRK(ii+1) = ysolRK(ii) + h/6*(k1 + 2*k2 + 2*k3 + k4);
    
    %{
     % RK2
        k1 = gg(t(ii),ysolRK(ii));
        k2 = gg(t(ii)+h/2,ysolRK(ii) + h*k1/2);
        k3 = gg(t(ii)+h,ysolRK(ii) + h*k2);
        ysolRK(ii+1) = ysolRK(ii) + h/3*(k1 + k2 + k3);
    %}
end
% multistep (stable and unstable) and plot

% Adams Bashforth
for ii = 1:N-3
    ymult(ii+4) = ymult(ii+3) + h/24*(55*gg(t(ii+3),ymult(ii+3))...
        -59*gg(t(ii+2),ymult(ii+2))...
        +37*gg(t(ii+1),ymult(ii+1))...
        -9*gg(t(ii),ymult(ii)));
end
% unstable multistep
for ii = 1:N-1
    ymultun(ii+2) = -4*ymultun(ii+1) +5*ymultun(ii)...
        + h*(4*gg(t(ii+1),ymultun(ii+1))-2*gg(t(ii),ymultun(ii)));
end
%{
% another unstable multistep
for ii = 1:N-1
ymultun(ii+2) = 3*ymultun(ii+1) - 2*ymultun(ii)...
+ h*((13/12)*gg(t(ii+2),ymultun(ii+2))...
    -(5/3)*gg(t(ii+1),ymultun(ii+1))...
    -(5/12)*gg(t(ii),ymultun(ii)));
end
%}

% plot solutions
MS = 'Markersize'; LW = 'linewidth'; FS = 'fontsize'; CO = 'Color'; TEX = 'interpreter'; tex = 'latex';
lw = 2; ms = 12; fs = 14;
plot(t,y_ee,'.-',t,y_ie,'.-',t,y_mp,'.-',tref, yref,'k--')
xlabel('x',FS,fs),ylabel('y',FS,fs)
h_legend = legend('Explicit Euler','Implicit Euler','Midpoint', 'Exact',FS,fs)
set(h_legend,FS,15,'Location','Best');
grid on

%export_figmy(['lec16_soln']) % YN's code for saving figure (ignore)

%% plot error
figure
plot(t,yex-y_ee,'.--',t,yex-y_ie,'.--',t,yex-y_mp,'.--',t,yex-ysolRK,'.--',t,yex-ymult,'--')
h_legend = legend('Explicit Euler','Implicit Euler','Midpoint', 'RK4','Adams-Bash',FS,fs)
set(h_legend,FS,15,'Location','Best');
xlabel('x',FS,fs),ylabel('error',FS,fs)
shg, grid on

%export_figmy(['lec16_err'])

%% unstable method
figure,
semilogy(t,abs(yex-ymultun),'.-')
title('error with unstable method!')
xlabel('x',FS,fs),ylabel('error',FS,fs)
grid on

%alignfigs, export_figmy(['lec16_errunstable'])

%% order of convergence

y0 = 1; T = 8;
Ns = 2.^(3:10);
ip = 1;erree = []; errie = []; errmp = []; errRK = []; errmult = [];

for N = Ns
    ip = ip+1;
    t = linspace(0,T,N+1);
    h = t(2)-t(1);
    
    yex = deval(yode45,t); % exact vals
    
    %pre-allocate arrays
    xx = zeros(size(t));y_ee = xx; y_ie = xx; y_mp = xx;
    ysolRK = xx; ymult = xx;
    %assign initial value
    y_ee(1) = y0; y_ie(1) = y0; y_mp(1) = y0; ysolRK(1) = y0;
    ymult(1:4) = yex(1:4);
    
    for ii = 1:N
        %explicit Euler
        y_ee(ii+1) = (1 + h*g(t(ii)))*y_ee(ii);
        %implicit Euler
        y_ie(ii+1) = y_ie(ii)/(1 - h*g(t(ii+1)));
        %implicit Midpoint Rule
        y_mp(ii+1) = (1 + h*g(t(ii)+h/2)/2)*y_mp(ii)/(1 - h*g(t(ii)+h/2)/2);
        
        k1 = gg(t(ii),ysolRK(ii));
        k2 = gg(t(ii)+h/2,ysolRK(ii) + h*k1/2);
        k3 = gg(t(ii)+h/2,ysolRK(ii) + h*k2/2);
        k4 = gg(t(ii)+h,ysolRK(ii) + h*k3);
        ysolRK(ii+1) = ysolRK(ii) + h/6*(k1 + 2*k2 + 2*k3 + k4);
        
        %{
        k1 = gg(t(ii),ysolRK(ii));
        k2 = gg(t(ii)+h/2,ysolRK(ii) + h*k1/2);
        k3 = gg(t(ii)+h,ysolRK(ii) + h*k2);
        ysolRK(ii+1) = ysolRK(ii) + h/3*(k1 + k2 + k3);
        %}
    end
    erree = [erree norm(yex-y_ee,inf)];
    errie = [errie norm(yex-y_ie,inf)];
    errmp = [errmp norm(yex-y_mp,inf)];
    errRK = [errRK norm(yex-ysolRK,inf)];
    
    
    % Adams Bashforth
    for ii = 1:N-3
        ymult(ii+4) = ymult(ii+3) + h/24*(55*gg(t(ii+3),ymult(ii+3))...
            -59*gg(t(ii+2),ymult(ii+2))...
            +37*gg(t(ii+1),ymult(ii+1))...
            -9*gg(t(ii),ymult(ii)));
    end
    errmult = [errmult norm(yex-ymult,inf)];
end

clf
NNs = 1./Ns;
loglog(NNs,erree,'.-',NNs,errie,'.-',NNs,errmp,'.-',NNs,errRK,'.-',NNs,errmult,'.-','markersize',12)
xlabel('h',FS,fs),ylabel('error',FS,fs)
grid on, shg
h_legend = legend('e-Euler','i-Euler','midpoint', 'R-K4','Adam-B4',FS,fs);
set(h_legend,FS,15,'Location','Best');
set(gca,FS,fs)
hold on
loglog(NNs,2*NNs/NNs(1),'k--',NNs,NNs.^2/(NNs(1)^2),'k--',NNs,NNs.^4/(NNs(1)^4),'k--','HandleVisibility','off')
text(NNs(end-2),NNs(end-2)/NNs(1),'O(h)',FS,fs)
text(NNs(end-2),(NNs(end-2)/NNs(1))^2,'O(h^2)',FS,fs)
text(NNs(end-2),(NNs(end-2)/NNs(1))^4,'O(h^4)',FS,fs)
yy = ylim; ylim([min(errRK) yy(2)])

% export_figmy(['lec16_convorders'])