Currently, when I have solution to an ode, say y(x)=sin(x)+_C1 and have some initial condition, then to solve for _C1,  I manually substitute the solution into the IC and replace each x by x0 and replace each derivative manually and so on.

This is because I could not find automatic way to do this. Using another software, it is possible to automate this by writing the solution using the  y -> Function[{x}, ...] syntax. But in Maple, I was not sure how to do the same.

Here is a simple made up example. 

sol:= y(x) = sin(x)+_C1;
IC := a*D(y)(x0)+c*y(x0)= b*y0+exp((D@@2)(y)(x0));

The goal is to replace the solution (which is function y(x)) into the IC, and automatically replace all its derivatives and replace x by x0 then solve for _C1 from the equation that results.

Now I do this manually like this

eval(IC,[ y(x0)=eval(rhs(sol),x=x0), 
          D(y)(x0)=eval(diff(rhs(sol),x),x=x0), 
         (D@@2)(y)(x0)= eval(diff(rhs(sol),x$2),x=x0) ])

which gives

But this is too much work.

Using the other software, I can do the above much more easily like this

sol = y -> Function[{x}, Sin[x] + C[1]]
ic = a*y'[x0] + c*y[x0] == b*y0 + Exp[y''[x0]];
ic /. sol

I looked at algsubs, dchange, or making the solution as function instead, and so on but could not emulate the y -> Function[{x}, Sin[x] + C[1]] method in Maple.

What would be similar method in Maple to do the above automatically?  May be there is already builtin function in Maple?

Variation of (a) Skin friction ∂W/∂Z​, (b) Heat Transfer ∂θ/∂Z​, and (c) Mass Transfer ∂ϕ/∂Z​ for γ=10.0, Pr=7.0, ε=1.0, Nt=0.4, and Nb=0.2.

how to get this values by solving the PDE by using the pdsolve method 

 # Equation 37 (continuity): Simplified

    pde1 := diff(U(X,Z), X) + diff(W(X,Z), Z) = 0;

    # Equation 38 (momentum): Simplified

    pde2 := U(X,Z)*diff(U(X,Z), X) + W(X,Z)*diff(U(X,Z), Z) = 

           diff(U(X,Z), Z, Z) - theta(X,Z) - phi(X,Z);

    # Equation 39 (energy): Simplified

    pde3 := U(X,Z)*diff(theta(X,Z), X) + W(X,Z)*diff(theta(X,Z), Z) = 

           (1/Pr)*diff(theta(X,Z), Z, Z) + 

           Nb*diff(phi(X,Z), Z)*diff(theta(X,Z), Z) + 

           Nt*diff(theta(X,Z), Z)^2;

    # Equation 40 (concentration): Simplified

    pde4 := U(X,Z)*diff(phi(X,Z), X) + W(X,Z)*diff(phi(X,Z), Z) = 

           (1/Sc_val)*diff(phi(X,Z), Z, Z) + 

           (Nt/Nb)*diff(theta(X,Z), Z, Z);

    # Boundary conditions

    bcs := {

        # At Z = 0

        U(X,0) = 0,

        W(X,0) = 0,

        theta(X,0) = 1,

        phi(X,0) = 1,

        # At Z = Z_max

        U(X,Z_max) = 0,

        theta(X,Z_max) = 0,

        phi(X,Z_max) = 0,

    };