function xdot=shootsci(t,x,flag,par)
%this is a sub-program for the OPV model.  It contains the differential
%equations and the components used for applying Newton's Method.  One would
%never call this program.  It is only used by the function mainsci.m.

%import parameter values in par 
%the state variables are x1-susceptible, x2-vaccinated, x3-wild-virus carriers, x4-removed
%p1, p2, p3, p4 are the costate variables from optimal control theory
%m is both the birth and death rate
%a is the rate a person carrying the wild-virus contacts other individuals
%b is similar to a, but for those carrying the vaccine virus
%c and d are the loss of vaccine virus and wild-virus infectivity (resp.)
%sw is the switching function: it tells p when to turn on

%make things easier to input
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
p1=x(5);
p2=x(6);
p3=x(7);
p4=x(8);
%for Newton's Method
n11=x(9);
n12=x(10);
n13=x(11);
n14=x(12);
n21=x(13);
n22=x(14);
n23=x(15);
n24=x(16);
n31=x(17);
n32=x(18);
n33=x(19);
n34=x(20);
n41=x(21);
n42=x(22);
n43=x(23);
n44=x(24);
o11=x(25);
o12=x(26);
o13=x(27);
o14=x(28);
o21=x(29);
o22=x(30);
o23=x(31);
o24=x(32);
o31=x(33);
o32=x(34);
o33=x(35);
o34=x(36);
o41=x(37);
o42=x(38);
o43=x(39);
o44=x(40);

m=par(1);
a=par(2);
b=par(3);
c=par(4);
d=par(5);
A=par(6);   %weighted parameters from Hamiltonian
B=par(7);


%define a switching function that turns control on and off (from optimal
%control theory)
sw=m*(p1-p2);
if sw < 0 
    p=.9;
else if sw > 0 
        p=0;
    else p=.5;  % would be nice to have smooth function, but this is bang-bang
    end
end


% these are the partial derivatives of the state variable derivatives to be used later
%fij is dfi/dxj
f11=-a*x3-b*x2-m;
f12=-b*x1;
f13=-a*x1;
f14=0;
f21=b*x2;
f22=b*x1-c-m;
f23=0;
f24=0;
f31=a*x3;
f32=0;
f33=a*x1-d-m;
f34=0;
f41=0;
f42=c;
f43=d;
f44=-m;
%fijk is d^2fi/dxjdxk  used to compute variational equation for Newton's
%method
f111=0;
f112=-b;
f113=-a;
f114=0;
f121=-b;
f122=0;
f123=0;
f124=0;
f131=-a;
f132=0;
f133=0;
f134=0;
f141=0;
f142=0;
f143=0;
f144=0;
f211=0;
f212=b;
f213=0;
f214=0;
f221=b;
f222=0;
f223=0;
f224=0;
f231=0;
f232=0;
f233=0;
f234=0;
f241=0;
f242=0;
f243=0;
f244=0;
f311=0;
f312=0;
f313=a;
f314=0;
f321=0;
f322=0;
f323=0;
f324=0;
f331=a;
f332=0;
f333=0;
f334=0;
f341=0;
f342=0;
f343=0;
f344=0;
%all second derivatives of x4 are zero 


%gi=-1*[p1 p2 p3 p4]*[f1i f2i f3i f4i]'; partial of H wrt state(i)
g1=-1*[p1 p2 p3 p4]*[f11 f21 f31 f41]';
g2=-1*[p1 p2 p3 p4]*[f12 f22 f32 f42]'-A*x2;
g3=-1*[p1 p2 p3 p4]*[f13 f23 f33 f43]'-B*x3;
g4=-1*[p1 p2 p3 p4]*[f14 f24 f34 f44]';
%gik=dgi/dxk = [p1 p2 p3 p4]*[d^2f1/dxidxk d^2f2/dxidxk d^2f2/dxidxk d^2f4/dxidxk]' 
%second partial of H wrt state
g11=-1*[p1 p2 p3 p4]*[f111 f211 f311 0]';
g12=-1*[p1 p2 p3 p4]*[f112 f212 f312 0]';
g13=-1*[p1 p2 p3 p4]*[f113 f213 f313 0]';
g14=-1*[p1 p2 p3 p4]*[f114 f214 f314 0]';
g21=-1*[p1 p2 p3 p4]*[f121 f221 f321 0]';
g22=-1*[p1 p2 p3 p4]*[f122 f212 f312 0]'-A;
g23=-1*[p1 p2 p3 p4]*[f123 f223 f323 0]';
g24=-1*[p1 p2 p3 p4]*[f124 f224 f324 0]';
g31=-1*[p1 p2 p3 p4]*[f131 f231 f331 0]';
g32=-1*[p1 p2 p3 p4]*[f132 f232 f332 0]';
g33=-1*[p1 p2 p3 p4]*[f133 f233 f333 0]'-B;
g34=-1*[p1 p2 p3 p4]*[f134 f234 f334 0]';
g41=-1*[p1 p2 p3 p4]*[f141 f241 f341 0]';
g42=-1*[p1 p2 p3 p4]*[f142 f242 f342 0]';
g43=-1*[p1 p2 p3 p4]*[f143 f243 f343 0]';
g44=-1*[p1 p2 p3 p4]*[f144 f244 f344 0]';

