Question: How do I solve a system of higher order ODE in Mapple?

Hi, 

I am trying to solve a system of 3 ODE of order 4. 

Is there any template to do it? 

When I try it, I get a message:

Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 10, got 8
What is wrong with it?

> restart;
> with*plots; blt := 5;
> lambda := .1; delta := .1;
                                     0.1
                                     0.1
> gamma2 := .1;
                                     0.1
> lambda2 := .1; lambda[T] := .7; lambda[C] = .1;
> Pr := 2;
> Nb := 0.1e-1;
> Nt := 0.1e-1;
> Sc := 0.1e-1;
> M := .5; E = -.1;
> N := 0.1e-1; n = .1;
> sigma[1] := .1;
> alpha[0] := 0; omega := 0;
                                      0
                                      0
> Eq1 := diff(f(eta), eta, eta, eta)+gamma2*((diff(f(eta), eta, eta))*(diff(f(eta), eta, eta))-f(eta)*(diff(f(eta), eta, eta, eta, eta)))-(1+lambda2)*((diff(f(eta), eta))*(diff(f(eta), eta))-f(eta)*(diff(f(eta), eta, eta)))-(1+lambda2)*M*sin(omega)^2*(diff(f(eta), eta))+(1+lambda2)(lambda*(1+lambda[T]*theta(eta))*theta(eta)+lambda*N*(1+lambda[C]*phi(eta))*phi(eta))*cos(alpha[0]) = 0;
                                                           2
/  d   /  d   /  d         \\\       /  d   /  d         \\ 
|----- |----- |----- f(eta)||| + 0.1 |----- |----- f(eta)|| 
\ deta \ deta \ deta       ///       \ deta \ deta       // 

                                                                           2
                /  d   /  d   /  d   /  d         \\\\       /  d         \ 
   - 0.1 f(eta) |----- |----- |----- |----- f(eta)|||| - 1.1 |----- f(eta)| 
                \ deta \ deta \ deta \ deta       ////       \ deta       / 

                /  d   /  d         \\          
   + 1.1 f(eta) |----- |----- f(eta)|| + 1.1 = 0
                \ deta \ deta       //          
> Eq2 := 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))) = 0;
/  d   /  d             \\            /  d             \
|----- |----- theta(eta)|| + 2 f(eta) |----- theta(eta)|
\ deta \ deta           //            \ deta           /

                                                                       2    
          /  d             \ /  d           \        /  d             \     
   + 0.02 |----- theta(eta)| |----- phi(eta)| + 0.02 |----- theta(eta)|  = 0
          \ deta           / \ deta         /        \ deta           /     
> Eq3 := diff(phi(eta), eta, eta)+Sc*f(eta)*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta, eta))/Nb-Sc*sigma[1]*(1+delta*theta(eta))^n*exp(-E/(1+delta*theta(eta))) = 0;
       /  d   /  d           \\               /  d           \
       |----- |----- phi(eta)|| + 0.01 f(eta) |----- phi(eta)|
       \ deta \ deta         //               \ deta         /

                        /  d   /  d             \\
          + 1.000000000 |----- |----- theta(eta)||
                        \ deta \ deta           //

                                      n    /          E         \    
          - 0.001 (1 + 0.1 theta(eta))  exp|- ------------------| = 0
                                           \  1 + 0.1 theta(eta)/    
> bcs1 := f(0) = 0, (D(f))(0) = 1, (D(f))(blt) = 0, (D(D(f)))(blt) = 0, theta(0) = 1, theta(blt) = 0, Nb*(D(phi))(0)+Nt*(D(theta))(0) = 0, phi(blt) = 0;
   f(0) = 0, D(f)(0) = 1, D(f)(5) = 0, @@(D, 2)(f)(5) = 0, theta(0) = 1, 

     theta(5) = 0, 0.01 D(phi)(0) + 0.01 D(theta)(0) = 0, phi(5) = 0
> NULL;
> L := [.1, .5, 1.0, 5.0, 10.0];
                         [0.1, 0.5, 1.0, 5.0, 10.0]
> 
> for k to 5 do R := dsolve(eval({Eq1, Eq2, Eq3, bcs1}, M = L[k]), [f(eta), theta(eta), phi(eta)], numeric, method = bvp[midrich], maxmesh = 4096, abserr = 0.1e-1, output = listprocedure); X1 || k := rhs(R[3]); X2 || k := rhs(R[4]); X3 || k := rhs(R[5]); Y1 || k := rhs(R[6]); Y2 || k := -rhs(R[7]); Z1 || k := rhs(R[8]); Z2 || k := -rhs(R[9]) end do;
Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 10, got 8
> R;
                                      R

> print([(X1 || (1 .. 5))(0)]);
                  [X11(0), X12(0), X13(0), X14(0), X15(0)]
> print([(X2 || (1 .. 5))(0)]);
                  [X21(0), X22(0), X23(0), X24(0), X25(0)]
> 
> print([(Y1 || (1 .. 5))(0)]);
                  [Y11(0), Y12(0), Y13(0), Y14(0), Y15(0)]
> print([(Y2 || (1 .. 5))(0)]);
                  [Y21(0), Y22(0), Y23(0), Y24(0), Y25(0)]
> print([(Z1 || (1 .. 5))(0)]);
                  [Z11(0), Z12(0), Z13(0), Z14(0), Z15(0)]
> print([(Z2 || (1 .. 5))(0)]);
                  [Z21(0), Z22(0), Z23(0), Z24(0), Z25(0)]
> NULL;
> NULL;
> plot([X1 || (1 .. 5)], 0 .. blt, labels = [eta, (D(f))(eta)], thickness = 1, color = black);
Warning, unable to evaluate the functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

> plot([Y1 || (1 .. 5)], 0 .. blt, labels = [eta, theta(eta)], thickness = 1, color = black);
Warning, unable to evaluate the functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

> 
> plot([Z1 || (1 .. 5)], 0 .. blt, labels = [eta, phi(eta)], thickness = 1, color = black);
Warning, unable to evaluate the functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

>