rungekutta_2.m

来自「1．编写欧拉前差、后差、梯形公式。 2．编写二阶、三阶龙格库塔法通用程序。 」· M 代码 · 共 23 行

function [yt,th]=rungekutta_2(f,y0,t0,td,h)   %our book page 36 function 3_54
y=y0;
t=t0;

yt(:,1)=double(y0);
th(1)=double(t0);

for i=2:fix((td-t0)/h)+1
fn=double(eval(f));                             %solve fn
y=yt(:,i-1)+fn'.*h/2;                           %obtain yn+0.5hfn
t=th(i-1)+h/2;                                  %obtain tn+0.5h
f_=double(eval(f));                             %solve f(tn+0.5h,yn+0.5hfn)

yt(:,i)=yt(:,i-1)+f_'.*h;                       %use yn,f,solve yn+1

th(i)=th(i-1)+h;                                %next time point
y=yt(:,i);                                      %this two is used to solve the next fn
t=th(i);                                        %as above 
end

%n=length(f);                                    %state space function's output
%C=zeros(1,n);C(1)=1;
%out=C*yt;