你好!代码供参考
function main %主函数
odefun = @(t,y) y^2 + t;
y0 = 1;
tspan = [0, 0.5];
[t,y] = euler(odefun, tspan, y0);
plot(t,y)
end
function [t, y] = euler(odefun, tspan, y0)
% odefun ode函数
% tspan: 求解时间范围
% x0:初值
t = tspan(:);
if(numel(tspan)==2)
t = linspace(tspan(1), tspan(2), 1001)';
end
dy0 = odefun(t(1), y0);
y = zeros(numel(t), numel(dy0));
y(1,:) = y0(:);
for i = 2:1:numel(t)
y(i,:) = y(i-1,:) + odefun(t(i-1),y(i-1,:))*(t(i)-t(i-1));
end
end