function xdot=shootsci(t,x,flag,par)

% This is the function that contains the differential equations for solving
% the IPV model optimally as well the components for running Newton's
% method.  This function is used by 'mainIPV.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);
x5=x(5);
p1=x(6);
p2=x(7);
p3=x(8);
p4=x(9);
p5=x(10);
%for Newton's Method
n11=x(11);
n12=x(12);
n13=x(13);
n14=x(14);
n15=x(15);
n21=x(16);
n22=x(17);
n23=x(18);
n24=x(19);
n25=x(20);
n31=x(21);
n32=x(22);
n33=x(23);
n34=x(24);
n35=x(25);
n41=x(26);
n42=x(27);
n43=x(28);
n44=x(29);
n45=x(30);
n51=x(31);
n52=x(32);
n53=x(33);
n54=x(34);
n55=x(35);
o11=x(36);
o12=x(37);
o13=x(38);
o14=x(39);
o15=x(40);
o21=x(41);
o22=x(42);
o23=x(43);
o24=x(44);
o25=x(45);
o31=x(46);
o32=x(47);
o33=x(48);
o34=x(49);
o35=x(50);
o41=x(51);
o42=x(52);
o43=x(53);
o44=x(54);
o45=x(55);
o51=x(56);
o52=x(57);
o53=x(58);
o54=x(59);
o55=x(60);

m=par(1);
a=par(2);
b=par(3);
c=par(4);
f=par(5);
g=par(6);   
h=par(7);
A=par(8);   %weighted parameters for Hamiltonian
B=par(9);

%define a switching function that turns control on and off (from optimal
%control theory)
sw=m*p1+a*m*(p2-p5);
if sw < 0 
    p=.90;
else if sw > 0 
        p=0;
    else p=.5;  % would be nice to have a 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=-b*(x3+g*x4)-m;
f12=0;
f13=-b*x1;
f14=-b*g*x1;
f15=0;
f21=0;
f22=-b*(x3+g*x4)*f-m;
f23=-b*f*x2;
f24=-b*g*f*x2;
f25=0;
f31=b*(x3+g*x4);
f32=0;
f33=b*x1-c-m;
f34=b*g*x1;
f35=0;
f41=0;
f42=b*(x3+g*x4)*f;
f43=b*f*x2;
f44=b*g*f*x2-c/h-m;
f45=0;
f51=0;
f52=0;
f53=c;
f54=c/h;
f55=-m;

%fijk is d^2fi/dxjdxk  used to compute variational equation for 
%Newton's method

f111=0;
f112=0;
f113=-b;
f114=-b*g;
f115=0;
f121=0;
f122=0;
f123=0;
f124=0;
f125=0;
f131=-b;
f132=0;
f133=0;
f134=0;
f135=0;
f141=-b*g;
f142=0;
f143=0;
f144=0;
f145=0;
f151=0;
f152=0;
f153=0;
f154=0;
f155=0;
f211=0;
f212=0;
f213=0;
f214=0;
f215=0;
f221=0;
f222=0;
f223=-b*f;
f224=-b*g*f;
f225=0;
f231=0;
f232=-b*f;
f233=0;
f234=0;
f235=0;
f241=0;
f242=-b*g*f;
f243=0;
f244=0;
f245=0;
f251=0;
f252=0;
f253=0;
f254=0;
f255=0;
f311=0;
f312=0;
f313=b;
f314=b*g;
f315=0;
f321=0;
f322=0;
f323=0;
f324=0;
f325=0;
f331=b;
f332=0;
f333=0;
f334=0;
f335=0;
f341=b*g;
f342=0;
f343=0;
f344=0;
f345=0;
f351=0;
f352=0;
f353=0;
f354=0;
f355=0;
f411=0;
f412=0;
f413=0;
f414=0;
f415=0;
f421=0;
f422=0;
f423=b*f;
f424=b*f*g;
f425=0;
f431=0;
f432=b*f;
f433=0;
f434=0;
f435=0;
f441=0;
f442=b*f*g;
f443=0;
f444=0;
f445=0;
f451=0;
f452=0;
f453=0;
f454=0;
f455=0;

