Как решить ODE одновременно, используя метод средней точки?

Question

Справочная информация Вопрос

Это мой файл функций средней точки, в котором я тоже не уверен. Я отредактировал оригинальный файл функции средней точки, чтобы удовлетворить этот вопрос, но я не знаю, правильно ли он или нет.

function [t,a] = midpoint2(dadt,tspan,z0,h)
% [t,y] = midpoint2(dydt,tspan,y0,h)

% uses midpoint method to solve an ODE
%
% INPUTS:
%  - dadt = function handle of the ODE, f(t,a)
%  - tspan = [<initial value>, <final value>] of independent variable
%  - z0 = initial value of dependent variable
%  - h = step size
% OUTPUTS:
%  - t = vector of time
%  - a = vector of solution for degree

% Input Validation: tspan
if ~(tspan(2)>tspan(1))
    error('upper limit must be greater than lower')
end

% define t vector
t = (tspan(1):h:tspan(2))';
n = length(t);
% if necessary, add an additional t so that range goes up to tspan(2)

% Preallocate a to improve efficiency
a = z0*ones(n,1);
% Implement Midpoint method
for i = 1:n-1
    t_half = a(i)+h/2*dadt(t(i),a(i));
    a(i+1) = a(i)+h*(dadt(t(i)+h/2,t_half));
end