function [t,y] = simpson(h, f, t0, y0, T, A)

t = t0:h:T;
N = length(t)-1;
d = size(y0,1);
y = zeros(d, N+1);
y(:,1) = y0;
[dummy,Yh]=ode45(f,t(1:2),y0);
Y(:,2)=Yh(end,:)';


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

end