Hi I've been trying to solve these set of PDEs below and have been at it for very long

My codes are below

I tried solving the first equation by using:

sys_ode := 2*(diff(T(eta), eta, eta, eta))+T(eta)*(diff(T(eta), eta, eta)) = 0

ics := T(0) = 0, (D(T))(0) = 1, (D(T))(20) = 0

Digits := 10

sol1 := dsolve({ics, sys_ode}, numeric, output = operator)

{q(eta) = rhs(sol1[2](eta)), w(eta) = rhs(sol1[3](eta))}

Then inputting into the subsequent two equations:

PDE1 := eval([2/P . (diff(g(x, eta), eta, eta))+q(eta)*(diff(g(x, eta), eta))-g(x, eta)*w(eta) = 2*x*w(eta)*(diff(g(x, eta), x)), 2/S . (diff(phi(x, eta), eta, eta))+q(eta)*(diff(phi(x, eta), eta)) = 2*x*w(eta)*(diff(phi(x, eta), x))])

subBC1 := -phi(x, 0)*exp(g(x, 0)*sqrt(x)/(1+varepsilon*g(x, 0)*sqrt(x)))

subBC2 := alpha . ((phi(x, 0))(sqrt(x)))(exp(g(x, 0)*sqrt(x)/(1+varepsilon*g(x, 0)*sqrt(x))))

BC := {diff(g(0, eta), eta, eta) = 0, diff(phi(0, eta), eta, eta) = 0, g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, (D[2](g))(x, 0) = subBC1, (D[2](phi))(x, 0) = subBC2}

P = 1

S = 1

pds := pdsolve(PDE, {BC}, numeric, spacestep = .25)

but always end up with :Error, (in pdsolve/numeric/process_PDEs) number of dependent variables and number of PDE must be the same

I know my BC conditions might probably have some major errors too but i really cant proceed on cos i always end up with this same error. I really hope anyone would be able to help me on this 

help please