Question: Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

  restart:
  ODES := (1+E)*(diff(f(eta), eta, eta, eta))+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2-E*delta*(diff(f(eta), eta, eta))^2*(diff(f(eta), eta, eta, eta))+M*(diff(f(eta), eta)-1)+Upsilon*(1-(diff(f(eta), eta))^2) = 0,
          diff(theta(eta), eta, eta)+(3/4)*R*(diff((1+(B-1)*theta(eta))^3*(diff(theta(eta), eta)), eta))+Pr*f(eta)*(diff(theta(eta), eta))+Pr*Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Pr*Nt*(diff(theta(eta), eta))*(diff(theta(eta), eta))+Pr*Df*(diff(phi(eta), eta, eta))+Pr*Ec*(diff(f(eta), eta, eta))^2+Q*theta(eta) = 0,
           diff(phi(eta), eta, eta)+Sc*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta, eta))/Nb-Sc*kappa*phi(eta)+Sc*Sr*(diff(theta(eta), eta, eta)) = 0:
  bcs:= f(0) = 0,
        D(f)(0) = 1,
        D(f)(10) = 0,
       
        
        (D(theta))(0) = -.2*(1-theta(0)),
        theta(10) = 0,
        phi(0) = 1,
        phi(10) = 0:
  params:=[E=EV,Nb = 3.9, lambda = .01, Sc = 1, Nt = .5, delta = 0.9, Pr = 7., M = mv, R =5., B = 6., Df = 1., Ec  =5., kappa = 1., Sr = 3., Upsilon =- 3., Q = 3.]:
  EVals:=[-1, -3, -7]:
  Mvals:= [3, 5, 7]:
  ans:=Matrix( numelems(Mvals)*numelems(EVals)+1, 5):
  ans[1,..]:= < 'M' | 'E' | expr1 | expr2 |expr3>:
  for k from 1 by 1 to 3 do
      mv:= Mvals[k]:
      for j from 1 by 1 to 3 do
          EV:=EVals[j]:
          sol:=dsolve( eval([ODES, bcs], params), numeric, output=listprocedure);
          ans[3*(k-1)+j+1,..]:= `<|>`([ mv,
                                        EV,
                                        R__e^(-0.5)*sh= eval( -diff(phi(eta), eta),
                                                              [sol[], params[]]
                                                            )(0),
                                        R__e^(-0.5)*NU= eval( -diff(theta(eta), eta),
                                                              [sol[], params[]]
                                                            )(0),
                                        R__e^(0.5)*C__f=eval( (1+0.5)*(diff(f(eta), eta,eta)-1/3*diff(f(eta), eta)^3),[sol[], params[]])(0)]);od:od:ans;
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging
               Matrix(%id = 18446744074195107958)
 for k from 1 by 1 to 3 do
      plot( [ seq( [ seq( [ans[j,1], rhs(ans[j,2+k]) ], j=i..10,3 ) ], i=2..4 ) ],
            color=[red, green, blue],
            labels=[typeset(M), typeset( lhs(ans[2,2+k]) )],
            labelfont=[times, bold, 20],
            legend=[typeset(E=0.01),typeset(E=0.1),typeset(E=0.2)],
            legendstyle=[font=[times, bold, 20]],
            title=typeset( ans[1,2+k], " versus ", M, " parameterized by ", E),
            titlefont=[times, bold, 24]
          )
 od;
Error, invalid input: rhs received 0, which is not valid for its 1st argument, expr