Matlab欧拉法求解非线性常微分方程

请教一下
这个怎么用Matlab编程求解这个非线性常微分方程
包括显式欧拉法，隐式欧拉法，两步欧拉法，改进欧拉法
包括图像

clear
clc
close all 
%%  仿真步长 h=0.01 时
Hfun = @(t,x) [ x(1)-x(2)-x(1)*x(3);x(1)+x(2)-x(2)*x(3); x(1)*x(1)+x(2)*x(2)-x(3)];   % 一阶微分方程导数表达式 
% 参数
t = [0 12];     % 时间范围
h = 0.01;       % 时间步长
x0 = [2 0 1];   % 初始状态 
% 显式欧拉法求解
[T1,X1] = ODE_ExplicitEuler( Hfun,t,h,x0 );
% 隐式欧拉法求解
[T2,X2] = ODE_ImplicitEuler( Hfun,t,h,x0 );
% 两步欧拉法求解
[T3,X3] = ODE_2OrderEuler( Hfun,t,h,x0 );
% 改进欧拉法求解
[T4,X4] = ODE_ImprovedEuler( Hfun,t,h,x0 );
% Matlab自带ode45求解
[T5,X5] = ode45( Hfun,t,x0 ); 
% 绘图对比
figure
subplot(311)
plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1),T4,X4(:,1),T5,X5(:,1))
xlabel('Time(s)')
ylabel('x_1')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
subplot(312)
plot(T1,X1(:,2),T2,X2(:,2),T3,X3(:,2),T4,X4(:,2),T5,X5(:,2))
xlabel('Time(s)')
ylabel('x_2')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
subplot(313)
plot(T1,X1(:,3),T2,X2(:,3),T3,X3(:,3),T4,X4(:,3),T5,X5(:,3))
xlabel('Time(s)')
ylabel('x_3')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
% 显示欧拉
function [T,X,dX] = ODE_ExplicitEuler( 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
        dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
        h = T(k) - T(k-1);
        X(k,:) = X(k-1,:) + h*dX(k-1,:);                        
    end
    dX(N,:) = Hfun( T(N),X(N,:) );
end
% 隐式欧拉
function [T,X,dX] = ODE_ImplicitEuler( 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);
        for jj = 1:m
           Ind =  zeros(1,m);         % 状态选择向量
           Ind(jj) = 1;        % 选择的向量位置  
           fh = @(x) h*Ind*Hfun(T(k),[X(k-1,1:jj-1)';x; X(k-1,jj+1:m)']) + X(k-1,jj)' - x;    
           X(k,jj) = fzero( fh,X(k-1,jj) ) ;   
        end                  
    end
    dX(N,:) = Hfun( T(N),X(N,:) );
end

function [T,X,dX] = ODE_2OrderEuler( 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
        dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
        h = T(k) - T(k-1);
        if k == 2
            X(k,:) = X(k-1,:) + h*dX(k-1,:);    
        else
            X(k,:) = X(k-2,:) + 2*h*dX(k-1,:);   
        end
    end
    dX(N,:) = Hfun( T(N),X(N,:) );
end
function [T,X,dX] = ODE_ImprovedEuler( 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
        dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
        h = T(k) - T(k-1);
        Xp = X(k-1,:) + h*dX(k-1,:);
        dXp = Hfun( T(k),Xp );
        X(k,:) = X(k-1,:) + (h/2)*(dX(k-1,:)+dXp');
    end
    dX(N,:) = Hfun( T(N),X(N,:) );
end