Why ODE15s does not work like an Euler's method with fixed point?

When i run this code for Euler´s method, it runs Ok, but when i try to use ODE15s it shows a wrong answer. It is ODE15s well defined?. I'll appreciate your help.

function []=cr_NRTL(X0,x1,x2,l,L1)

global X1 X2 Xea estado
x1inicial=x1;
x2inicial=x2;

sol_daes=2; %1 Euler, 2 ODE
tspan=[0 l];
opts=odeset('NormControl','on','AbsTol',1e-10,'RelTol',1e-5);

for k=1:2
    x1=x1inicial;
    x2=x2inicial;
    for i=1:L1
        X1 = x1; X2 = x2;
        
        %Para solucionar la envolvente
        options = optimset('Display','off');
        lb=[0,0,0,0,0,0,0,0,300];
        ub=[1,1,1,1,1,1,1,1,400];
        Xea=lsqnonlin(@funcion,X0,lb,ub,options);
        
        X0=Xea;
        
        for j=1:9
            xea(i,j)=Xea(j);
        end
        
        %Para evaluar las residuales del ELV fuera de la envolvente
        if k==1 && (Xea(7)>0.9999 || Xea(8)>0.9999)
            if Xea(7)>1 || Xea(8)>1
                Xea
            end
            cr_NRTL2(x1,x2,l);
            break
        else
            
            xglobala=[x1; x2];
            
            %Para solucionar las residuales del ELLV
            switch sol_daes
                
                case 1 %Euler
                    switch k
                        case 1
                            tao=0.02;
                            dxdt=cambioELLV(xglobala);
                            xsa=[xglobala(1)+dxdt(1)*tao;
                                xglobala(2)+dxdt(2)*tao];
                        case 2
                            tao=0.1;
                            dxdt=cambioELLV(xglobala);
                            xsa=[xglobala(1)-dxdt(1)*tao;
                                xglobala(2)-dxdt(2)*tao];
                    end
                    
                case 2 %ODE
                    switch k
                        case 1
                            [~,dxdt]=ode15s(@cambioELLV,tspan,xglobala,opts);
                        case 2
                            [~,dxdt]=ode15s(@cambioELLV,-tspan,xglobala,opts);
                    end
                    xsa=dxdt;
                    
            end
            
            x1=xsa(1);
            x2=xsa(2);
            
            %                     if k==2
            for j=1:2
                xsalida(i,j)=xglobala(j);
                %                             end
            end
            
        end
        
        %             end
        %             if k ==2
        xe = xea;
        %             end
    end
    toc
    % Grafica
    hold on
    if estado==1
        plot(xe(:,1),xe(:,2),'r')
    end
    plot(xsalida(:,1),xsalida(:,2),'b');
end
end

[darova]

The code is too long. Do you have shorter code?

[darova]

Is it possible to see this?

[Cesar García Echeverry]

% function dxdt=cambioELLV(x) %for Euler
function dxdt=cambioELLV(~,x) %For ODE

global Xea
y=vlle(Xea);


dxdt=[x(1)-y(1);
    x(2)-y(2)];

end

When i use Euler method it works ok, if you have an idea, i can share with you all the codes by e-mail. Thank you.

[darova]

I can't handle it. It's to complicated

sorry

[Cesar García Echeverry]

Ok, thanks.