function [Y,T] = expEuler(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;

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