Question: Solving ODEs first order and second order

Hello all,

I'm a new user, hopefully you can assist me as I have work on this problem for quite sometime, I keep getting error of Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system

I copy my work here. If anyone need the original worksheet I can send this immadiately. Please help and I would like to say Thank you.

R:= 8.31434;                             # R is the gas constant (J mol^-1 K^-1) 
F:= 96485;                               # F is the Faraday constant (A s mol^-1)   
alpha1:= 0.5;                            # alpha is the transfer coefficient (l) 
alpha2:= 0.5; 
T:= 298.15;                              # T is the absolute temperature (K)
j0H:= 0.87;                            # j0 is the exchange current density (A m^-2)
j0Cl2:= 0.00001;                         # j0 is the exchange current density (A m^-2)
  
                                   8.31434
                                    96485
                                     0.5
                                     0.5
                                   298.15
                                    0.87
                                   0.00001
Variables
Appliedpotential:= -1.0:  
d1:= 0.0085;                             # d1 is the thickness of cathode boundry layer (m)
d2:= 0.02;                              # d2 is the distance cathode boundry layer and membrane (m)   
d3:= 0.0018;                            # d3 is the thickness of the membrane (m)  
d4:= 0.02;                              # d4 is the distance between membrane and anode boundry layer(m) 
d5:= 0.00125;                             # d5 is the thickness of anode boundry layer(m)
n:=1;n2:= 2;z1:=1;z2:=2;                            # n is the number of electrons
alpha:= 0.5;
nu:= 10^(-5);                            # nu is the kinematic viscosity (m/s)
d:= 0.000015;                            # d is the diameter of graphite pores (m)      
v:= 0.002;                               # v is the volume (m^3)  
Ac:=0.0063;                               # Ac is the area of cathode (m^2)   
Am:=0.0064;                                # Am is area of the membrane (m^2)
Aa:= Ac;           
                                   0.0085
                                    0.02
                                   0.0018
                                    0.02
                                   0.00125
                                      1
                                      2
                                      1
                                      2
                                     0.5
                                     1  
                                   ------
                                   100000
                                  0.000015
                                    0.002
                                   0.0063
                                   0.0064
                                   0.0063
Concentrations
Bulk concentration of chloride, zinc(II) and nickel(II) ions in solution in mol/m3
cCl:= 5000;co2:=1e-3:
cZn2:= 30;
cNi2:= 18;pcl2:=1:ph2:=1:po2:=1;
                                    5000
                                     30
                                     18
                                      1
Equilibrium Electrode Potentials
ee1 (SHE / V) is the reversible potential for the reduction of zinc(II) an nickel(I) ions
eNi2 := evalf(-0.245+R*T/(2*F)*ln(10)*log10(abs(y1(x))));
eZn2 := evalf(-0.76+R*T/(2*F)*ln(10)*log10(abs(y2(x))));
eH:= evalf(-0.887-R*T/(2*F)*ln(10)*log10(ph2)+R*T/F*ln(10)*log10(y3(x)));
eCl2:= evalf(1.396+(0.059/2)*log10(0.001*(abs(y4(x))))-0.059*log10(0.006));
eo2:=1.229+R*T/(2*F)*ln(10)*log10(po2)+R*T/F*ln(10)*log10(abs(y3(x)));
ehocl:=1.494+R*T/(2*F)*ln(10)*log10(abs(C3(x)))+R*T/(2*F)*ln(10)*log10(abs(C5(x)))-0.0296*log10(abs(y4(x)));
                     -0.245 + 0.01284614433 ln(|y1(x)|)
                      -0.76 + 0.01284614433 ln(|y2(x)|)
                      -0.887 + 0.02569228866 ln(y3(x))
                1.527089076 + 0.01281168722 ln(0.001 |y4(x)|)
                      1.229 + 0.02569228866 ln(|y3(x)|)
        1.494 + 0.01284614433 ln(|C3(x)|) + 0.01284614433 ln(|C5(x)|)

             0.0296 ln(|y4(x)|)
           - ------------------
                   ln(10)     
Rate Coefficients
Diffusion Coefficient (Electrochemical Systems Third Edition  John Newman, Karen E Thomas-Alyea)
dCl:=2.03e-9:dCl2:=1.83e-9:dhocl:=1.28e-9:do2:=2e-9:dNi2:=5e-10:dZn2:=0.71e-9:doh:=5.26e-9:nu:=1e-6:dh:=9.312e-9:dNa:=1.334e-9:dso4:=1.065e-9:pH:=3;
                  
                                      3
