function [Y,T] = expMPR(f,tspan,yinit,h)

t0     = tspan(1);
tfinal = tspan(2);

N = ceil((tfinal-t0)/h);

Y = zeros(size(yinit,1),N+1);
Y(:,1) = yinit;

T = linspace(t0,tfinal,N+1);
h = (tfinal-t0)/N;

% Euler-step
Y(:,2) = Y(:,1) + h*f(T(1),Y(:,1));

for k = 2:N
  Y(:,k+1) = Y(:,k-1) + 2*h*f(T(k),Y(:,k));
end