If I1, I2, I3 and I4 are "knowns", you shouldn't assign values to them as you did.  I think you want something like this.  I'll call the values i1, i2, i3, i4.

> solve({I1=i1,I2=i2,I3=i3,I4=i4},{M,a,b,theta},AllSolutions);

Unfortunately this doesn't produce a result, so I'll try Groebner basis methods.

> with(Groebner);
G:= map(numer,subs(cos(theta)=c,sin(theta)=s,convert([I1-i1,I2-i2,I3-i3,I4-i4,s^2+c^2-1],trig)));
B:= Basis(G,plex(a,b,c,s,M,i1,i2,i3,i4));
B[1];

i4*i3^2*i2^2-i3^2*i2^2*i1-i2^2*i1*M^2+i2*i1^2*M^2-i4^2*i3*i1^2+i4^2*i2*i1^2+i4*i3^2*M^2+i3*i2^2*M^2*I+i4^2*i1*M^2*I+i4^2*i3^2*i1*I-I*i4*i2^2*i1^2-I*i4*i1^2*M^2-I*i4^2*i3^2*i2-I*i3^2*i2*M^2+i3*i2^2*i1^2*I-i4^2*i3*M^2+(1-I)*i4*i3^2*i2*i1+(-1+I)*i4*i3*i2^2*i1+(1-I)*i4*i3*i2*i1^2+(-2+2*I)*i4*i3^2*i2*M-(1+I)*i4*i3*i2^2*M-(1+I)*i4*i3^2*i1*M+(1+I)*i3^2*i2*i1*M+(1+I)*i4*i2^2*i1*M+(2-2*I)*i3*i2^2*i1*M+(1+I)*i4*i3*i1^2*M+(-2+2*I)*i4*i2*i1^2*M-(1+I)*i3*i2*i1^2*M+(1-I)*i4*i3*i2*M^2+(-1+I)*i4*i3*i1*M^2+(1-I)*i4*i2*i1*M^2+(-1+I)*i3*i2*i1*M^2+(-1+I)*i4^2*i3*i2*i1+(1+I)*i4^2*i3*i2*M+(2-2*I)*i4^2*i3*i1*M-(1+I)*i4^2*i2*i1*M

This isquadratic in M, and has two solutions in term of i1,i2,i3,i4 (which I won't show here because they're rather complicated).

> solve(%,M);

I'm sorry, but it doesn't make sense.  The solution of a differential equation depends on all the values of the functions on an interval, not only on the values at 4 points.

I'm sorry, but it doesn't make sense.  The solution of a differential equation depends on all the values of the functions on an interval, not only on the values at 4 points.

The equations can be obtained (in Firefox: I'm sure Internet Explorer is similar but maybe with slighly different names) by right-clicking on the images, choosing "Properties" and looking at "Alternate Text".

> Eq1:=55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*(2.048000000*10^15*a+7.680000000*10^12*b+2.560000000*10^10*c+6.40000000*10^7*d)/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-1.177600000*10^12*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-1.766400000*10^10*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0;
  Eq2:= 55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*(7.680000000*10^12*a+2.880000000*10^10*b+9.60000000*10^7*c+2.400000*10^5*d)/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-4.416000000*10^9*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-44160000*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0;
  Eq3:=  55200.0*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^.5*(2.560000000*10^10*a+9.600000000*10^7*b+320000*c+800*d)/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c))-14720000*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2*(4.800000000*10^5*a+1200*b+2*c))-73600*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c)^2)+18400*(1+d^2)^1.5/(c^2*d) = 0;
  Eq4:= 55200.0*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^.5*(6.400000000*10^7*a+2.400000000*10^5*b+800*c+2*d)/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c))-36800*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2*(4.800000000*10^5*a+1200*b+2*c))-55200.0*(1+d^2)^.5/c+18400*(1+d^2)^1.5/(c*d^2) = 0;
  sys:= convert({Eq1, Eq2, Eq3, Eq4}, rational); 
  S:= [solve(sys)];