%all second derivatives of x5 are zero 


%gi=-1*[p1 p2 p3 p4]*[f1i f2i f3i f4i]'; partial of H wrt state(i)
g1=-1*[p1 p2 p3 p4 p5]*[f11 f21 f31 f41 f51]';
g2=-1*[p1 p2 p3 p4 p5]*[f12 f22 f32 f42 f52]';
g3=-1*[p1 p2 p3 p4 p5]*[f13 f23 f33 f43 f53]'-A*x3;
g4=-1*[p1 p2 p3 p4 p5]*[f14 f24 f34 f44 f54]'-B*x4;
g5=-1*[p1 p2 p3 p4 p5]*[f15 f25 f35 f45 f55]';
%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 p5]*[f111 f211 f311 f411 0]';
g12=-1*[p1 p2 p3 p4 p5]*[f112 f212 f312 f412 0]';
g13=-1*[p1 p2 p3 p4 p5]*[f113 f213 f313 f413 0]';
g14=-1*[p1 p2 p3 p4 p5]*[f114 f214 f314 f414 0]';
g15=-1*[p1 p2 p3 p4 p5]*[f115 f215 f315 f415 0]';
g21=-1*[p1 p2 p3 p4 p5]*[f121 f221 f321 f321 0]';
g22=-1*[p1 p2 p3 p4 p5]*[f122 f212 f312 f322 0]';
g23=-1*[p1 p2 p3 p4 p5]*[f123 f223 f323 f323 0]';
g24=-1*[p1 p2 p3 p4 p5]*[f124 f224 f324 f324 0]';
g25=-1*[p1 p2 p3 p4 p5]*[f125 f225 f325 f325 0]';
g31=-1*[p1 p2 p3 p4 p5]*[f131 f231 f331 f431 0]';
g32=-1*[p1 p2 p3 p4 p5]*[f132 f232 f332 f432 0]';
g33=-1*[p1 p2 p3 p4 p5]*[f133 f233 f333 f433 0]'-A;
g34=-1*[p1 p2 p3 p4 p5]*[f134 f234 f334 f434 0]';
g35=-1*[p1 p2 p3 p4 p5]*[f135 f235 f335 f435 0]';
g41=-1*[p1 p2 p3 p4 p5]*[f141 f241 f341 f441 0]';
g42=-1*[p1 p2 p3 p4 p5]*[f142 f242 f342 f442 0]';
g43=-1*[p1 p2 p3 p4 p5]*[f143 f243 f343 f443 0]';
g44=-1*[p1 p2 p3 p4 p5]*[f144 f244 f344 f444 0]'-B;
g45=-1*[p1 p2 p3 p4 p5]*[f145 f245 f345 f445 0]';
g51=-1*[p1 p2 p3 p4 p5]*[f151 f251 f351 f451 0]';
g52=-1*[p1 p2 p3 p4 p5]*[f152 f252 f352 f452 0]';
g53=-1*[p1 p2 p3 p4 p5]*[f153 f253 f353 f453 0]';
g54=-1*[p1 p2 p3 p4 p5]*[f154 f254 f354 f454 0]';
g55=-1*[p1 p2 p3 p4 p5]*[f155 f255 f355 f455 0]';