xdot(1)=(1-p)*m-a*x1.*x3-b*x1.*x2-m*x1;  %differential equations
xdot(2)=p*m+b*x1.*x2-c*x2-m*x2;
xdot(3)=a*x1.*x3-d*x3-m*x3;
xdot(4)=d*x3+c*x2-m*x4;
xdot(5)=-[p1 p2 p3 p4]*[f11 f21 f31 f41]';
xdot(6)=-[p1 p2 p3 p4]*[f12 f22 f32 f42]'-A*x2;
xdot(7)=-[p1 p2 p3 p4]*[f13 f23 f33 f43]'-B*x3;
xdot(8)=-[p1 p2 p3 p4]*[f14 f24 f34 f44]';
xdot(9) =[f11 f12 f13 f14]*[n11 n21 n31 n41]'; 
xdot(10)=[f11 f12 f13 f14]*[n12 n22 n32 n42]';
xdot(11)=[f11 f12 f13 f14]*[n13 n23 n33 n43]';
xdot(12)=[f11 f12 f13 f14]*[n14 n24 n34 n44]';
xdot(13)=[f21 f22 f23 f24]*[n11 n21 n31 n41]'; 
xdot(14)=[f21 f22 f23 f24]*[n12 n22 n32 n42]';
xdot(15)=[f21 f22 f23 f24]*[n13 n23 n33 n43]';
xdot(16)=[f21 f22 f23 f24]*[n14 n24 n34 n44]';
xdot(17)=[f31 f32 f33 f34]*[n11 n21 n31 n41]'; 
xdot(18)=[f31 f32 f33 f34]*[n12 n22 n32 n42]';
xdot(19)=[f31 f32 f33 f34]*[n13 n23 n33 n43]';
xdot(20)=[f31 f32 f33 f34]*[n14 n24 n34 n44]';
xdot(21)=[f41 f42 f43 f44]*[n11 n21 n31 n41]'; 
xdot(22)=[f41 f42 f43 f44]*[n12 n22 n32 n42]';
xdot(23)=[f41 f42 f43 f44]*[n13 n23 n33 n43]';
xdot(24)=[f41 f42 f43 f44]*[n14 n24 n34 n44]'; %BEGIN NEWTON!!!!
xdot(25)=[g11 g12 g13 g14]*[n11 n21 n31 n41]' - [f11 f21 f31 f41]*[o11 o21 o31 o41]'; 
xdot(26)=[g11 g12 g13 g14]*[n12 n22 n32 n42]' - [f11 f21 f31 f41]*[o12 o22 o32 o42]';  
xdot(27)=[g11 g12 g13 g14]*[n13 n23 n33 n43]' - [f11 f21 f31 f41]*[o13 o23 o33 o43]';
xdot(28)=[g11 g12 g13 g14]*[n14 n24 n34 n44]' - [f11 f21 f31 f41]*[o14 o24 o34 o44]';
xdot(29)=[g21 g22 g23 g24]*[n11 n21 n31 n41]' - [f12 f22 f32 f42]*[o11 o21 o31 o41]';
xdot(31)=[g21 g22 g23 g24]*[n12 n22 n32 n42]' - [f12 f22 f32 f42]*[o12 o22 o32 o42]';
xdot(31)=[g21 g22 g23 g24]*[n13 n23 n33 n43]' - [f12 f22 f32 f42]*[o13 o23 o33 o43]';
xdot(32)=[g21 g22 g23 g24]*[n14 n24 n34 n44]' - [f12 f22 f32 f42]*[o14 o24 o34 o44]';
xdot(33)=[g31 g32 g33 g34]*[n11 n21 n31 n41]' - [f13 f23 f33 f43]*[o11 o21 o31 o41]';
xdot(34)=[g31 g32 g33 g34]*[n12 n22 n32 n42]' - [f13 f23 f33 f43]*[o12 o22 o32 o42]';
xdot(35)=[g31 g32 g33 g34]*[n13 n23 n33 n43]' - [f13 f23 f33 f43]*[o13 o23 o33 o43]';
xdot(36)=[g31 g32 g33 g34]*[n14 n24 n34 n44]' - [f13 f23 f33 f43]*[o14 o24 o34 o44]';
xdot(37)=[g41 g42 g43 g44]*[n11 n21 n31 n41]' - [f14 f24 f34 f44]*[o11 o21 o31 o41]';
xdot(38)=[g41 g42 g43 g44]*[n12 n22 n32 n42]' - [f14 f24 f34 f44]*[o12 o22 o32 o42]';
xdot(39)=[g41 g42 g43 g44]*[n13 n23 n33 n43]' - [f14 f24 f34 f44]*[o13 o23 o33 o43]';
xdot(40)=[g41 g42 g43 g44]*[n14 n24 n34 n44]' - [f14 f24 f34 f44]*[o14 o24 o34 o44]';

xdot=xdot';