# MATLAB: Confusion on a new schema for solving ode

solve ode numerically new schema

Hello all, recently i've been working on an algorithm for solving ode, called QSS family, which was found by Prof.Ernesto Kofman. I want to run the algorithm on matlab code, but i'm not sure if the code is right. Can anyone help me to tell if it's right?
The code below is an example to solve the ode:, with and all equal 0 at .
clcclearticdelta_q=1e-4;q1=0;q2=0;x1=0;x2=0;t=0;delta_t=0;tfinal=20;A=zeros(0,3);n=0;while (t<tfinal)    Dx1=q2;    Dx2=1-3*q1-4*q2;    if (Dx1>0)        delta_x1=delta_q+q1-x1;    else        delta_x1=delta_q-q1+x1;    end    if (Dx2>0)       delta_x2=delta_q+q2-x2;   else       delta_x2=delta_q-q2+x2;   end       delta_t1=delta_x1/abs(Dx1);       delta_t2=delta_x2/abs(Dx2);   if (delta_t1<delta_t2)       delta_t=delta_t1;       t=t+delta_t;       x1=x1+Dx1*delta_t;       x2=x2+Dx2*delta_t;       q1=x1;   else       delta_t=delta_t2;       t=t+delta_t;       x1=x1+Dx1*delta_t;       x2=x2+Dx2*delta_t;       q2=x2;   end   n=n+1;   A(n,1)=t;   A(n,2)=x1;   A(n,3)=x2;%    A(n,4)=q1;%    A(n,5)=q2;endA; figure(1) plot(A(:,1),A(:,2),A(:,1),A(:,3));toc       

function euler1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clcclearticdelta_q=1e-4;q1=0;q2=0;x1=0;x2=0;t=0;delta_t=0;tfinal=20;A=zeros(0,3);n=0;while (t<tfinal)    Dx1=q2;    Dx2=1-3*q1-4*q2;    if (Dx1>0)        delta_x1=delta_q+q1-x1;    else        delta_x1=delta_q-q1+x1;    end    if (Dx2>0)        delta_x2=delta_q+q2-x2;    else        delta_x2=delta_q-q2+x2;    end    delta_t1=delta_x1/abs(Dx1);    delta_t2=delta_x2/abs(Dx2);    if (delta_t1<delta_t2)        delta_t=delta_t1;        t=t+delta_t;        x1=x1+Dx1*delta_t;        x2=x2+Dx2*delta_t;        q1=x1;    else        delta_t=delta_t2;        t=t+delta_t;        x1=x1+Dx1*delta_t;        x2=x2+Dx2*delta_t;        q2=x2;    end    n=n+1;    A(n,1)=t;    A(n,2)=x1;    A(n,3)=x2;    %    A(n,4)=q1;    %    A(n,5)=q2;endfigure(1)close(1)figure(1)subplot(1,3,1)hold onaxis([0,20 -0.05 0.35])plot(A(:,1),A(:,2),'b')plot(A(:,1),A(:,3),'r');grid ontitle('HS')tocwhos%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t=linspace(0,20,501);x0=[0 0];whos[t,x]=ode45(@fun,t,x0);whossubplot(1,3,2)hold onaxis([0,20 -0.05 0.35])plot(t,x(:,1),'k')plot(t,x(:,2),'k')grid ontitle('edo45')subplot(1,3,3)hold onaxis([0,20 -0.05 0.35])plot(A(:,1),A(:,2),'b')plot(A(:,1),A(:,3),'r');plot(t,x(:,1),'k.')plot(t,x(:,2),'k.')grid ontitle('HS, edo45')    function [dxdt]=fun(~,x)        x1p=x(2);        x2p=1-3*x(1)-4*x(2);        dxdt=[x1p; x2p];    endend