clear; 
clc;

%系统龙格库塔法
[t,h] = ode45(@test_fun,[0 40],[12 4 0]);
plot3(h(:,1),h(:,2),h(:,3));
grid on;

%自定义龙格库塔法
[t1,h1]=runge_kutta(@test_fun,[12 4 0],0.01,0,40);
figure;
plot3(h1(1,:),h1(2,:),h1(3,:),\'r\')
grid on;

%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
function [x,y]=runge_kutta(ufunc,y0,h,a,b)
n=floor((b-a)/h);       %步数
x(1)=a;                 %时间起点
y(:,1)=y0;              %赋初值，可以是向量，但是要注意维数
for i=1:n               %龙格库塔方法进行数值求解    
    x(i+1)=x(i)+h;    
    k1=ufunc(x(i),y(:,i));  
    k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2);    
    k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2);   
    k4=ufunc(x(i)+h,y(:,i)+h*k3);   
    y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;
end

%构造微分方程
function dy=test_fun(t,y)
a = 16;
b = 4;
c = 45;

dy=[a*(y(2)-y(1));
    c*y(1)-y(1)*y(3)-y(2);
    y(1)*y(2)-b*y(3)];

 clear; 
clc;
%自定义龙格库塔法
%P = [0.1, 0.25, 0.25, 0.25, 0.25];
P = [1.9, 0.8, 0.8, 0.8, 0.8];
[t ,h ]=runge_kutta(@test_fun,P,0.1,0,72); 
plot(t,h(1,:));
hold on;
plot(t,h(1,:));
plot(t,h(2,:));
plot(t,h(3,:));
plot(t,h(4,:));
plot(t,h(5,:));
hold off;

%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
function [x,y]=runge_kutta(ufunc,y0,h,a,b)
n=floor((b-a)/h);       %步数 
%x = zeros(1,n)           %时间起点
x(1) = a;
y(:,1)=y0;              %赋初值，可以是向量，但是要注意维数
for i=1:n               %龙格库塔方法进行数值求解
    x(i+1) = x(i)+h;     
    k1 = ufunc(x(i),y(:,i));  
    k2 = ufunc(x(i)+h/2,y(:,i)+h*k1/2);    
    k3 = ufunc(x(i)+h/2,y(:,i)+h*k2/2);   
    k4 = ufunc(x(i)+h,y(:,i)+h*k3);   
    y(:,i+1) = y(:,i)+h*(k1+2*k2+2*k3+k4)/6;
end

%构造微分方程
function dy=test_fun(t,P)  
vs = 0.76, vm = 0.65, vd = 0.95;
k1 = 1.9, k2 = 1.3 ,ks = 0.38;
M = P(1), P0 = P(2), P1 = P(3), P2 = P(4), PN = P(5);
V1 =3.2 , V2 = 1.58, V3 = 5 , V4 = 2.5;
K1 = 1, K2 = 2, K3 = 2, K4 = 2, Kd = 0.2, Km = 0.5;
n = ;
dy= [vs*(K1^n/(K1^n+PN^n))-vm*M/(Km+M);
    ks*M-V1*P0/(K1+P0)+V2*P1/(K2+P1);
    V1*P0/(K1+P0)-V2*P1/(K2+P1)-V3*P1/(K3+P1)+V4*P2/(K4+P2);
    V3*P1/(K3+P1)-V4*P2/(K4+P2)-k1*P2+k2*PN-vd*P2/(Kd+P2);
    k1*P2-k2*PN];