NDSolve second-order differential eqns with boundary condition

Questions
,
1

Change the yo (3rd boundary condition: y[0]) to meet the target(4th boundary condition): funy’[1] + pey*bd*funy[1]
=0, and get the answer yo (0<yo<1)
Initial yo=0.3
but,the program cannot run, please give me some help, many thans------------------------------------Clear[fx, fy, cxo, cy1, h, pey, noofc, hox, ans];

(*Input data for the calculation of the outlet concentrations*)
{
{Flowrate of Cu Raff (m3/h), fx = 5}
} ;
{
{Flowrate of barren solvent (m3/hr), fy = 1}
} ;
{
{Subscript[U, 3] Subscript[O, 8] concentration in Cu Raff (mg/L),
cxo = 267}
} ;
{
{Subscript[U, 3] Subscript[O, 8]
     concentration in barren solvent (mg/L), cy1 = 130}
} ;
{
{Effective column height (m), h = 9}
} ;
{
{Peclet number of organic phase , pey = .7}
} ;
{
{no of compartments, noofc = 13}
} ;
{
{Height of transfer unit (m), hox = 1.65855537601783531}
} ;

(* Equilibrium constants represented by Langmuir model, \
cx^*=(a*cy)/(b*cy+1) *)
a = 0.00982089325801922;
b = -0.0002962283775423278;

(* Press once you complete the input section *)

nox = h/hox;
uxuy = fx/fy; dc = h/noofc; bd =
h/dc; SuperStar = (a/cxo - b)^-1; SuperStar = (a*cy1)/(
b*cy1 + 1);
b1 = cxo - SuperStar; b2 = SuperStar - cy1; a1 =
a*(SuperStar - cy1); a2 = a*cy1; a3 = b*(SuperStar - cy1);
a4 = (b*cy1) + 1;
yo = .3; parta :=
Module[{useless},
ini = NDSolve[{-x’ - (nox/
         b1)*(b1*x + SuperStar - (a1*y + a2)/(
          a3*y + a4)) == 0,
     y’’ +
       pey*bd*y’ + ((nox*pey*bd*uxuy)/
         b2)*(b1*x + SuperStar - (a1*y + a2)/(
          a3*y + a4)) == 0, y’[0] == 0, x[0] == 1,
     y[0] == yo}, {x, y}, {z, 0, 1}]; funx = ini[[1, 1, 2]];
funy = ini[[1, 2, 2]]; funy’[1] + pey*bd*funy[1]];
try = FindRoot[parta == 0, {ans, 0, 1}, MaxIterations → 350,
AccuracyGoal → 6];
Plot[{funx, funy}, {z, 0, 1}, Frame → True,
FrameLabel → {"Dimensionless length (Z) ",
“Dimensionless conc (X,Y)”}, PlotStyle → AbsoluteThickness[2],
PlotRange → All, ImageSize → 350]
Print["y0 = " , ans, " mg/L "]
Print["Concentration in Raff = " ,
funx[1]*(cxo - SuperStar) + SuperStar, " mg/L "]
Print["Concentration in Loaded Solvent = " ,
funy[0]*(SuperStar - cy1) + cy1, " mg/L "]The 1st and 2nd boundary conditions are:y’[0] == 0, x[0] == 1

2

Many parts of this worry me.

In[1]:= Clear[fx, fy, cxo, cy1, h, pey, noofc, hox, ans];
(*Input data for the calculation of the outlet concentrations*)
(*Flowrate of Cu Raff (m3/h)*)fx = 5;
(*Flowrate of barren solvent (m3/hr)*)fy = 1;
(*Subscript[U,3] Subscript[O,8] concentration in Cu Raff (mg/L)*)cxo = 267;
(*Subscript[U,3] Subscript[O,8] concentration in barren solvent(mg/L)*)cy1 = 130;
(*Effective column height (m)*)h = 9;
(*Peclet number of organic phase*)pey = .7;
(*no of compartments*)noofc = 13;
(*Height of transfer unit (m)*)hox = 1.65855537601783531;
(*Equilibrium constants represented by Langmuir model,cx^*=(a*cy)/(b*cy+1)*)
a = 0.00982089325801922;
b = -0.0002962283775423278;
(*Press once you complete the input section*)
nox = h/hox;
uxuy = fx/fy; dc = h/noofc; bd = h/dc; SuperStar = (a/cxo - b)^-1; SuperStarcx1 = (a*cy1)/(b*cy1 + 1);
b1 = cxo - SuperStarcx1; b2 = SuperStar - cy1; a1 = a*(SuperStar - cy1); a2 = a*cy1; a3 = b*(SuperStar - cy1);
a4 = (b*cy1) + 1;
(*yo=.3;*)
parta[yo_?NumericQ] := Module[{}, 
  ini = NDSolve[{-x'[z] - (nox/b1)*(b1*x[z] + SuperStarcx1 - (a1*y[z] + a2)/(a3*y[z] + a4)) == 0, 
     y''[z]+pey*bd*y'[z]+((nox*pey*bd*uxuy)/b2)*(b1*x[z]+SuperStarcx1-(a1*y[z]+a2)/(a3*y[z]+a4)) == 0, 
     y'[0] == 0, x[0] == 1, y[0] == yo}, {x[z], y[z]}, {z, 0, 1}];
  funx[z_] = ini[[1, 1, 2]];
  funy[z_] = ini[[1, 2, 2]];
  funy'[1] + pey*bd*funy[1]
  ]

In[17]:= lowz = 0.; highz = 1.;
Do[midz = (lowz + highz)/2.0;
  If[Sign[parta[lowz]] == Sign[parta[midz]], lowz = midz, highz = midz], {30}
  ];
ans = midz;

(*try=FindRoot[parta\[Equal]0,{ans,0,1},MaxIterations\[Rule]350,AccuracyGoal\[Rule]6];*)

In[20]:= Plot[{funx[z], funy[z]}, {z, 0, 1}, Frame -> True, 
 FrameLabel -> {"Dimensionless length (Z) ", "Dimensionless conc (X,Y)"},
 PlotStyle -> AbsoluteThickness[2], PlotRange -> All, ImageSize -> 350]

Out[20]= ...PlotSnipped...

In[21]:= Print["y0 = ", ans, " mg/L "]

During evaluation of In[21]:= y0 = 0.458824 mg/L 

In[22]:= Print["Concentration in Raff = ", funx[1]*(cxo - SuperStarcx1) + SuperStarcx1, " mg/L "]

During evaluation of In[22]:= Concentration in Raff = 3.36841 mg/L 

In[23]:= Print["Concentration in Loaded Solvent = ", funy[0]*(SuperStar - cy1) + cy1, " mg/L "]

During evaluation of In[23]:= Concentration in Loaded Solvent = 1448.16 mg/L 

The 1 st and 2 nd boundary conditions are : y'[0] == 0, x[0] == 1

Please test this carefully to find and correct more errors and misunderstandings