> 
restart;
with(plots): with(DEtools): with(LinearAlgebra): with(Statistics): with(CurveFitting): with(Optimization):

> 
iV := Vector([300870, 83142, 155094, 146394, 177041, 94876, 150116, 158735, 228170, 144377, 167341, 44166, 136422, 125252, 161935, 138512], datatype = float);


(1) 

(2) 
> 
tV := Vector([seq(j, j = 1 .. len)])


(3) 
> 
pp := plot(tV, iV, style = point):

> 
eq1 := diff(B(T), T) = M__hbeta__1*psi*(B(T)+C(T))*H(T)sigma*psi*beta__1*E(T)*DD(T)mu__r*B(T)


(4) 
> 
eq2 := diff(C(T), T) = beta__1*psi*(B(T)+C(T))*H(T)(alpha+xi+mu__r)*C(T)


(5) 
> 
eq3 := diff(DD(T), T) = alpha*C(T)(varpi+xi+mu__r)*DD(T)


(6) 
> 
eq4 := diff(E(T), T) = varpi*DD(T)(Ggamma+mu__r)*E(T)


(7) 
> 
eq5 := diff(F(T), T) = Ggamma*E(T)+sigma*psi*beta__1*F(T)*H(T)mu__r*F(T)


(8) 
> 
eq6 := diff(G(T), T) = M__cpsi*beta__o*G(T)*DD(T)mu__b*G(T)


(9) 
> 
eq7 := diff(H(T), T) = psi*beta__o*G(T)*DD(T)mu__b*H(T)


(10) 
> 
B0 := .50:
C0 := .30:
DD0 := 0.21:
E0 := 0.14:
F0 := 0.7:
G0 := 0.45:
H0 := 0.14:

> 
sys := {eq1, eq2, eq3, eq4, eq5, eq6, eq7, B(0) = B0, C(0) = C0, DD(0) = DD0, E(0) = E0, F(0) = F0, G(0) = G0, H(0) = H0}


(11) 
> 
FreeParams := convert(indets(sys, name) minus {T}, list)


(12) 
> 
N := dsolve(
sys,
numeric,
parameters = FreeParams,
output = listprocedure
);


(13) 
> 
ParamValues := [M__h = .352, beta__o = 0.79e3, beta__1 = 0.198e2, mu__r = 0.96e2, sigma = .7902, alpha = .65, psi = 0.2e2, xi = 0.5e1, Ggamma = .7, M__c = .144, mu__b = 0.68e1, varpi = .134];


(14) 
> 
# check if the params in ParamValues are those in FreeParams
evalb( {FreeParams[]} = {lhs~(ParamValues)[]} );
{FreeParams[]} ; {lhs~(ParamValues)[]} ;


(15) 
> 
# instanciate N:
N(parameters=subs(ParamValues, FreeParams)):

> 
# example :
b := subs(N, B(T)):
plot(b(T), T=tV[1]..tV[len]);

> 
ip := subs(N, i(t)): # I don't understand what you want to do

> 
SSQ := proc(A, B)
subs({alpha=A, beta__1=B}, FreeParams);
subs(ParamValues, %);
N(parameters=%):
add((b(tV(j))iV[j])^2, j=1..len) # I replace arbitrarily ip by B to show you how to proceed
end proc:


(16) 
> 
plot3d('SSQ(A, B)', A=0 .. 0.2e3, B=0 .. .1, axes = boxed, labels = [beta, alpha, SSQ])

