Question: How do I solve system of differential equations in finite difference method?

How do  I solve system of differential equations in finite difference method or finite element method?

eq1 := (diff(f(x), x, x, x))*(a*beta*f(x)^2-1)+(diff(f(x), x))^2-2*a*beta*f(x)*(diff(f(x), x))*(diff(f(x), x, x))+(diff(f(x), x))*(M+k[1])-(diff(f(x), x, x))*f(x)-(alpha*theta(x)+delta*phi(x))/rho = 0;

eq2 := -(diff(theta(x), x, x))*K[SB]*(Df-(Rd+k[hnf]/k[bf])/Pr)+N[t]*K[SB]*(diff(theta(x), x))^2-N[b]*(diff(theta(x), x))*(diff(phi(x), x))-(diff(f(x), x))*(diff(theta(x), x))-lambda*theta(x)-mu*Ec*(M*(diff(f(x), x))^2+(diff(f(x), x, x))^2) = 0;

eq3 := diff(phi(x), x, x)+Le*Sr*(diff(theta(x), x, x))+Le*f(x)*(diff(phi(x), x)) = 0;

ics := f(0) = 0, (D(f))(0) = 0, theta(0) = 1, phi(0) = 1;

bcs := (D(f))(100) = 0, theta(100) = 0, phi(100) = 0;

Parameters1 := rho = 2063.905, k[hnf] = .29942, k[bf] = .2520, mu = .38694, a = .1, beta = 5, k[1] = 2.0, M = 10, alpha = 20, delta = 20, K[SB] = .5, Df = 3, Pr = 1.2, Rd = 5, N[t] = 1.2, N[b] = 1.0, lambda = 1.5, Ec = 5, Le = .1, Sr = .1;