### I am writing a code for an Optimal control using pontryagins maximum principle. I have supplied all neccessary boundary conditions yet it fails to detect one of the conditions. It gives an error message whenever I run it.

I need guidance. Thank you

restart;
with(plots);
r := 3; r[1] := 3; k := 10; a := 0.2e-1; b := 0.1e-1; c := 0.1e-1; beta := 0.3e-1; alpha := 0.3e-1; m := 0.5e-1;
z := 40; q := 5; p := 100; T := 3;
sigma := 0.1e-1; k[1] := 10; rho := 0.5e-1;

u[1] := min(max(0, z), 1); z := (a*m*k*lambda[2](t)*x(t)*y(t)-lambda[1](t)*r*(1+b*x(t)+c*y(t))*x(t)*x(t))/(z*k*(1+b*x(t)+c*y(t))); u[2] := min(max(0, q), 1); q := -lambda[1](t)*beta*x(t)*s(t)/q; u[3] := min(max(0, p), 1); p := -(r[1]*lambda[3](t)*s(t)*s(t))/(p*k[1]);
NULL;
sys := diff(x(t), t) = r*x(t)*(1-(1-u[1])*x(t)/k)-a*m*x(t)*y(t)/(1+b*x(t)+c*y(t))-beta*(1-u[2])*x(t)*s(t), diff(y(t), t) = -alpha*y(t)+a*m*x(t)*y(t)/(1+b*x(t)+c*y(t)), diff(s(t), t) = sigma*s(t)+r[1]*s(t)*(1-(1-u[3])*s(t)/k[1])-rho*s(t)*y(t), diff(lambda[1](t), t) = -lambda[1](t)*(r-2*r*(1-u[1])*x(t)/k-a*y(t)*(1+c*y(t))/((1+b*x(t)+c*y(t)) . (1+b*x(t)+c*y(t)))-beta*(1-u[2])*s(t))-lambda[2](t)*a*m*(1-u[1])*(1+c*y(t))*y(t)/((1+b*x(t)+c*y(t)) . (1+b*x(t)+c*y(t))), diff(lambda[2](t), t) = -lambda[1](t)*a*x(t)*(1+b*x(t))/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t)))+lambda[2](t) . (-alpha+(a*m*(1-u[1]) . (1+b*x(t)))*x(t)/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t))))+lambda[3](t)*rho*s(t), diff(lambda[3](t), t) = lambda[1](t)*beta*(1-u[2])*x(t)-lambda[1](t)*(r[1]-2*r[1]*(1-u[3])*s(t)/k[1]-sigma-rho*y(t)), x(0) = 10, y(0) = 20, s(0) = 10, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0;
p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = .1);
Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 7, got 6
p2o := odeplot(p1, [t, y(t)], 0 .. 2, numpoints = 100, labels = ["Time (months)", " "*`badbiomass""spatina"""`], labeldirections = [horizontal, vertical], style = line, color = red, axes = boxed);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
r := 3; r[1] := 3; k := 10; a := 0.2e-1; b := 0.1e-1; c := 0.1e-1; beta := 0.3e-1; alpha := 0.3e-1; m := 0.5e-1; sigma := 0.1e-1; k[1] := 10; rho := 0.5e-1;
z := 40; q := 5; p := 100; T := 3;

fun := diff(x(t), t) = r*x(t)*(1-x(t)/k)-a*m*x(t)*y(t)/(1+b*x(t)+c*y(t))-beta*x(t)*s(t), diff(y(t), t) = -alpha*y(t)+a*m*x(t)*y(t)/(1+b*x(t)+c*y(t)), diff(s(t), t) = sigma*s(t)+r[1]*s(t)*(1-s(t)/k[1])-rho*s(t)*y(t), x(0) = 10, y(0) = 20, s(0) = 10;
p2 := dsolve({fun}, type = numeric);
Error, (in f) unable to store 'HFloat(10.0)*r*(1-HFloat(10.0)/k)-HFloat(3.1538461538461537)' when datatype=float[8]
p2i := odeplot(p2, [t, y(t)], 0 .. 2, numpoints = 100, labels = ["Time(months)", " bad biomass"], labeldirections = [horizontal, vertical], style = line, axes = boxed, color = blue);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
plots[display](p2i, p2o);
Error, (in plots:-display) expecting plot structure but received: p2i

I am carrying out a research in dynamical system to which end I need to do optimal control. I have coded the control equations but its not displaying any result. I need guidance. The code is shown below....

restart;
with(plots);
r := 3; r[1] := 3; k := 10; a := 0.2e-1; b := 0.1e-1; c := 0.1e-1; beta := 0.3e-1; alpha := 0.3e-1; m := 0.5e-1;
z := 40; q := 5; p := 100; T := 3;
sigma := 0.1e-1; k[1] := 10; rho := 0.5e-1;

u[1] := min(max(0, z), 1); z := (a*m*k*lambda[2](t)*x(t)*y(t)-lambda[1](t)*r*(1+b*x(t)+c*y(t))*x(t)*x(t))/(z*k*(1+b*x(t)+c*y(t))); u[2] := min(max(0, q), 1); q := -lambda[1](t)*beta*x(t)*s(t)/q; u[3] := min(max(0, p), 1); p := -(r[1]*lambda[3](t)*s(t)*s(t))/(p*k[1]);
NULL;
sys := diff(x(t), t) = r*x(t)*(1-(1-u[1])*x(t)/k)-a*m*x(t)*y(t)/(1+b*x(t)+c*y(t))-beta*(1-u[2])*x(t)*s(t), diff(y(t), t) = -alpha*y(t)+a*m*x(t)*y(t)/(1+b*x(t)+c*y(t)), diff(s(t), t) = sigma*s(t)+r[1]*s(t)*(1-(1-u[3])*s(t)/k[1])-rho*s(t)*y(t), diff(lambda[1](t), t) = -lambda[1](t)*(r-2*r*(1-u[1])*x(t)/k-a*y(t)*(1+c*y(t))/((1+b*x(t)+c*y(t)) . (1+b*x(t)+c*y(t)))-beta*(1-u[2])*s(t))-lambda[2](t)*a*m*(1-u[1])*(1+c*y(t))*y(t)/((1+b*x(t)+c*y(t)) . (1+b*x(t)+c*y(t))), diff(lambda[2](t), t) = -lambda[1](t)*a*x(t)*(1+b*x(t))/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t)))+lambda[2](t) . (-alpha+(a*m*(1-u[1]) . (1+b*x(t)))*x(t)/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t))))+lambda[3](t)*rho*s(t), diff(lambda[3](t), t) = lambda[1](t)*beta*(1-u[2])*x(t)-lambda[1](t)*(r[1]-2*r[1]*(1-u[3])*s(t)/k[1]-sigma-rho*y(t)), x(0) = 100, y(0) = 200, s(0) = 100, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0;
p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = .1);
 

I am working on a project in dynamical system and I am to apply Pontryagins maximum principle of control. i have coded the control equations but it shows an error message "Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 7, got 6". What could be the possible problem?