matlab - MATLAB の非線形微分方程式系

Question

私は MATLAB の初心者で、ルンゲ クッタ アルゴリズムを使用してフリードマン方程式を解決しようとしていました。

ご存じない場合は、フリードマン方程式は次の形式になります。

膨張宇宙のフリードマン方程式

ここで、曲率kは の値で与えられます-1, 0 or 1。

また、圧力の状態方程式は次のようpに与えられます。

圧力の状態方程式

だから、私はルンゲクッタアルゴリズムを持っています：

    function [t,x,y] =rk_2_1(f,g,t0,tf,x0,y0,n)
  h=(tf-t0)/n;
  t=t0:h:tf;
  x=zeros(n+1,1); %reserva memoria para n+1 element(i)os del vect(i)or x(i)
  y=zeros(n+1,1);
  x(1)=x0; y(1)=y0;

  for i=1:n
    k1=h*f(t(i),x(i),y(i));
    l1=h*g(t(i),x(i),y(i));
    k2=h*f(t(i)+h/2,x(i)+k1/2,y(i)+l1/2);
    l2=h*g(t(i)+h/2,x(i)+k1/2,y(i)+l1/2);
    k3=h*f(t(i)+h/2,x(i)+k2/2,y(i)+l2/2);
    l3=h*g(t(i)+h/2,x(i)+k2/2,y(i)+l2/2);
    k4=h*f(t(i)+h,x(i)+k3,y(i)+l3);
    l4=h*g(t(i)+h,x(i)+k3,y(i)+l3);

    x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
    y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;        
  end
end

、しかし、私は3つの問題があります:

1. 関数の密度を記号的に表す方法がわかりません\rho。
2. g=@(t,x,y)時間に関する a の導関数があるため、2 番目の関数 を定義する方法がわかりません。
3. 最後に、 はどうなりk=sqrt(-1)ますか? その結果は象徴的に必要ですが、結果全体が数値的に必要だからです。

私の質問が基本的なもので申し訳ありませんが、これを行う方法がわからず、アドバイスや助けが必要です.

ありがとうございます ：）