function [t,y] = backeuler_sys(f, lh, ic, stepsize)
%Usage: [t,y] = backeulers_sys(f, lh, ic, stepsize)
%Backward/Implicit Euler's Method implementation for ODE system y' = f(t,y) 
%hence y can be a vector. Here lh is the interal lh[1] <= t <= lh[2] over 
%which the solution is computed and y(lh[1]) = ic, where ic is the vector
%of initial conditions at time t = lh[1]. The step size is in argument 
%"stepsize". If the number of steps from initial to final time is not an 
%integer, the last step is incremental of proper length. Solution times 
%are in t, solution values in y.
%THE FUNCTION f SHOULD RETURN A ROW VECTOR.

%Prep for marching out in time
finalT = lh(2);
N = floor((finalT-lh(1))/stepsize); %Round down
T1 = lh(1)+N*stepsize; %Corresponding time
Nf = N+1; %Number of steps we'll take, unless T1 = finalT
if abs(T1-finalT) < 1.0e-8 %To some tolerance
    Nf = N;
end

%Make that ic is a row vector.
ic = ic(:)';

%Set up arrays for solution return
N = length(ic);
t = zeros(Nf+1,1);
y = zeros(Nf+1,N);
t(1) = lh(1); %Initial time
y(1,:) = ic; %Initial value

%March out in time
for i=2:Nf
    t(i) = t(i-1) + stepsize; %Update time
    
    %Extrapolate solution. The equation to solve is
    %ynew = yold + h*f(told,ynew). Use fsolve for solving.
    %Define f2 as appropriate function to zero.
    f2 = @(w) y(i-1,:) + stepsize*f(t(i-1),w) - w;
    %evalc used to suppress lots of annoying output from fsolve
    [~,y(i,:)] = evalc('fsolve(f2,y(i-1,:))');
end

%Final, possibly incremental step
t(Nf+1) = finalT; %Update time
stepsize = finalT-t(Nf);
    
%Extrapolate solution. The equation to solve is
%ynew = yold + h*f(told,ynew).
f2 = @(w) y(Nf,:) + stepsize*f(t(Nf),w) - w;
[~,y(Nf+1,:)] = evalc('fsolve(f2,y(Nf,:))');