function task3_2()    
    time = 0:1:200; % ODEシミュレーションの時間区間
    x_ini = 0; % xの初期値
 【ここを埋める】 % パラメータ
    figure;
    hold on;
    for x_out = 0:0.5:10
        param = [ alpha , beta , rho , x_out ]; % パラメータベクトル
        for y_ini = 0:0.1:1.0
            s0 = [ x_ini , y_ini ]; % 初期値ベクトル
            % ODEシミュレーション
            [~, time_course] = ode45(@(t, x) ODE(t,x,param), time, s0);
            y_end = time_course(end,2); % 最終時刻のyの値を格納
            % 縦軸 yの最終値、横軸 x_outでプロット
            plot(x_out, y_end, 'r+', 'markersize', 10, 'linewidth', 2);
        end
    end
	xlabel('TMGout (Xout)'); % 横軸のラベル名
	ylabel('Final state of LacY-GFP (x1)'); % 縦軸のラベル名

    for x_out = 0:0.5:10   
        syms x y;
        % 代数方程式を解く
        [~,ysol] = vpasolve([【ここを埋める】,x^2/(rho+x^2)-y==0],[x,y]);
        for i=1:size(ysol)
            if(isreal(ysol(i))) % 実数解のみを表示させる
                plot(x_out, ysol(i), 'ko', 'markersize', 10, 'linewidth', 2);
            end
        end
    end 
    
end

function dsdt = ODE(~,s,param) % ODE(t,s,param)という関数の定義
 【ここを埋める】
end