% This script illustrates the use of Euler's method
% for the problem 
% y' = -y + 2exp(-t)cos(2t),   y(0)=0
% on the interval [0,10]
% with exact solution
% y(t) = exp(-t)sin(2t)
clear all
clc
close all
scrnsz = get(0, 'Screensize');
f = @(t,y) -y+2*exp(-t).*cos(2*t);
t0=0; y0=0;    % initial condition
tmax = 10;
Nvalues = [20 100 500];
for N=Nvalues
    h=tmax/N;
    [t,y] = Euler(t0,y0,f,h,N);
    figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
    hold on 
    title(sprintf('Approximate Euler solution with %d points',N)) 
    xlabel('t')
    ylabel('y')
    plot(t,y,'b','LineWidth',2); 
    plot(t,y,'ro','LineWidth',2); 
    legend('Euler soluton','Euler points',-1);
    pause
    title(sprintf('Approximate Euler and exact solutions with %d points',N)) 
    tt=linspace(0,10,201);
    plot(tt,exp(-tt).*sin(2*tt),'g','LineWidth',2);  % exact solution
    legend('Euler soluton','Euler points', 'Exact Solution',-1);
    hold off
    pause
    figure(length(Nvalues)+1)
    hold on 
    title('Logs of Errors') 
    xlabel('t')
    ylabel('log(|Euler-exact|)')
    err = log(abs(y(2:end)-exp(-t(2:end)).*sin(2*t(2:end))));
    plot(t(2:end),err,'LineWidth',2)
    maxerr = max(abs(y-exp(-t).*sin(2*t)));
    if N>Nvalues(1)
        orderestimate = log(oldmax/maxerr)/log(oldh/h)
    end
    oldmax = maxerr;
    oldh = h;
    if N<Nvalues(end)
        pause
    end
end
    