非剛性常微分方程的數值解法通常會用四階龍格庫塔演算法，其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)];

得到很經典的洛倫茲吸引子，結果如下：