function y = euler_implicit(h, f, t0, y0, T, A)
%EULER_EXPLICIT Summary of this function goes here
%   Detailed explanation goes here

t = t0:h:T;
N = length(t)-1;
d = size(y0,1);
y = zeros(d, N+1);
y(:,1) = y0;

if nargin < 6,
    options = optimset('Display','off');
    for n = 1:N,
        x0 = y(:,n) + h * f( t(n), y(:,n));
        fun = @(x) x - y(:,n) - h * f(t(n+1),x);
        y(:, n+1) = fsolve(fun,x0,options);
    end
else
    for n = 1:N,
        y(:, n+1) = (eye(d) - h * A(t(n+1))) \ y(:,n);
    end  
end

end