Rate Coefficients
j0Cl:=1e-5:bcCl:=38:baCl:=77:bahocl:=39:bchocl:=39:bco:=19:bao:=38:bch:=19:bah:=19:baNi:19:bcNi:19:baZn:=50:bcZn:50:f:=25:
j0o2:=1e-8:j0hocl:=1e-6:
j0H:=0.87:j0Cl2:=1e-5:
j0Ni:=1e-4:j0Zn:=3.2e-7:
j(x):=j2(x)+j3(x)+j4(x)+j5(x):

Ionic Mobility,Ui
Error, missing operator or `;`
> mu[Ni2] := dNi2/(R*T); mu[Zn2] := dZn2/(R*T); mu[Cl2] := dCl2/(R*T); mu[H] := dh/(R*T); mu[Na] := dNa/(R*T); mu[Cl] := dCl/(R*T); mu[OH] := doh/(R*T); mu[SO4] := dso4/(R*T);
                                            -13
                              2.017007023 10  
                                            -13
                              2.864149973 10  
                                            -13
                              7.382245705 10  
                                            -12
                              3.756473879 10  
                                            -13
                              5.381374738 10  
                                            -13
                              8.189048516 10  
                                            -12
                              2.121891389 10  
                                            -13
                              4.296224961 10  
Diffusion Limiting Current Desnity
jL is the diffusion limiting current (A/m2)

