matlab代码实现Euler公式求解常微分初值问题

可以使用以下的 Matlab 代码来实现 Euler 公式求解常微分初值问题：

function [t, y] = euler(f, tspan, y0, N)
% f: 常微分方程右端项，tspan: 时间区间，y0: 初值，N: 离散数
t0 = tspan(1);
tN = tspan(2);
t = linspace(t0, tN, N + 1);
h = (tN - t0) ./ N;
y = zeros(size(y0, 1), N + 1);
y(:, 1) = y0;
for n = 1:N
  y(:, n + 1) = y(:, n) + h * f(t(n), y(:, n));
end
end

这里的 `f` 是一个函数句柄，它应该接受一个时间和状态向量并返回右端项。 `tspan` 是一个长度为 2 的数组，包含时间区间的起始点和结束点。 `y0` 是初值，应该是一个列向量。 `N` 是离散数，决定了时间步长。函数的输出是一个时间向量 `t` 和一个状态矩阵 `y`，其中 `t` 是一个长度为 `N+1` 的向量，包含所有时间步的时间点，而 `y` 是一个大小为 `(m, N+1)` 的矩阵，其中 `m` 是状态向量的长度。

调用示例：

% 定义右端项
f = @(t, y) [y(2); -y(1)];
% 时间区间
tspan = [0 10];
% 初值
y0 = [1; 0];
% 离散数
N = 100;
% 求解
[t, y] = euler(f, tspan, y0, N);
% 绘图
plot(t, y(1, :), t, y(2, :));

这里的右端项是二阶常微分方程 $y'' + y = 0$ 的矢量化版本。它返回状态向量的导数，即 $[y_2, -y_1]$，因为 $y_2 = y'$，$y_1 = y$，所以 $y_2' = -y_1$。求解结果将绘制 $y$ 和 $y'$ 在时间上的图像。