xb is undefined.
If you define it symbolically, then you run into the problem that dsolve() cannot handle two constraints on the same function. When dsolve is able to work at all, it resolves functions to the general form f(x)=expression in x + constant . In your case as you also have diff(f,x,x) then f' will be of the form f'(x) = expression in x + constant1 and f(x) will be of the form f(x) = expression in x + x * constant1 + constant0 . Those one-point boundary constants are what is being constrained by the initial conditions. After you resolve constant0 according to the boundary f(0)=whatever then dsolve() does not have any further manipulations available to try to solve a second constraint at the same level.
This does not stop you from leaving out the second boundary condition for the dsolve() and then substituting into the resulting functions to pin down any available freedoms.
The logical structure is like this:
Da = 10; Rd = 0.5; Tw = 1.5; C1 = 1; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr =7;
syms x xb f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x) f4(x) g4(x)
df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
eqn0 = [ d3f0 == 0, d2g0 == 0 ];
cond0a = [ f0(0) == 0, d2f0(0) == -2/C1, g0(0) == 1];
cond0b = [df0(xb) == 0, g0(xb) == 0];
F0 = dsolve(eqn0,cond0a);
cond0bsubs = subs(cond0b,{f0,g0},{F0.f0, F0.g0});
constant0_names = setdiff(symvar(cond0bsubs),xb);
constant1_vals = solve(cond0bsubs,constant0_names);
f0(x) = subs(F0.f0, constant1_vals);
g0(x) = subs(F0.g0, constant1_vals);
df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
df1 = diff(f1,x); d2f1 = diff(df1,x); d3f1 = diff(d2f1,x);
dg1 = diff(g1,x); d2g1 = diff(dg1,x);
eqn1 = [ a1*d3f1 + a2*(f0*d2f0 - df0^2 ) - Da*df0 == 0,
(1+Rd*A1*(g0*(Tw-1)+1)^3)*d2g1 + Pr*B1*( f0*dg0 - 2*df0*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*dg0) == 0];
cond1a = [ f1(0) == 0, d2f1(0) == 0, g1(0) == 0];
cond1b = [ df1(xb) == 0, g1(xb) == 0 ];
F1 = dsolve( eqn1, cond1a);
disp(F1.f1)
disp(F1.g1)
cond1bsubs = subs(cond1b, {f1,g1}, {F1.f1, F1.g1});
constant1_names = setdiff(symvar(cond1bsubs),xb);
constant1_vals = solve(cond1bsubs,constant1_names);
f1(x) = subs(F1.f1, constant1_vals);
g1(x) = subs(F1.g1, constant1_vals);
This gets through finding F1, but F1.g1 is given in a form that involves int() of a symsum() of an expression involving root() of a cubic. MATLAB does not offer any way to resolve the root() of the cubic, even though it would be possible to take it down to exact solutions.
It is a nuisance and takes a bunch of custom work to resolve the symsum(). I have my system attempting to resolve the int() at the moment, but it is slow.
In theory I could write a function to automatically resolve the symsum() in this case, but it would not be simple.