jLCl2:=1.554*2*F*(dCl2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y4(x);
jLNi:=1.554*2*F*(dNi2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y1(x);
jLZn:=1.554*2*F*(dZn2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y2(x);
jLH:= 0.0193:
                                        (1/2)     
                          4486.505500 25      y4(x)
                                        (1/2)     
                          1889.096518 25      y1(x)
                                        (1/2)     
                          2386.600407 25      y2(x)
Charge Balance(Electroneutrality)
[H+] + 2* [Ni 2+] + 2* [Zn 2+] + [Na+] = [Cl-] + [OH-]+2*[SO42-]
eq1:=y3(x)+2*y1(x)+2*y2(x)+y7(x)-y4(x)-y5(x)-2*y6(x)=0;
       y3(x) + 2 y1(x) + 2 y2(x) + y7(x) - y4(x) - y5(x) - 2 y6(x) = 0
>
Distribution of electrode and electrolyte phase potential, depending on j(x)
Sigma[e] is the particle phase's conductivity and Sigma[s] is the electrolyte phase's conductivity (S/m)
Phi[e] is the electrode phase potential (V), Phi[s] is the electrolyte phase potential (V), x is the distance in the direction of current flow normalised with respect to L
kappa[s](x)=evalf(F^2*((abs(2^2)*mu[Ni2]*y1(x))+(abs(2^2)*mu[Zn2]*y2(x))+(abs(1^2)*mu[H]*y3(x))+(abs(1^2)*mu[Na]*y7(x))-(abs(-1^2)*mu[Cl]*y4(x))-(abs(-1^2)*mu[OH]*y5(x))-(abs(-2^2)*mu[SO4]*y6(x))));
eq1:=y3(x)+2*y1(x)+2*y2(x)+y7(x)-y4(x)-y5(x)-2*y6(x)=0;
kappa[e]:=7.407*10^5;

ODE1:=(1000*(1-lambda)+1)*diff(phi[e](x),x)=j(x)/(kappa[e]);

ODE2:=(1000*(1-lambda)+1)*diff(phi[s](x),x)=-j(x)/(F^2*((abs(2^2)*mu[Ni2]*y1(x))+(abs(2^2)*mu[Zn2]*y2(x))+(abs(1^2)*mu[H]*y3(x))+(abs(1^2)*mu[Na]*y7(x))-(abs(-1^2)*mu[Cl]*y4(x))-(abs(-1^2)*mu[OH]*y5(x))-(abs(-2^2)*mu[SO4]*y6(x))));

kappa[s](x) = 0.007510813947 y1(x) + 0.01066535580 y2(x) + 0.03497034973 y3(x)

   + 0.005009712903 y7(x) - 0.007623476159 y4(x) - 0.01975344069 y5(x)

   - 0.01599803371 y6(x)
       y3(x) + 2 y1(x) + 2 y2(x) + y7(x) - y4(x) - y5(x) - 2 y6(x) = 0
                                            5
                               7.40700000 10
                            / d           \                         
       (1001 - 1000 lambda) |--- phi[e](x)| = 0.000001350074254 j2(x)
                            \ dx          /                         

          + 0.000001350074254 j3(x) + 0.000001350074254 j4(x)

          + 0.000001350074254 j5(x)
                       / d           \     (j2(x) + j3(x) + j4(x) + j5(x))/
  (1001 - 1000 lambda) |--- phi[s](x)| = -                                
                       \ dx          /                                    

    /           /              -13                       -12     
    \9309355225 \8.068028092 10    y1(x) + 1.145659989 10    y2(x)

                     -12                       -13     
     + 3.756473879 10    y3(x) + 5.381374738 10    y7(x)

                     -13                       -12     
     - 8.189048516 10    y4(x) - 2.121891389 10    y5(x)

                     -12      \\
     - 1.718489984 10    y6(x)//

 

Relationship between the concentration c and the local current density
Eq2 is the H material balance, Eq3 is the Ni(II) material balance, Eq4 is the Zn(II) material balance, Eq5 is the Cl2 material balance
>
ODE3:=(1000*(1-lambda)+1)*v*n*F*((dh*diff(y3(x),x,x))+(z1*mu[H]*F*y3(x)*F*diff(phi[s](x),x,x))+(z1*mu[H]*F*diff(y3(x),x)*diff(phi[s](x),x)))=Ac*j2(x);

                                   /        -9 / d  / d       \\
      192.970 (1001 - 1000 lambda) |9.312 10   |--- |--- y3(x)||
                                   \           \ dx \ dx      //

                               / d  / d           \\
         + 0.03497034973 y3(x) |--- |--- phi[s](x)||
                               \ dx \ dx          //

                         -7 / d       \ / d           \\              
         + 3.624433822 10   |--- y3(x)| |--- phi[s](x)|| = 0.0063 j2(x)
                            \ dx      / \ dx          //              
ODE4:=(1000*(1-lambda)+1)*v*n2*F*((dNi2*diff(y1(x),x,x))+(z2*mu[Ni2]*y1(x)*F*diff(phi[s](x),x,x))+(z2*mu[Ni2]*F*diff(y1(x),x)*diff(phi[s](x),x)))=Ac*j3(x);  
         
                                   /    -10 / d  / d       \\
      385.940 (1001 - 1000 lambda) |5 10    |--- |--- y1(x)||
                                   \        \ dx \ dx      //

                         -8       / d  / d           \\
         + 3.892218452 10   y1(x) |--- |--- phi[s](x)||
                                  \ dx \ dx          //

                         -8 / d       \ / d           \\              
         + 3.892218452 10   |--- y1(x)| |--- phi[s](x)|| = 0.0063 j3(x)
                            \ dx      / \ dx          //              
ODE5:=(1000*(1-lambda)+1)*v*n2*F*((dZn2*diff(y2(x),x,x))+(z2*mu[Zn2]*y2(x)*F*diff(phi[s](x),x,x))+(z2*mu[Zn2]*F*diff(y2(x),x)*diff(phi[s](x),x)))=Ac*j4(x);   

                                   /      -10 / d  / d       \\
      385.940 (1001 - 1000 lambda) |7.1 10    |--- |--- y2(x)||
                                   \          \ dx \ dx      //

                         -8       / d  / d           \\
         + 5.526950203 10   y2(x) |--- |--- phi[s](x)||
                                  \ dx \ dx          //

                         -8 / d       \ / d           \\              
         + 5.526950203 10   |--- y2(x)| |--- phi[s](x)|| = 0.0063 j4(x)
                            \ dx      / \ dx          //              
ODE6:=(1000*(1-lambda)+1)*v*n2*F*((dCl2*diff(y4(x),x,x))+(z2*mu[Cl]*y4(x)*F*diff(phi[s](x),x,x))+(z2*mu[Cl]*F*diff(y4(x),x)*diff(phi[s](x),x)))=Ac*j5(x);

                                   /       -9 / d  / d       \\
      385.940 (1001 - 1000 lambda) |1.83 10   |--- |--- y4(x)||
                                   \          \ dx \ dx      //

                         -7       / d  / d           \\
         + 1.580240692 10   y4(x) |--- |--- phi[s](x)||
                                  \ dx \ dx          //

                         -7 / d       \ / d           \\              
         + 1.580240692 10   |--- y4(x)| |--- phi[s](x)|| = 0.0063 j5(x)
                            \ dx      / \ dx          //              
>
>
> ODE7 := (1000*(1-lambda)+1)*v*n*F*(doh*(diff(y5(x), x, x))+z1*mu[OH]*y5(x)*F*(diff(phi[s](x), x, x))+z1*mu[OH]*F*(diff(y5(x), x))*(diff(phi[s](x), x))) = 0;
                                       /       -9 / d  / d       \\
          192.970 (1001 - 1000 lambda) |5.26 10   |--- |--- y5(x)||
                                       \          \ dx \ dx      //

                             -7       / d  / d           \\
             + 2.047306907 10   y5(x) |--- |--- phi[s](x)||
                                      \ dx \ dx          //

                             -7 / d       \ / d           \\   
             + 2.047306907 10   |--- y5(x)| |--- phi[s](x)|| = 0
                                \ dx      / \ dx          //   
> ODE8 := (1000*(1-lambda)+1)*v*n*F*(dso4*(diff(y6(x), x, x))+z2*mu[SO4]*y6(x)*F*(diff(phi[s](x), x, x))+z2*mu[SO4]*F*(diff(y6(x), x))*(diff(phi[s](x), x))) = 0;
                                      /        -9 / d  / d       \\
         192.970 (1001 - 1000 lambda) |1.065 10   |--- |--- y6(x)||
                                      \           \ dx \ dx      //

                            -8       / d  / d           \\
            + 8.290425307 10   y6(x) |--- |--- phi[s](x)||
                                     \ dx \ dx          //

                            -8 / d       \ / d           \\   
            + 8.290425307 10   |--- y6(x)| |--- phi[s](x)|| = 0
                               \ dx      / \ dx          //   
> ODE9 := (1000*(1-lambda)+1)*v*n*F*(dNa*(diff(y7(x), x, x))+z1*mu[Na]*y7(x)*F*(diff(phi[s](x), x, x))+z1*mu[Na]*F*(diff(y7(x), x))*(diff(phi[s](x), x))) = 0;
                                      /        -9 / d  / d       \\
         192.970 (1001 - 1000 lambda) |1.334 10   |--- |--- y7(x)||
                                      \           \ dx \ dx      //

                            -8       / d  / d           \\
            + 5.192219416 10   y7(x) |--- |--- phi[s](x)||
                                     \ dx \ dx          //

                            -8 / d       \ / d           \\   
            + 5.192219416 10   |--- y7(x)| |--- phi[s](x)|| = 0
                               \ dx      / \ dx          //   

Butler-Volmer Relationship between reaction current density and overpotential
j2(x) is the current due to hydrogen evolution, j3(x) is the current due to Ni(II) reduction,  j4(x) is the current due to Zn(II) reduction, j5(x) is the current due to Cl2 reduction
j2(x):=-j0H*exp(-17.85*(phi[e](x)-phi[s](x)-eH));
j3(x):=(j0Ni*(-exp(-(phi[e](x)-phi[s](x)-eNi2)*alpha*F/(R*T))))/(1+j0Ni/jLNi*exp(-alpha*F*(phi[e](x)-phi[s](x)-eNi2)/(R*T)));j4(x):=(j0Zn*(-exp(-(phi[e](x)-phi[s](x)-eZn2)*alpha*F/(R*T))))/(1+j0Zn/jLZn*exp(-alpha*F*(phi[e](x)-phi[s](x)-eZn2)/(R*T)));
j5(x):=(j0Cl2*(-exp(-(phi[e](x)-phi[s](x)-eCl2)*alpha*F/(R*T))))/(1+j0Cl2/jLCl2*exp(-alpha*F*(phi[e](x)-phi[s](x)-eCl2)/(R*T)));


           -0.87 exp(-17.85 phi[e](x) + 17.85 phi[s](x) - 15.83295

              + 0.4586073526 ln(y3(x)))
 - (0.0001 exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 4.767967604

                                //      1   /              -9   (1/2)    
    + 0.2500000000 ln(|y1(x)|))) |1 + ----- \2.117414310 10   25      exp(
                                 \    y1(x)                              
 -19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 4.767967604

                               \\
    + 0.2500000000 ln(|y1(x)|))/|
                                /
  /      -7                                                                
- \3.2 10   exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 14.79043012

                              \//      1   /              -12   (1/2)    
   + 0.2500000000 ln(|y2(x)|))/ |1 + ----- \5.363277388 10    25      exp(
                                \    y2(x)                               
-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 14.79043012

                              \\
   + 0.2500000000 ln(|y2(x)|))/|
                               /
- (0.00001 exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) + 29.71882140

                                     //      1   /              -11   (1/2)
   + 0.2493294270 ln(0.001 |y4(x)|))) |1 + ----- \8.915624868 10    25     
                                      \    y4(x)                           

  exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) + 29.71882140

                                    \\
   + 0.2493294270 ln(0.001 |y4(x)|))/|
                                     /
Numerical solution of the differential equation system using Maple's dsolve routine
> dsys := {eq1, ODE1, ODE2, ODE3, ODE4, ODE5, ODE6, ODE7, ODE8, ODE9, y1(d1) = cNi2, y2(d1) = cZn2, y3(d1) = 5000, y4(d1) = 5000, y5(d1) = 96, y6(d1) = 1000, y7(d1) = 2000, phi[e](0) = Appliedpotential, phi[s](0) = 0};

> sol10 := dsolve(dsys, numeric, continuation = lambda, abserr = 0.1e-4, maxmesh = 8192);

Please Wait...