非剛性常微分方程的數值解法通常會用四階龍格庫塔算法，其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)];
得到很經典的洛倫茲吸引子，結果如下：
