function [tvec,yvec] = EulerSystem(t0,y0,f,h,N)
% [tvec,yvec] = EulerSystem(t0,y0,y1,f,h,N)
% Euler 1st-order method
% Inputs
% t0,y0: initial condition
% f:     right-hand side function f(t,y)
% h:     stepsize
% N:     number of steps
% Outputs
% tvec: vector of t values
% yvec: vector of corresponding y values

d = length(y0);
yvec = zeros(N+1,d);
tvec = linspace(t0,t0+N*h,N+1)';
yvec(1,:) = y0';
for n=1:N
   fvalue1 = feval(f,tvec(n),yvec(n,:));
   yvec(n+1,:) = yvec(n,:)+h*fvalue1';
end