function
R = EulerM
(
a, b, M, y_0
)
h =
(
b - a
)
/ M;
Y =
zeros
(
1
, M+
1
)
;
Y
(
1
)
= y_0;
T = a:h:b;
for
i
=
1
:M
Y
(
i
+
1
)
= Y
(
i
)
+ h * f
(
T
(
i
)
, Y
(
i
)
)
;
Y
(
i
+
1
)
= Y
(
i
)
+ h /
2
*
(
f
(
T
(
i
)
, Y
(
i
)
)
+ f
(
T
(
i
+
1
)
, Y
(
i
+
1
)
)
)
;
end
R = T' Y';