This returns three solutions: two of them have a and b arbitrary, one has c arbitrary.  All involve complex numbers.  There are no real solutions.

The equations can be obtained (in Firefox: I'm sure Internet Explorer is similar but maybe with slighly different names) by right-clicking on the images, choosing "Properties" and looking at "Alternate Text".

> Eq1:=55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*(2.048000000*10^15*a+7.680000000*10^12*b+2.560000000*10^10*c+6.40000000*10^7*d)/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-1.177600000*10^12*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-1.766400000*10^10*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0;
  Eq2:= 55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*(7.680000000*10^12*a+2.880000000*10^10*b+9.60000000*10^7*c+2.400000*10^5*d)/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-4.416000000*10^9*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-44160000*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0;
  Eq3:=  55200.0*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^.5*(2.560000000*10^10*a+9.600000000*10^7*b+320000*c+800*d)/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c))-14720000*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2*(4.800000000*10^5*a+1200*b+2*c))-73600*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c)^2)+18400*(1+d^2)^1.5/(c^2*d) = 0;
  Eq4:= 55200.0*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^.5*(6.400000000*10^7*a+2.400000000*10^5*b+800*c+2*d)/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)*(4.800000000*10^5*a+1200*b+2*c))-36800*(1+(3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2)^1.5/((3.200000000*10^7*a+1.200000000*10^5*b+400*c+d)^2*(4.800000000*10^5*a+1200*b+2*c))-55200.0*(1+d^2)^.5/c+18400*(1+d^2)^1.5/(c*d^2) = 0;
  sys:= convert({Eq1, Eq2, Eq3, Eq4}, rational); 
  S:= [solve(sys)];

This returns three solutions: two of them have a and b arbitrary, one has c arbitrary.  All involve complex numbers.  There are no real solutions.

There are many ways to handle delay differential equations.    You might look at www.mapleprimes.com/forum/delaylogistic or groups.google.ca/group/sci.math.symbolic/browse_frm/thread/1f04240aa7b3a42c/fcbe7928da2ad9d1

There are many ways to handle delay differential equations.    You might look at www.mapleprimes.com/forum/delaylogistic or groups.google.ca/group/sci.math.symbolic/browse_frm/thread/1f04240aa7b3a42c/fcbe7928da2ad9d1

Another way is using a Groebner basis. 

> with(Groebner):
   G:= convert([seq(lhs(EQ||i),i=1..12)],rational);
   B:= Basis(G, plex(A1, A2, A3, B1, B2, B3, C2, D2, D3, P, a3, alfa));

In the (rather lengthy) result, the first element is a cubic polynomial in alfa which has roots
29516.12162, -3.424535+116.6084978*I, -3.424535-116.6084978*I.  So I'm quite confident that solve is not omitting any solutions.
 

Another way is using a Groebner basis. 

> with(Groebner):
   G:= convert([seq(lhs(EQ||i),i=1..12)],rational);
   B:= Basis(G, plex(A1, A2, A3, B1, B2, B3, C2, D2, D3, P, a3, alfa));

In the (rather lengthy) result, the first element is a cubic polynomial in alfa which has roots
29516.12162, -3.424535+116.6084978*I, -3.424535-116.6084978*I.  So I'm quite confident that solve is not omitting any solutions.
 

> map2({`=`}, m, %);

> map2({`=`}, m, %);

When using algebraic commands such as solve, it's sometimes better to use exact rationals rather than floats.

> sysr:= convert(sys,rational);
  solve(sysr);
  V:= [allvalues(%)];
  evalf[20](V);

[{A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = .18998177911873809335e-3, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = -.14623998857559515179, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = 29516.121625530464137, P = 304.45797935695207268, a3 = 1522.4817244768608202, alfa = 29516.121625530464137}, {A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = -.95966523985645933810e-4+.22776888089692050129e-4*I, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = .1855552024200438738e-2+.23441729108477519857e-1*I, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = -3.42453468052661868+116.60850395991041032*I, P = -153.79250638725309906+36.501423220660336746*I, a3 = -10.75314511305216855+122.12808630347007648*I, alfa = -3.42453468052661868+116.60850395991041032*I}, {A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = -.95966523985645933810e-4-.22776888089692050129e-4*I, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = .1855552024200438738e-2-.23441729108477519857e-1*I, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = -3.42453468052661868-116.60850395991041032*I, P = -153.79250638725309906-36.501423220660336746*I, a3 = -10.75314511305216855-122.12808630347007648*I, alfa = -3.42453468052661868-116.60850395991041032*I}]

There are three solutions, but only one of them is real.

> remove(has,%,I);

[{A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = .18998177911873809335e-3, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = -.14623998857559515179, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = 29516.121625530464137, P = 304.45797935695207268, a3 = 1522.4817244768608202, alfa = 29516.121625530464137}]

This is the only real solution.

When using algebraic commands such as solve, it's sometimes better to use exact rationals rather than floats.

> sysr:= convert(sys,rational);
  solve(sysr);
  V:= [allvalues(%)];
  evalf[20](V);

[{A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = .18998177911873809335e-3, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = -.14623998857559515179, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = 29516.121625530464137, P = 304.45797935695207268, a3 = 1522.4817244768608202, alfa = 29516.121625530464137}, {A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = -.95966523985645933810e-4+.22776888089692050129e-4*I, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = .1855552024200438738e-2+.23441729108477519857e-1*I, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = -3.42453468052661868+116.60850395991041032*I, P = -153.79250638725309906+36.501423220660336746*I, a3 = -10.75314511305216855+122.12808630347007648*I, alfa = -3.42453468052661868+116.60850395991041032*I}, {A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = -.95966523985645933810e-4-.22776888089692050129e-4*I, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = .1855552024200438738e-2-.23441729108477519857e-1*I, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = -3.42453468052661868-116.60850395991041032*I, P = -153.79250638725309906-36.501423220660336746*I, a3 = -10.75314511305216855-122.12808630347007648*I, alfa = -3.42453468052661868-116.60850395991041032*I}]

There are three solutions, but only one of them is real.

> remove(has,%,I);

[{A1 = -163.33387743488749552, A2 = -165.10415276719380192, A3 = .18998177911873809335e-3, B1 = 1.5300000000000000000, B2 = 17.718654803163617943, B3 = -.14623998857559515179, C2 = -8.5552387203254616659, D2 = 7.4813451968363820566, D3 = 29516.121625530464137, P = 304.45797935695207268, a3 = 1522.4817244768608202, alfa = 29516.121625530464137}]

This is the only real solution.

The requirement that f_C(t) >= 0 means that f_all(t) >= p*f_B(t) for all t, i.e. that p is bounded by the infimum of f_all(t)/f_B(t).  If, as in your case, f_all and f_B are normal with the same variance sigma^2 but different means mu_all and mu_B, that infimum is actually 0, because f_all(t)/f_B(t) = exp(((t - mu_B)^2 - (t - mu_all)^2)/(2*sigma^2))
= const * exp(t*(mu_all - mu_B)/sigma^2) -> 0 as t -> +infinity or -infinity, depending on the sign of mu_all - mu_B. 

The requirement that f_C(t) >= 0 means that f_all(t) >= p*f_B(t) for all t, i.e. that p is bounded by the infimum of f_all(t)/f_B(t).  If, as in your case, f_all and f_B are normal with the same variance sigma^2 but different means mu_all and mu_B, that infimum is actually 0, because f_all(t)/f_B(t) = exp(((t - mu_B)^2 - (t - mu_all)^2)/(2*sigma^2))
= const * exp(t*(mu_all - mu_B)/sigma^2) -> 0 as t -> +infinity or -infinity, depending on the sign of mu_all - mu_B.