急!!求使用matlab求解三個一階ODE聯立方程

問題:用matlab求解三個一階ODE聯立方程(初始值問題)以下為matlab 兩個副程式及主程式：function rs= rk4syskzpi6(X0

a

b

m)%to return the approximation values of X(t)%by RK4 method

X0 is the initial%t is in[a

b] n=length(X0); rs=zeros(n

m) ; K1=zeros(n

1) ;K2=zeros(n

1) ; K3=zeros(n

1) ;K4=zeros(n

1) ; t=a ; X=X0 ;X2=zeros(n

1);X3=X2; X4=X2 ; h=(b-a)/m ; for k=1:m%compute K1

K2

K3and K4K1=fxsyskzpi6(t

X); for i=1:nX2(i)=X(i) 1/2*h*K1(i); endK2=fxsyskzpi6(t h/2

X2); for i=1:nX3(i)=X(i) 1/2*h*K2(i); endK3=fxsyskzpi6(t h/2

X3); for i=1:nX4(i)=X(i) h*K3(i); endK4=fxsyskzpi6(t h

X4);%compute X(t h) for i=1:nX(i)=X(i) h/6*(K1(i) 2*K2(i) 2*K3(i) K4(i)); rs(i

k)=X(i) ; end %for i t=t h ; end % for k%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function fX=fxsyskzpi6(t

X)%fX(1)=dI/dy、X(1)=I%fX(2)=d(eta)/dy、X(2)=eta%fX(3)=d(C/C*)/dy、X(3)=C/C*fX=zeros(length(X)

1);fX(1)=(1/(9.5*10^3))*X(3)*exp(X(2)); fX(2)=2*X(1); fX(3)=-0.3737*(1-2*X(1));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%a=0; b= 1; delt=0.01;k=(b-a)/delt;m=100;X0=[0;8.9432; 1];t=a:delt:b;Rs =rk4syskzpi6(X0

a

b

m);%%為求解di/dx之值%C=Rs(1

:);D = zeros(1

(k-1));%定義 D是一個1

k-1的零矩陣;for n=1:(k-1);D(n)=[(C(1

n 1)-C(1

n))/delt];%di/dx值end%E=D