Здравствуйте.
Прошу помощи в разборе задачи: нужно решить ОДУ 2 степени в матлаб, вида

, только вручную, модифицированным (с пересчетом) методом Эйлера.
Для начала нужно сделать обыкновенный гармонический осциллятор, чтобы убедиться, что все понимаю правильно. Самое главное- вынести правую часть в отдельный m-файл. И то ли лыжи не едут, то ли как обычно все, но гармонический осциллятор не получается.
Уравнение использую вида

Правую часть пытаюсь сделать сразу, на будущее, зависящей от

и

.
Вот код:
файл euler_s.m
clc;
clear all;
sec = 1;
h = 0.2; %шаг интегрирования
zp = zeros(n, sec/h+1); %predictor
zc = zeros(n, sec/h+1); %corrector 2st func.
yp = zeros(n, sec/h+1); %predictor
yc = ones(n, sec/h+1); %corrector 1st func.
for t = 2:sec/h+1
%prediction
yp(t) = yc(t-1) + h*zc(t-1);
zp(t) = zc(t-1) + h*Fun(zc(t-1),t-1);
%correction
yc(t) = yc(t-1)+h*(zc(t-1) + zp(t))/2
zc(t) = zc(t-1)+h*(Fun(zc(t-1),t-1) + Fun(zp(t),t))/2;
end
plot([1:sec/h+1],yc),grid;
Здесь уравнение 2 степени с помощью известной замены становится системой ОДУ 1 степени:

Далее файл правой части Fun.m:
function = Fun(r,t);
Fun = -r;
Вот правая часть как раз для гармонического осциллятора. График показывает горизонтальную прямую на высоте начального условия. Не понимаю ошибки у себя, прошу помощи.
Спасибо.