Generalization needed in dsolve code

Pr = 1;
        ODE = @(x,y) [y(2);  y(3);  -(1/2)*y(3)*y(1);  y(5);  - (Pr/2)*y(1)*y(5);];
        BC = @(ya,yb)[ya(1);ya(2);ya(4)-1;yb(2)-1; yb(4);];
        xa = 0; xb = 5; xn = linspace(xa,xb,100);    x = xn;
        solinit = bvpinit(x,[0 1 0 1 0]);        sol = bvp5c(ODE,BC,solinit);  S = deval(sol,x);
        fV = deval(sol,0);    p1 = fV(3); q1 = fV(5); 
        
Pr = sym('Pr');x = sym('x');f0(x) = sym('f0(x)');  g0(x) = sym('g0(x)');  
eqn0 = [ diff(f0,3) == 0, diff(g0,2) == 0];
cond0 = [f0(0) == 0, subs(diff(f0),0) == 0, subs(diff(f0),5) == 1, g0(0) == 1, g0(5) == 0];
F0 = dsolve(eqn0,cond0);     f0 = F0.f0; g0 = F0.g0;   % disp([f0,g0])

f1(x) = sym('f1(x)');  g1(x) = sym('g1(x)');  
eqn1 = [ diff(f1,3) + (1/2)*f0*diff(f0,2) == 0, diff(g1,2) + (Pr/2)*f0*diff(g0) == 0];
cond1 = [f1(0) == 0, subs(diff(f1),0) == 0, subs(diff(f1),5) == 0, g1(0) == 0, g1(5) == 0];
F1 = dsolve(eqn1,cond1);     f1 = F1.f1; g1 = F1.g1;   %disp([f1,g1])

f2(x) = sym('f2(x)');  g2(x) = sym('g2(x)');  
eqn2 = [ diff(f2,3) + (1/2)*(f0*diff(f1,2) + f1*diff(f0,2)) == 0, diff(g2,2) + (Pr/2)*(f0*diff(g1) + f1*diff(g0)) == 0];
cond2 = [f2(0) == 0, subs(diff(f2),0) == 0, subs(diff(f2),5) == 0, g2(0) == 0, g2(5) == 0];
F2 = dsolve(eqn2,cond2);     f2 = F2.f2; g2 = F2.g2;   %disp([f2,g2])

f3(x) = sym('f3(x)');  g3(x) = sym('g3(x)');  
eqn2 = [ diff(f3,3) + (1/2)*(f0*diff(f2,2) + f1*diff(f1,2) + f2*diff(f0,2)) == 0, diff(g3,2) + (Pr/2)*(f0*diff(g2) + f1*diff(g1) + f2*diff(g0)) == 0];
cond2 = [f3(0) == 0, subs(diff(f3),0) == 0, subs(diff(f3),5) == 0, g3(0) == 0, g3(5) == 0];
F3 = dsolve(eqn2,cond2);     f3 = F3.f3; g3 = F3.g3;   %disp([f3,g3])

f = f0 + f1 + f2 + f3;    f = collect(f,x);
g = g0 + g1 + g2 + g3;    g = collect(g,x);  g = subs(g,Pr,1);

