% This script illustrates Euler's method 
% for a non-Lipschitz problem 
% y' = ln(3)(y-floor(y)-3/2),   y(0)=0
% on the interval [0,8]
% with exact solution
% y(t) = -n +1/2(1-3^(t-n)),  n <= t <= n+1, n=0,1,...
clear all
clc
close all
scrnsz = get(0, 'Screensize');
f = @(t,y) log(3)*(y-floor(y)-3/2);
t0=0; y0=0;    % initial condition
tmax = 8;
Nvalues = [80 800];
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,8,201);
    nn=floor(tt);
    plot(tt,-nn+1/2*(1-3.^(tt-nn)),'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)+floor(t(2:end))-1/2*(1-3.^(t(2:end)-floor(t(2:end))))));
    plot(t(2:end),err,'LineWidth',2)
    maxerr = max(abs(y+floor(t)-1/2*(1-3.^(t-floor(t)))));
    if N>Nvalues(1)
        orderestimate = log(oldmax/maxerr)/log(oldh/h)
    end
    oldmax = maxerr;
    oldh = h;
    if N<Nvalues(end)
        pause
    end
end
    