1. matlab 新建.m文件，编写龙格-库塔法求解函数

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
n=floor((b-a)/h); %求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值，可以是向量，但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end

2.另外再新建一个.,m文件，定义要求解的常微分方程函数
function dx=fun1(t,x)
dx =zeros(2,1);%初始化列向量
dx(1) =0.08574*x(2)-1.8874*x(1)-20.17;
dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;

3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程

[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20);  %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的龙格库塔函数')
grid

运行步骤3文件即可得到结果，F1为估测值

或者可以调用matlab自带函数ode45求解

方法如下：

[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
grid