% euler_backward Euler approximation for ODE initial value problem
%               Backward Euler method
% [t, u] = euler_backward(f, t0, t1, u0, n) returns a column vector of time
% points and a column vector of function values, 
% f = f(t, u) is the RHS, 
% t0 is the initial time, 
% t1 is the final time, 
% u0 is the initial function value, 
% n is the number of steps of the numerical method
%
% Lei Zhang, 18-09-2012
function [t,y]=euler_backward(f, fy, t0, t1, y0, numSteps)

if nargin == 0
    test_euler_backward();
    return
end

% calculate time step dt
dt = (t1 - t0)/numSteps; 

% initialize t and y as column vectors
t = linspace(t0, t1, numSteps+1);
y = zeros(size(t)); y(1) = y0;

for k = 1 : numSteps
  Y = y(k) + dt * f(t(k), y(k));
  
  % Newton iteration
  [Y, isConverged]= newton4euler(f, fy, t(k+1), y(k), Y, dt);

  if ~ isConverged
    error(['back_euler failed to converge at step ', num2str(k)])
  end

  y(k+1) = Y;
end
end

function test_euler_backward

choice = 1;
switch choice
    case 1
        f = @(t, u)pi*cos(pi*t);
        fy = @(t, u)0;
    case 2
        f = @(t, u)(t + u);
        fy = @(t, u)1;
    case 3
        f = @(t, u)exp(t);
        fy = @(t, u)0;
    case 4
        f = @(t, u)-20*u + 20*t + 21;
        fy = @(t, u)-20;
    case 5
        f = @(t, u)u.^2.*(t-t.^3);
        fy = @(t, u)u.*(t-t.^3);
    case 6
        f = @(t, u)t^2 - u;
        fy = @(t, u)-1;
    case 7
        f = @(t, u)u.^(0.2);
        fy = @(t, u)0.2*u.^(-0.8);
end
t0 = -1/2;
t1 = 3;
u0 = 1;

n = 2^16;
[t_ex, u_ex] = euler_backward(f, fy, t0, t1, u0, n);
figure
plot(t_ex, u_ex)
hold

p = 3:10;
e = zeros(size(p));
for k = 1:length(p)
    n(k) = 2^p(k); 
    [t, u] = euler_backward(f, fy, t0, t1, u0, n(k));
    if k == 3
        plot(t, u, 'rx-')
    end
    e(k) = abs(u(end) - u_ex(end));
end
figure
loglog(n, e); 
P = polyfit(log(e), log(n), 1) 
hold on
loglog(n, 1./n, 'Color','r');
xlabel('n'); ylabel('error');
text(n(end), 1/n(end), 'n^{-1}');
end