非刚性常微分方程的数值解法通常会用四阶龙格库塔算法，其matlab函数对应ode45。
对于dy/dx = f(x,y)，y(0)=y0。
其四阶龙格库塔公式如下：
对于通常计算，四阶已经够用，四阶以上函数f(x,y)计算工作量大大增加而精度提高较慢。
下面以龙格库塔法解洛伦兹方程为例：
matlab代码如下：
main.m：
1 clear all;
2 close all;
3 clc;
4
5 %系统龙格库塔法
6 [t,h] = ode45(@test_fun,[0 40],[12 4 0]);
7 plot3(h(:,1),h(:,2),h(:,3));
8 grid on;
9
10 %自定义龙格库塔法
11 [t1,h1]=runge_kutta(@test_fun,[12 4 0],0.01,0,40);
12 figure;
13 plot3(h1(1,:),h1(2,:),h1(3,:),'r')
14 grid on;
runge_kutta.m（函数参考网络）：
1 %参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
2 function [x,y]=runge_kutta(ufunc,y0,h,a,b)
3 n=floor((b-a)/h); %步数
4 x(1)=a; %时间起点
5 y(:,1)=y0; %赋初值，可以是向量，但是要注意维数
6 for i=1:n %龙格库塔方法进行数值求解
7 x(i+1)=x(i)+h;
8 k1=ufunc(x(i),y(:,i));
9 k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2);
10 k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2);
11 k4=ufunc(x(i)+h,y(:,i)+h*k3);
12 y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;
13 end
test_fun（洛伦兹方程）：
1 %构造微分方程
2 function dy=test_fun(t,y)
3 a = 16;
4 b = 4;
5 c = 45;
6
7 dy=[a*(y(2)-y(1));
8 c*y(1)-y(1)*y(3)-y(2);
9 y(1)*y(2)-b*y(3)];
得到很经典的洛伦兹吸引子，结果如下：
