restart;
Digits := 5;
with(ODETools);
with(student);
with(plots);
inf := 5;
with(LinearAlgebra);
equ1 := (1+2*n)*f(eta)*(diff(theta(eta), eta))/(1+3*n)-(diff(theta(eta), `$`(eta, 2))) = 0;
equ2 := ((1+n)*(diff(f(eta), eta))^2/(1+3*n)-(1+2*n)*f(eta)*(diff(f(eta), eta, eta))/(1+3*n))/Bo+(diff(f(eta), `$`(eta, 3)))^n-theta(eta) = 0;
Bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(inf) = 0, theta(0) = 1, theta(inf) = 0;
Bo := 1; n := 2;
SolP1 := dsolve({Bcs, equ1, equ2}, numeric);
Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system
SolP1(0);
SolP1(0)

restart;
digit := 15;
F[0] := fw;
F[1] := delta*F[2]+epsilon;
F[2] := A;
T[0] := B;
T[1] := -1;
a := 0.01;
m := 1;
M := 1;
lambda := 1;
b := 1;

fw := 1;
epsilon := 1;
delta := 1;
Pr := 6;
Ec := 1;
R := 1;
for k from 0 to 7 do F[k+3] := (a*(sum((i+1)*(i+2)*F[i+2]*(k-i+1)*T[k-i+1], i = 0 .. k))-m(1+a*T[k]-(sum((i+1)*F[i+1]*(k-i+1)*F[k-i+1], i = 0 .. k))+a*(sum(sum(T[i]*(r-i+1)*F[r-i+1]*(k-r+1)*F[k-r+1], i = 0 .. r), r = 0 .. k)))-M*(1-(k+1)*F[k+1]+a*T[k]-a*(sum(T[i]*(k-i+1)*F[k-i+1], i = 0 .. k)))-((m+1)*(1/2))*(sum((i+1)*(i+2)*F[i+2]*F[k-i], i = 0 .. k)+a*(sum(sum(T[i]*F[k-r]*(r-i+1)*(r-i+2)*F[r-i+2], i = 0 .. r), r = 0 .. k)))-lambda*(T[k]+a*(sum(T[i]*T[k-i], i = 0 .. k))))/((k+1)*(k+2)*(k+3)); T[k+2] := (-(1+b/Pr)*(k+1)*T[k+1]-b*(sum(T[i]*(k-i+1)*T[k-i+1], i = 0 .. k))-M*Ec*(1-(2*(k+1))*(k+2)*F[k+2])-Ec*(1+M)*(sum((i+1)*(i+2)*F[i+2]*(k-i+1)*(k-i+2)*F[k-i+2], i = 0 .. k))-(1-2*m)*(sum(T[i]*(k-i+1)*F[k-i+1], i = 0 .. k))-((m+1)*(1/2))*(sum(F[i]*(k-i+1)*T[k-i+1], i = 0 .. k)))/((4/3)*R*(k+1)*(k+2)) end do;
FF := sum(F[j]*y^j, j = 0 .. 5);
TT := sum(T[j]*y^j, j = 0 .. 5);
with(numapprox); pade(FF, y, [4, 4]);
pade(diff(FF, y), y, [4, 4]);
pade(TT, y, [4, 4]);
solve({limit(pade(TT, y, [4, 4]), y = infinity) = 0, limit(pade(diff(FF, y), y, [4, 4]), y = infinity) = 1}, [A, B]);