xdot(1)=m*(1-p)-b*(x3+g*x4).*x1-m*x1;  %differential equations
xdot(2)=p*m*(1-a)-b*(x3+g*x4)*f.*x2-m*x2;
xdot(3)=b*(x3+g*x4).*x1-(c+m)*x3;
xdot(4)=b*(x3+g*x4)*f.*x2-(c/h+m)*x4;
xdot(5)=c*(x3+x4/h)-m*x5+m*a*p;
xdot(6)= -[p1 p2 p3 p4 p5]*[f11 f21 f31 f41 f51]';
xdot(7)= -[p1 p2 p3 p4 p5]*[f12 f22 f32 f42 f52]';
xdot(8)= -[p1 p2 p3 p4 p5]*[f13 f23 f33 f43 f53]'-A*x3;
xdot(9)= -[p1 p2 p3 p4 p5]*[f14 f24 f34 f44 f54]'-B*x4;
xdot(10)=-[p1 p2 p3 p4 p5]*[f15 f25 f35 f45 f55]';
xdot(11)=[f11 f12 f13 f14 f15]*[n11 n21 n31 n41 n51]'; 
xdot(12)=[f11 f12 f13 f14 f15]*[n12 n22 n32 n42 n52]';
xdot(13)=[f11 f12 f13 f14 f15]*[n13 n23 n33 n43 n53]';
xdot(14)=[f11 f12 f13 f14 f15]*[n14 n24 n34 n44 n54]';
xdot(15)=[f11 f12 f13 f14 f15]*[n15 n25 n35 n45 n55]';
xdot(16)=[f21 f22 f23 f24 f25]*[n11 n21 n31 n41 n51]'; 
xdot(17)=[f21 f22 f23 f24 f25]*[n12 n22 n32 n42 n52]';
xdot(18)=[f21 f22 f23 f24 f25]*[n13 n23 n33 n43 n53]';
xdot(19)=[f21 f22 f23 f24 f25]*[n14 n24 n34 n44 n54]';
xdot(20)=[f21 f22 f23 f24 f25]*[n15 n25 n35 n45 n55]';
xdot(21)=[f31 f32 f33 f34 f35]*[n11 n21 n31 n41 n51]'; 
xdot(22)=[f31 f32 f33 f34 f35]*[n12 n22 n32 n42 n52]';
xdot(23)=[f31 f32 f33 f34 f35]*[n13 n23 n33 n43 n53]';
xdot(24)=[f31 f32 f33 f34 f35]*[n14 n24 n34 n44 n54]';
xdot(25)=[f31 f32 f33 f34 f35]*[n15 n25 n35 n45 n55]';
xdot(26)=[f41 f42 f43 f44 f45]*[n11 n21 n31 n41 n51]'; 
xdot(27)=[f41 f42 f43 f44 f45]*[n12 n22 n32 n42 n52]';
xdot(28)=[f41 f42 f43 f44 f45]*[n13 n23 n33 n43 n53]';
xdot(29)=[f41 f42 f43 f44 f45]*[n14 n24 n34 n44 n54]';
xdot(30)=[f41 f42 f43 f44 f45]*[n15 n25 n35 n45 n55]';
xdot(31)=[f51 f52 f53 f54 f55]*[n11 n21 n31 n41 n51]'; 
xdot(32)=[f51 f52 f53 f54 f55]*[n12 n22 n32 n42 n52]';
xdot(33)=[f51 f52 f53 f54 f55]*[n13 n23 n33 n43 n53]';
xdot(34)=[f51 f52 f53 f54 f55]*[n14 n24 n34 n44 n54]';
xdot(35)=[f51 f52 f53 f54 f55]*[n15 n25 n35 n45 n55]';   %LET THE NEWTONING BEGIN!!!!
xdot(36)=[g11 g12 g13 g14 g15]*[n11 n21 n31 n41 n51]' - [f11 f21 f31 f41 f51]*[o11 o21 o31 o41 o51]'; 
xdot(37)=[g11 g12 g13 g14 g15]*[n12 n22 n32 n42 n52]' - [f11 f21 f31 f41 f51]*[o12 o22 o32 o42 o52]';  
xdot(38)=[g11 g12 g13 g14 g15]*[n13 n23 n33 n43 n53]' - [f11 f21 f31 f41 f51]*[o13 o23 o33 o43 o53]';
xdot(39)=[g11 g12 g13 g14 g15]*[n14 n24 n34 n44 n54]' - [f11 f21 f31 f41 f51]*[o14 o24 o34 o44 o54]';
xdot(40)=[g11 g12 g13 g14 g15]*[n15 n25 n35 n45 n55]' - [f11 f21 f31 f41 f51]*[o15 o25 o35 o45 o55]';
xdot(41)=[g21 g22 g23 g24 g25]*[n11 n21 n31 n41 n51]' - [f12 f22 f32 f42 f52]*[o11 o21 o31 o41 o51]';
xdot(42)=[g21 g22 g23 g24 g25]*[n12 n22 n32 n42 n52]' - [f12 f22 f32 f42 f52]*[o12 o22 o32 o42 o52]';
xdot(43)=[g21 g22 g23 g24 g25]*[n13 n23 n33 n43 n53]' - [f12 f22 f32 f42 f52]*[o13 o23 o33 o43 o53]';
xdot(44)=[g21 g22 g23 g24 g25]*[n14 n24 n34 n44 n54]' - [f12 f22 f32 f42 f52]*[o14 o24 o34 o44 o54]';
xdot(45)=[g21 g22 g23 g24 g25]*[n15 n25 n35 n45 n55]' - [f12 f22 f32 f42 f52]*[o15 o25 o35 o45 o55]';
xdot(46)=[g31 g32 g33 g34 g35]*[n11 n21 n31 n41 n51]' - [f13 f23 f33 f43 f53]*[o11 o21 o31 o41 o51]';
xdot(47)=[g31 g32 g33 g34 g35]*[n12 n22 n32 n42 n52]' - [f13 f23 f33 f43 f53]*[o12 o22 o32 o42 o52]';
xdot(48)=[g31 g32 g33 g34 g35]*[n13 n23 n33 n43 n53]' - [f13 f23 f33 f43 f53]*[o13 o23 o33 o43 o53]';
xdot(49)=[g31 g32 g33 g34 g35]*[n14 n24 n34 n44 n54]' - [f13 f23 f33 f43 f53]*[o14 o24 o34 o44 o54]';
xdot(50)=[g31 g32 g33 g34 g35]*[n15 n25 n35 n45 n55]' - [f13 f23 f33 f43 f53]*[o15 o25 o35 o45 o55]';
xdot(51)=[g41 g42 g43 g44 g45]*[n11 n21 n31 n41 n51]' - [f14 f24 f34 f44 f54]*[o11 o21 o31 o41 o51]';
xdot(52)=[g41 g42 g43 g44 g45]*[n12 n22 n32 n42 n52]' - [f14 f24 f34 f44 f54]*[o12 o22 o32 o42 o52]';
xdot(53)=[g41 g42 g43 g44 g45]*[n13 n23 n33 n43 n53]' - [f14 f24 f34 f44 f54]*[o13 o23 o33 o43 o53]';
xdot(54)=[g41 g42 g43 g44 g45]*[n14 n24 n34 n44 n54]' - [f14 f24 f34 f44 f54]*[o14 o24 o34 o44 o54]';
xdot(55)=[g41 g42 g43 g44 g45]*[n15 n25 n35 n45 n55]' - [f14 f24 f34 f44 f54]*[o15 o25 o35 o45 o55]';
xdot(56)=[g51 g52 g53 g54 g55]*[n11 n21 n31 n41 n51]' - [f15 f25 f35 f45 f55]*[o11 o21 o31 o41 o51]';
xdot(57)=[g51 g52 g53 g54 g55]*[n12 n22 n32 n42 n52]' - [f15 f25 f35 f45 f55]*[o12 o22 o32 o42 o52]';
xdot(58)=[g51 g52 g53 g54 g55]*[n13 n23 n33 n43 n53]' - [f15 f25 f35 f45 f55]*[o13 o23 o33 o43 o53]';
xdot(59)=[g51 g52 g53 g54 g55]*[n14 n24 n34 n44 n54]' - [f15 f25 f35 f45 f55]*[o14 o24 o34 o44 o54]';
xdot(60)=[g51 g52 g53 g54 g55]*[n15 n25 n35 n45 n55]' - [f15 f25 f35 f45 f55]*[o15 o25 o35 o45 o55]';


xdot=xdot';