figure(2),plot(xn,S(2,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bff^{\prime}(\eta)');hold on,fplot(diff(f),[0 5],'--','LineWidth',1.5)
   figure(4),plot(xn,S(4,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bf\theta(\eta)');hold on,fplot(g,[0 5],'--','LineWidth',1.5)
%%  I need to merge all the dsolve code into a single code to solve the problem given in attached pdf,
%%  NUMERICAL code is of Eqns (12) &  (14) but  dsolve  code is of Eqns (17) - (24).  
%% Any attempt will be a great work

[MINATI]

After submitting this code in this forum, I didn't get any Email.Why?

[Walter Roberson]

What did you expect the email to say?

The activity notices are only sent when someone else responds to you or edits your postings.

[MINATI PATRA]

ooh

Anyway,

Dear Walter Roberson

any idea if you can impose

[Walter Roberson]

I need to merge all the dsolve code into a single code

Okay, so why not

F = dsolve([eqn1, eqn2, eqn3, eqn4, cond1, cond2, cond3, cond4])

where you have been careful to give different variable names for variables that are not intended to be related ?

[MINATI PATRA]

If you can follow the PDF and Eqns. mentioned above, may be a loop or recurrence format can be created to solve.

N.B: First dsolve code and rest dsolve codes have different conditions:

[Walter Roberson]

Yes, perhaps a loop or recurrence formula could be created to solve, but you do not need those.

You have a fixed number of functions, f0, f1, f2, f3 . You have a fixed number of derivatives you are taking of them,

df0 = diff(f0,theta)
d2f0 = diff(df0, theta)
d3f0 = diff(d2f0, theta);

and so on for the other functions.

Then you can construct whatever equations are desired, such as

eqn1 = d3f0 + f0 = theta*(t-1)   %substitute actual equation here
cond01 = f0(0) == 1;
cond02 = df0(0) == 0;
cond03 = d2f0(0) == 0;
cond04 = d3f0(0) == 0;

and so on for the other equations.

Then

equations = [eqn1, eqn2, eqn3, eqn4, ...];
conds = [cond01, cond02, cond03, cond04, cond11, cond12, ...];

solution = solve(equations, conds);

[MINATI PATRA]

No Dear Walter Roberson, Actually I have posted a small part of my work. I need to find

f = f0 + f1 + ... + f30; or some thing like this type, which can't be possible to code explicitly. Thats why I need a recurrence type dsolve code to solve this.

Please help.

Thanks

[Walter Roberson]

Here is the technical trick you need:

num_f = 31;
syms x
funs = arrayfun(@(idx) str2sym(sprintf('f%d(x)', idx)), 0:num_f-1);
d1 = arrayfun(@(f)diff(f,x), funs);
d2 = arrayfun(@(f)diff(f,x), d1);
d3 = arrayfun(@(f)diff(f,x), d2);

Now you can proceed with other arrayfun that create your equations in vectorized form, like sum(funs) for f0 + f1 + ... f30

[MINATI PATRA]

@Dear Walter, how to include B.Cs (It is slightly different for f0 and other steps)

%% Please arrange this to run

Pr = 1;

ODE = @(x,y) [y(2); y(3); -(1/2)*y(3)*y(1); y(5); - (Pr/2)*y(1)*y(5);];

BC = @(ya,yb)[ya(1);ya(2);ya(4)-1;yb(2)-1; yb(4);];

xa = 0; xb = 5; xn = linspace(xa,xb,100); x = xn;

solinit = bvpinit(x,[0 1 0 1 0]); sol = bvp5c(ODE,BC,solinit); S = deval(sol,x);

fV = deval(sol,0); p1 = fV(3); q1 = fV(5);

Pr = sym('Pr');x = sym('x');f0(x) = sym('f0(x)'); g0(x) = sym('g0(x)');

eqn0 = [ diff(f0,3) == 0, diff(g0,2) == 0];

cond0 = [f0(0) == 0, subs(diff(f0),0) == 0, subs(diff(f0),5) == 1, g0(0) == 1, g0(5) == 0];

F0 = dsolve(eqn0,cond0); f0 = F0.f0; g0 = F0.g0; % disp([f0,g0])

num_f = 31; num_g = 31; syms x

funs = arrayfun(@(idx) str2sym(sprintf('f%d(x)', idx)), 0:num_f-1);

df1 = arrayfun(@(f)diff(f,x), funs); df2 = arrayfun(@(f)diff(f,x), df1); df3 = arrayfun(@(f)diff(f,x), df2);

guns = arrayfun(@(idx) str2sym(sprintf('g%d(x)', idx)), 0:num_g-1);

dg1 = arrayfun(@(g)diff(g,x), guns); dg2 = arrayfun(@(g)diff(g,x), dg1); dg3 = arrayfun(@(g)diff(g,x), dg2);

f = sum(funs); g = sum(guns);

figure(1),plot(xn,S(2,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bff^{\prime}(\eta)');hold on,fplot(diff(f),[0 5],'--','LineWidth',1.5)

figure(2),plot(xn,S(4,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bf\theta(\eta)');hold on,fplot(g,[0 5],'--','LineWidth',1.5)