The complex values must come from eq3, where you have (1 + delta*theta(eta))^n and where later n = 0.2.
(-2)^0.2; # answer: 0.9293164906 + 0.6751879524*I (the principal 5th root.).
You could then think of using surd (see the help). That turns out not to work in dsolve/numeric/bvp because the round disappears and you get surd( ... , 5. ), i.e a float instead of an integer as the second argument.
You can use eq3_new instead of eq3:
eq3_new:=subs((1 + delta*theta(eta))^n=abs(1 + delta*theta(eta))^n*signum(1 + delta*theta(eta)),eq3);
The problem is illustrated here:
surd(-2, round(1/0.2)); #Here it works because the integer 5 remains an integer
### The alternative
You run into singularity problems though, so there are obviously other problems with your system.
With these extra arguments to dsolve
you get the error:
Error, (in dsolve/numeric/bvp) matrix is singular
The exp term in eq3 gives problems here:
exp(-E/(1 + delta*theta(eta)));
plot(exp(-10/(1 + 10*theta)),theta=-0.5..0.5,thickness=3);