MATLAB里ode45命令用4阶龙格-库塔解一阶微分方程（常系数/变系数）都能解。
首先把微分方程化为dx/dt = f(x)的标准形式写成函数放在一个文件里，再写主程序用ode45调用这个文件。举例如下：
函数文件：方程的定义（一定要用函数名命名）
% di/dt = L-1u - L-1(WL+R)ifunction iprime = hw11_2(t,i)L = [.0072, 0, 0, .1, .0054, 0;0, .007, 0, 0, 0, .0026;0, 0, .001, 0, 0, 0;1.5*.1, 0, 0, 2.5, .125, 0;1.5*.0054, 0, 0, .125, .0068, 0;0, 1.5*.0026, 0, 0, 0, .0016];W = zeros(6,6);W(1,2) = -100*pi;W(2,1) = 100*pi;R = diag([.002,.002,.002,.4,.015,.015]);u = zeros(6,1);u(4) = 400;iprime = L\u - L\(W*L+R)*i;
主程序文件：调用ode45
% solve function hw11_2theta0 = 0;h = 1e-3; % stepi0 = zeros(1,6);i0(4) = 400/.4;ts = 0:h:1; % time span[t,i] = ode45('hw11_2',ts,i0); % solvethetaA = 100*pi*ts + theta0;thetaB = thetaA - pi/3*2;thetaC = thetaA + pi/3*2;n = length(ts);ia = zeros(n,1);ib = zeros(n,1);ic = zeros(n,1);% inverse Park transformfor j = 1:ntt = ts(j);thetaA = 100*pi*tt + theta0;thetaB = thetaA - 2*pi/3;thetaC = thetaA + 2*pi/3;ia(j) = cos(thetaA)*i(j,1) - sin(thetaA)*i(j,2) + i(j,3);ib(j) = cos(thetaB)*i(j,1) - sin(thetaB)*i(j,2) + i(j,3);ic(j) = cos(thetaC)*i(j,1) - sin(thetaC)*i(j,2) + i(j,3);endfigure(1);plot(ts,ia,'y',ts,ib,'g',ts,ic,'r');title(['three phase short circuit (\theta_0 = ',num2str(theta0),')']);xlabel('t/s');ylabel('i/A');legend('I_a','I_b','I_c');absia = abs(ia);[iamax,jamax] = max(absia);tamax = ts(jamax);
