Question: How to implement finite element method for the given equations and boundary conditions

  restart;

  local gamma:
  local GAMMA:

  odeSystem:= [ (1 + GAMMA)*diff(f(eta), eta$4) - S*(eta*diff(f(eta), eta$3) + 3*diff(f(eta), eta$2)
                +
                diff(f(eta), eta)*diff(f(eta), eta$2) - f(eta)*diff(f(eta), eta$3))
                -
                GAMMA*delta(2*diff(f(eta), eta$2)*diff(f(eta), eta$3)^2 + diff(f(eta), eta$2)^2*diff(f(eta), eta$4))
                -
                M^2*diff(f(eta), eta$2) = 0,

                (1 + (4*R)/3)*diff(theta(eta), eta$2) + Pr*S*(f(eta)*diff(theta(eta), eta)
                -
                eta*diff(theta(eta), eta) + Q*theta(eta)) = 0,

                diff(phi(eta), eta$2) + Sc*S*(f(eta)*diff(phi(eta), eta)
                -
                eta*diff(phi(eta), eta)) - Sc*gamma*phi(eta) = 0
              ];

  params:= [ S = 0.5, GAMMA = 0.1, delta = 0.1, gamma = 0.1, M = 1,
             Pr = 1, Ec = 0.2, Sc = 0.6, R = 1, Q = 1
           ];
  bcs :=[ f(0) = 0, (D@@2)(f)(0) = 0, f(1) = 1, D(f)(1) = 0, D(theta)(0) = 0,
          theta(1) = 1, phi(1) = 1, D(phi)(0) = 0
        ];