function ozbudak_fig2b
    time = 0:0.5:30;

    x_init = 0;   y_init = 0; 
    alpha = 0.1;  beta = 10;  rho = 25; x_out = 0; 
    
    %s0 = [ x_init , y_init ];    % Initial values
    param = [ alpha , beta , rho , x_out ]; % Constants
    
    figure(1)
    hold on
    
    while x_out <= 20;
        s0 = [ x_init , y_init ];    % Initial values
        [t,time_course] = ode15s(@(t,s) ODE(t,s,param),time,s0); 
        y_end = �@�@�@�@�@�@�@�@
        plot( x_out , y_end , 'ko');        
        x_out = x_out + 1;
        
        param = [ alpha , beta , rho , x_out ]; % Constants
        
    end
     
    while x_out >= 0
        �@�@�@�@�@�@�@�@
		�@�@�@�@�@�@�@�@
    end    
end

function dsdt = ODE(t,s,param)
    x = s(1);  y = s(2); 
    alpha = param(1);  beta = param(2); rho = param(3); x_out = param(4);
    
    dsdt(1,:) = - x + alpha * x_out + beta * y;
    dsdt(2,:) = x^2  / ( rho + x^2 ) - y ;
end
