function task3_1()    
    time = 0:1:100;
    x_ini = 0; 
　　【ここを埋める】% 各パラメータの値
    param = [ alpha , beta , rho , x_out ]; % パラメータ・ベクトル
    for y_ini = 0:0.33:1.0 % yの初期値を変動
        s0 = [ x_ini , y_ini ];
    	[t, time_course] = ode45(@(t, s) ODE(t,s,param), time, s0);
        figure; 
        hold on;
        plot(t, time_course(:,1), 'r'); % 横軸:時間、縦軸:xでplot
        plot(t, time_course(:,2), 'b'); % 横軸:時間、縦軸:yでplot
        ylim([0,10])
        xlabel('Time (t)'); % 横軸のラベル名
        ylabel('Cocentration'); % 縦軸のラベル名
        legend({'TMGin (x)','LacY (y)'}); % 各系列の説明
        title(join(["y(0) = ",y_ini])); % タイトルのラベル名
    end
end

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