function y = trapez(h, f, t0, y0, T, A)
%TRAPEZ 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,
        fun = @(x) x - y(:,n) - 0.5 * h * f(t(n),y(:,n)) - 0.5 * h * f(t(n+1),x);
        x0 = y(:,n) + h * f( t(n), y(:,n));
        y(:, n+1) = fsolve(fun,x0,options);
    end
else
    for n = 1:N,
        x = (eye(d) + 0.5 * h * A(t(n))) * y(:,n);
        y(:, n+1) = (eye(d) - 0.5 * h * A(t(n+1))) \ x;
    end  
end

end

