function [Y] = RK45(t,X,f,h)K1=f(t,X);K2=f(t+h/2,X+h/2*K1);K3=f(t+h/2,X+h/2*K2);K4=f(t+h,X+h*K3);Y=X+h/6*(K1+2*K2+2*K3+K4);end以上是4階龍格庫塔法的程式碼:自己寫函式，存為f.mfunction dxdt = f (t,x)dxdt(1)=exp(x(1)*sin(t))+x(2);dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f，x(2)是你的gdxdt=dxdt(:);end自己給出t0,x0,h的值（初始時間，初值，步長）如果求t0到t1的軌跡的話:給個例子如下t0=0;t1=5;h=0.02;x0=[-1;-1];T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;for j=1:length(T)-1 X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);endplot(T,X(1,:));hold on;plot(T,X(2,:),"r");具體引數自己設定
function [Y] = RK45(t,X,f,h)K1=f(t,X);K2=f(t+h/2,X+h/2*K1);K3=f(t+h/2,X+h/2*K2);K4=f(t+h,X+h*K3);Y=X+h/6*(K1+2*K2+2*K3+K4);end以上是4階龍格庫塔法的程式碼:自己寫函式，存為f.mfunction dxdt = f (t,x)dxdt(1)=exp(x(1)*sin(t))+x(2);dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f，x(2)是你的gdxdt=dxdt(:);end自己給出t0,x0,h的值（初始時間，初值，步長）如果求t0到t1的軌跡的話:給個例子如下t0=0;t1=5;h=0.02;x0=[-1;-1];T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;for j=1:length(T)-1 X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);endplot(T,X(1,:));hold on;plot(T,X(2,:),"r");具體引數自己設定