The method is ok, though could be more streamlined, for example:
f1=@(x) 5*x+50;
f2 =@(x,y) x*10 +10*y;
x1=0;
y1=0;
x2=0;
y2= 0;
xn=10;
h=0.2;
n = round(xn/h);
for i = 1:n
y1=y1+h*f1(x1);
x1=x1+h;
y2=y2+h*f2(x2,y2);
x2=x2+h;
end
fprintf('\n x1 y1 x2 y2 ');
fprintf('\n%4.3f %4.3f %4.3f %4.3f ',x1,y1,x2,y2);
Your functions (especially f2) show up the deficiencies in the Euler method (the true final value at x = 10, for f1 is 750, and for f2 is ~2.688*10^42).
