Matlab龙格库塔法求解非线性常微分方程

请教一下
这个怎么用Matlab编程求解这个非线性常微分方程
用龙格库塔法

调用这个方法

function [T,X,dX] = ODE_RK4( Hfun,t,h,x0 )
    T = t(1):h:t(2);
    % 计算
    N = length(T);
    x0 = x0(:);  
    x0 = x0';                  % 初值变为行向量  
    m = length(x0);            % 状态量维数
    X = zeros(N,m);            % 初始化状态量
    dX = zeros(N,m);           % 状态导数
    X(1,:) = x0;
    for k = 2:N
        h = T(k) - T(k-1);
        K1 = Hfun( T(k-1)     , X(k-1,:)' );    
        K2 = Hfun( T(k-1)+h/2 , X(k-1,:)'+h*K1/2 ); 
        K3 = Hfun( T(k-1)+h/2 , X(k-1,:)'+h*K2/2 ); 
        K4 = Hfun( T(k-1)+h   , X(k-1,:)'+h*K3 ); 
        X(k,:) = X(k-1,:)' + (h/6) * ( K1 + 2*K2 + 2*K3 + K4 );      
        dX(k-1,:) = (1/6) * ( K1 + 2*K2 + 2*K3 + K4 );
    end
    dX(N,:) = Hfun( T(N),X(N,:) );
end