oderk

function [t,xt]=oderk(f,tI,x0,n)
% 龙格库塔 求常微分方程初值问题数值解法， 本次实验我写的是四阶龙格库塔
% f 导数值 tI x的定义域 x0 初值
if nargin<4
    n=1;
end
% 定义域区间为2时，将他分为100个区间
if length(tI)==2
    t=linspace(tI(1),tI(2));
else
    t=tI;
end
switch n
    case 1
% 四阶龙格库塔公式
        k1=@(t,x,h)f(t,x);
        k2=@(t,x,h)f(t+h/2,x+h/2*k1(t,x,h));
        k3=@(t,x,h)f(t+h/2,x+h/2*k2(t,x,h));
        k4=@(t,x,h)f(t+h,x+h*k3(t,x,h));
        w=@(t,x,h)(k1(t,x,h)+2*k2(t,x,h)+2*k3(t,x,h)+k4(t,x,h))/6;
end
%非等距分割，求导
dt=diff(t);
%定义数组，保存xt
xt=zeros(length(x0),length(t));
%初值，保存零点值
xt(:,1)=x0;
%从t-t-1 开始迭代
for k=1:length(t)-1
    xt(:,k+1)=xt(:,k)+dt(k)*w(t(k),xt(:,k),dt(k));
end
end