function ozbudak_vector_field( )

    time = 0:0.1:60;

    x_init = 6; % TMG (something like IPTG)       
    y_init = 0; % Fluorescence of tagged LacY permease
        
    alpha = 0.1; 
    beta = 10; 
    rho = 25; 
    x_out = 5;
    
    param = [ alpha , beta , rho , x_out ]; % Constants
    s0 = [ x_init , y_init ];    % Initial values
      
    [t,time_course] = ode15s(@(t,s) ODE(t,s,param),time,s0); 
    figure(2)
    plot(t,time_course)
    plot_phase_plane(time_course , alpha , beta , rho, x_out);
    vector_field( alpha, beta , rho , x_out )
end

function dsdt = ODE(t,s,param)
    x = s(1); % TMG (something like IPTG)       
    y = s(2); % Fluorescence of tagged LacY permease
    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

function plot_phase_plane( time_course , alpha , beta , rho , x_out)
    x_time_course = time_course(:,1);
    y_time_course = time_course(:,2);

    figure(1);
    plot( x_time_course , y_time_course , 'k');
    hold on;
    xlabel('TMG')
    ylabel('LacY-GFP')
    
    y = 0:0.01:0.8;
    x = beta * y + alpha * x_out; % d TMG / dt == 0
    plot( x , y , 'r' );
    
    x = 0:0.1:10;
    y = x.^2 ./ ( rho + x.^2 ); % d LacY_GFP / dt == 0
    plot( x , y , 'b' );
end

function vector_field( alpha , beta , rho , x_out )    
	[x,y] = meshgrid(�@�@�@�@�@�@�@�@, �@�@�@�@�@�@�@�@);    
	DX = �@�@�@�@�@�@�@�@�@�@�@
	DY = �@�@�@�@�@�@�@�@�@�@�@
	figure(1)    
	quiver(�@�@�@�@�@�@�@, �@�@�@�@�@�@�@�@,�@�@�@�@�@�@�@�@,�@�@�@�@�@�@�@�@);
end
