%偏微分方程組N=4的龍格－庫塔方法求解
 %形式：求t在區間[a,b]內微分方程組：
 %     x'1(t)=f1(t,x1(t),...xn(t))
 %     x'2(t)=f2(t,x1(t),...xn(t))
 %     ...
 %     x'n(t)=fn(t,x1(t),...xn(t))
 % 初值條件：x1(a)=a1,x2(a)=a2...xn(a)=an,的近似解
 %_______________________________________________________
 %調用格式：先定義F函數,既定義方法如下(以三個方程組成的方程組為列：
 %function  Z=F(X)
 %          x=X(1);
 %          y=X(2);
 %          z=X(3);
 %          Z=zeros(1,2);  
 %          Z(1)=3*x+2*y+z; //Z(1)相當于表達式f1(t,x1(t)...)
 %          Z(2)=x+2*y-z;   //Z(2)相當于表達式f2(t,x1(t)...)
 %          Z(3)=x-3*z;     //Z(3)相當于表達式f3(t,x1(t)...)
 %然后直接調用[T,Z]=rks4(’F‘,a,b,Za,M)，其中F上面的函數
 %                  a合b是t的范圍，Za為初值向量【X1(a) X2(a)...】
 %                  M為需要把區間分為多少次進行代換即(b-a)/M
 %                  Z的一列為x1(t),第二列為x2(t)........
 %_______________________________________________________
 
 function [T,Z]=rks4(F,a,b,Za,M)
 h=(b-a)/M;
 T=zeros(1,M+2);
 Z=zeros(M+1,length(Za));
 T=a:h:b;
 Z(1,:)=Za;
 
 for j=1:M
     k1=h*feval(F,[T(j),Z(j,:)]);
     k2=h*feval(F,[T(j)+h/2,Z(j,:)+k1/2]);
     k3=h*feval(F,[T(j)+h/2,Z(j,:)+k2/2]);
     k4=h*feval(F,[T(j)+h,Z(j,:)+k3]);
     Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
 end