Nonlinear ODE minimization - boundary value problem

Hello, I'm trying to implement the system used in chemotherapy drug administration, suggested in the paper (see attachment) by L.G.dePillis & A.Radunskaya.

The equaitons are :

N_dot = r2*N*(1 -b2*N) -c4*T*N -a3*(1-exp(-u))*N;
T_dot = r1*T*(1 -b1*T) -c2*I*T -c3*T*N -a2*(1-exp(-u))*T;
I_dot = s +ro*I*T/(alpha +T) -c1*I*T -d1*I -a1*(1-exp(-u))*I;
u_dot = v - d2*u;

The objective function to be minimized at the final time (tf) is :

J(tf) = T(tf)

under the constraint :

N >= 0.75; % for every t

Also in p.13, there is an analysis of the Hamiltonian and the co-state variables p1,p2,p3,p4.

Their boundary values for tf are :

p1(tf) = 0;
p2(tf) = 1;
p3(tf) = 0;
p4(tf) = 0;

Then the control equation is p4 and the drug (v) is administered as follows :

if p4>0
    v = 0;
elseif p4<0
    v = 1; % maximum dosage
end

The paper says this is a two point boundary value problem.

How can I implement it in Matlab (or Simulink)? I think some information is missing to use bvp4c.

Thank you.