How to plot u1 and u2 in program ???

I want to plot graph u1 and u2. But I can't plot. And I'm not sure write u1 and u2 is a function are correct.

function optimal_control_of_the_customer_dynamics
close all
clear
solinit = bvpinit(linspace(0,7,100),[600 900 700 500 1000 300]);
sol = bvp4c(@marketing_ode,@marketing_bc,solinit);
x = linspace(0,7);
y = deval(sol,x);
sol = bvp4c(@marketing_ode1,@marketing_bc,solinit);
x = linspace(0,7);
y1 = deval(sol,x);
figure()
plot(x,y1(3,:),'LineWidth',2.5) ,hold on
plot(x,y(3,:),'LineWidth',2.5)
plot(x,y(2,:),'LineWidth',2.5)
plot(x,y1(2,:),'LineWidth',2.5)
title('Evolution of C and P')
legend('P w/o c.','P with c.','C with c.','C w/o c.','Location','Best')
figure()
plot(x,y(1,:),'LineWidth',2.5) ,hold on
plot(x,y1(1,:),'LineWidth',2.5)
axis([0 7 0 .02])
title('Evolution of R')
legend('R w/o c.','R with c.','Location','Best')

% -------------------------------------------------------------
% -------------------------- Function  ------------------------
% -------------------------------------------------------------

function dydx = marketing_ode1(x,y)
global N bet gam u1 u2 P1 P2 P3 k1 k2 k3 lam1 lam2 alp1 alp2;
N = y(1)+y(2)+y(3);
bet = 0.01+(0.99/(1+exp(-2*x+8)));
gam = 0.10;
lam1 = 0.002;
lam2 = 0.018;
alp1 = 0.05;
alp2 = 0.1;
k1 = 1;
k2 = 1.5;
k3 = 0.01;
u1=0;
u2=0;
P1 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(3)*(y(2)+y(3)))/N^2;
P2 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(3)*y(1))/N^2;
P3 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(1)*(y(2)+y(1)))/N^2;
dydx = [-lam2*y(1)+lam1*y(2)-gam*y(1)+alp1*u1*y(3)+(alp2*(bet+u2)*y(1)*y(3))/N
    -lam1*y(2)+lam2*y(1)-gam*y(2)+((1-alp2)*(bet+u2)*y(1)*y(3))/N+u1*(1-alp1)*y(3)
    (-(bet+u2)*y(1)*y(3))/N-u1*y(3)+gam*y(1)+gam*y(2)
    lam2*(y(4)-y(5))+gam*(y(4)-y(6))+P1
    lam1*(y(5)-y(4))+gam*(y(5)-y(6))-P2
    -k1+(y(6)-(alp1*y(4))-(1-alp1)*y(5))*u1*P3];

% -------------------------------------------------------------
% ---------------------- Function control ---------------------
% -------------------------------------------------------------

function dydx = marketing_ode(x,y)
global N bet gam u1 u2 P1 P2 P3 k1 k2 k3 lam1 lam2 alp1 alp2;
N = y(1)+y(2)+y(3);
bet = 0.01+(0.99/(1+exp(-2*x+8)));
gam = 0.10;
lam1 = 0.002;
lam2 = 0.018;
alp1 = 0.05;
alp2 = 0.1;
k1 = 1;
k2 = 1.5;
k3 = 0.01;
u1=min(max(0,((y(6)-(alp1*y(4))-(1-alp1)*y(5))*y(3))/(2*k2)),0.06);
u2=min(max(0,((y(6)-(alp2*y(4))-(1-alp2)*y(5))*y(3)*y(1))/(2*k3*N)),1.0);
P1 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(3)*(y(2)+y(3)))/N^2;
P2 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(3)*y(1))/N^2;
P3 = ((y(6)-(alp2*y(4))-(1-alp2)*y(5))*(bet+u2)*y(1)*(y(2)+y(1)))/N^2;
dydx = [-lam2*y(1)+lam1*y(2)-gam*y(1)+alp1*u1*y(3)+(alp2*(bet+u2)*y(1)*y(3))/N
    -lam1*y(2)+lam2*y(1)-gam*y(2)+((1-alp2)*(bet+u2)*y(1)*y(3))/N+u1*(1-alp1)*y(3)
    (-(bet+u2)*y(1)*y(3))/N-u1*y(3)+gam*y(1)+gam*y(2)
    lam2*(y(4)-y(5))+gam*(y(4)-y(6))+P1
    lam1*(y(5)-y(4))+gam*(y(5)-y(6))-P2
    -k1+(y(6)-(alp1*y(4))-(1-alp1)*y(5))*u1*P3];
% -------------------------------------------------------------
% -------------------------- Boundary -------------------------
% -------------------------------------------------------------

function res = marketing_bc(ya,yb)
res = [ya(1)-0.001
    ya(2)-0.009
    ya(3)-0.99
    yb(4)-0
    yb(5)-0
    yb(6)-0];

[KSSV]

But seems u1 and u2 are constants....?

[thiti prasertjitsun]

Oh!!!, I want u1 and u2 is a function but I'm not sure write is correct.