## 13579 Reputation

13 years, 105 days

## If I fix...

ypur latest worksheet so that it actually executes, I produced the attached. I had to fix loads of syntax and make a lot of guesses aboout your intentions. Somee of these guesses may be wrong.

I also had to use the DirectSearch package from the Maple Application centre in order to solve the reulting algebraic equations. You might be able to get a simple fsolve() to work with these equations by restricting the solution intervals for all variables. I don't have the time/patience to try this

 > restart; with(plots);
 (1)
 > A := 0.5; M := 0.5; R := 0.5; Pr := 0.72; Ec := 0.1; g := 0.1;
 (2)
 > ODE := diff(f(eta), eta, eta, eta) + f(eta)*diff(f(eta), eta, eta) - diff(f(eta), eta)^2 - M*diff(f(eta), eta) - A*diff(f(eta), eta) - A*eta*diff(f(eta), eta, eta)/2 = 0, (1 + R)*diff(Theta(eta), eta, eta) + Pr*(f(eta)*diff(Theta(eta), eta) - Theta(eta)*diff(f(eta), eta)) - Pr*A*(Theta(eta) + eta/2*diff(Theta(eta), eta)) + Pr*(Ec*diff(f(eta), eta, eta)^2 + g*Theta(eta)); BCS:=f(0)=0,D(f)(0)=1,Theta(0)=1,D(f)(10)=0,D(Theta)(10)=0; sol:=dsolve({ODE,BCS},numeric,output = listprocedure): p1:=odeplot(sol,[eta,f(eta)],0..9,color=green,linestyle=3,thickness=3,legend = ["dsolve"]): p2:=odeplot(sol, [eta, Theta(eta)], 0 .. 9, color = green, linestyle = 3, thickness = 3, legend = ["dsolve"]):
 (3)
 > d1:= (f[i+1]-f[i-1])/(2*h): d2:=(f[i+1]-2*f[i]+f[i-1])/h^2: d3:= (f[i+1]-3*f[i]+3*f[i-1]-f[i-2])/h^3: d4:= (Theta[i+1]-Theta[i-1])/(2*h): d5:=(Theta[i+1]-2*Theta[i]+Theta[i-1])/h^2: h := 0.5: ode1 := d3 + f[i]*d2 - d1^2 - M*d1 - A*d1 - A*eta*d2/2 = 0: ode2 := (1 + R)*d5 + Pr*(-d1*Theta[i] + d4*f[i]) - Pr*A*(Theta[i] + eta/2*d4) + Pr*(Ec*d2^2 + g*Theta[i]):
 > eta[0]:=0: f[0]:=0: Theta[0] := 1: f[-1]:=f[1]-2*h: f[10]:=f[9]: Theta[10] := Theta[9]: for i from 0 to 9 do     eta[i+1]:=eta[i]+h;     Eq[i]:=eval(ode1,eta=eta[i]);     eq[i]:=eval(ode2,eta=eta[i]); end do:
 > soln1:= DirectSearch:-SolveEquations( [Eq[0],Eq[1],Eq[2],Eq[3],Eq[4],Eq[5],Eq[6],Eq[7],Eq[8],Eq[9]])[3];  soln2:= DirectSearch:-SolveEquations( eval([eq[0],eq[1],eq[2],eq[3],eq[4],eq[5],eq[6],eq[7],eq[8],eq[9]], soln1))[3];  data1:=[seq([eta[0]+i*h,eval(f[i], soln1 )],i=0..9)]:  data2 := [seq([h*i + eta[0], eval(Theta[i], soln2)], i = 0 .. 9)]:
 (4)
 > p3:=plot([data1],color=red,linestyle=1,thickness=3,legend = ["FDM"]):  p4:=plot([data2],color = red, linestyle = 1, thickness = 3, legend = ["FDM"]):  display({p1, p3}, axes = boxed);  display({p2, p4}, axes = boxed);

## Hmmm...

You state

i want to copare runge kutta 4th order method solution and FEM SOLUTION  comparison

Basically, you can't! Runge-Kutta methods can only be applied to initial value problems, and you have a boundary value problem. The Maple solution you provide in your latest worksheet, using the command

 msol:=dsolve({subs(beta=0.5,ODE),BCS},numeric,continuation=lambda1):

is actually computed using a finite difference method built-in to Maple. (I have no intention of debugging your attempt at writing your own finite difference method)

## And...

what is the relevance of this latest post? Seems to be nothing to do with your original question?

So precisely what are you trying to achieve?

## Don't understand this comment...

The" shooting mehod" is a way of converting a BVP to an IVP, by appropriate manipulation of the boundary conditions - in general it would not be used just to solve a BVP.

Suggest you read the help at

?dsolve/numeric/BVP

## Try...

adding the following to my previous answer

` parse~(Explode(Select(IsDigit, sstr1)));`

## Question...

makes no sense. Try posting a worksheet whicvh illustrates the problem

## There are...

a variety of ways - but I have to admit that I cheated a little, by using Maple's geometry() package to ectract relevant ellipse parameters, like the centre and lengths of the major and minor axes - see the attached

 > restart;   with(geometry):   TWeq:=2.96736996560705*10^(-12)*p^2+1.31319840299485*10^(-13)*t^2-8.89549693662593*10^(-7)* p+8.53128393394231*10^(-7)*t-3.65558815509970*10^(-30)*p*t-1 = 0:   ellipse( Ell, TWeq, [p,t]):   detail(Ell);   draw(Ell);
 >

## Well...

If you upload your worksheet using the big green up-arrow in the Mapleprimes edit toolbar, then someone can probably help you.

## After fixing...

numerous (minor) typos, the attached produces all the plots you *seem* to want.

PS I threw out a lot of code whichyou don't seem to need.

 >
 > restart: with(plots): rf := 1050: kf := .52: cpf := 3617: `&sigma:f` := .8: sigma1 := 4.10*10^7: rs1 := 19300: ks1 := 318: cps1 := 129: sigma2 := 10^(-10): rs2 := 3970: ks2 := 40: cps2 := 765: sigma3 := 6.30*10^7: rs3 := 10500: ks3 := 429: cps3 := 235: sigma4 := 10^(-12): rs4 := 4250: ks4 := 8.9538: cps4 := 686.2: l := 1: la := 1: Qc := 1: Tp := 1: g := 1: Ec := 1: H := 3: Rd := 1: M := 2: Pr := 21: p1 := 0.25e-1: p2 := 0.25e-1: t := (1/3)*Pi: lambda0 := 3: lambda1 := 2: delta1 := 1: a1 := 1/(1-p1-p2)^2.5: a1 := 1/(1-p1-p2)^2.5: a2 := 1-p1-p2+p1*rs1/rf+p2*rs2/rf: a3 := 1-p1-p2+p1*cps1/cpf+p2*cps2/cpf: a4 := ((ks1*p1+ks2*p2)/(p1+p2)+2*kf+2*(ks1*p1+ks2*p2)-2*(p1+p2)*kf)/((ks1*p1+ks2*p2)/(p1+p2)+2*kf-ks1*p1-ks2*p2+(p1+p2)*kf): a5 := 1+3*((p1*sigma1+p2*sigma2)/`&sigma:f`-p1-p2)/(2+(p1*sigma1+p2*sigma2)/((p1+p2)*`&sigma:f`)-((p1*sigma1+p2*sigma2)/`&sigma:f`-p1-p2)): a6 := 1/(1-p1-p2)^2.5: a7 := 1-p1-p2+p1*rs3/rf+p2*rs4/rf: a8 := 1-p1-p2+p1*cps3/cpf+p2*cps4/cpf: a9 := ((ks3*p1+ks4*p2)/(p1+p2)+2*kf+2*(ks3*p1+ks4*p2)-2*(p1+p2)*kf)/((ks3*p1+ks4*p2)/(p1+p2)+2*kf-ks3*p1-ks4*p2+(p1+p2)*kf): a10 := 1+3*((p1*sigma3+p2*sigma4)/`&sigma:f`-p1-p2)/(2+(p1*sigma3+p2*sigma4)/((p1+p2)*`&sigma:f`)-((p1*sigma3+p2*sigma4)/`&sigma:f`-p1-p2)): OdeSys := (a1+6*delta1*(diff(U0(Y), Y))^2)*(diff(U0(Y), Y, Y))-a2*R*(diff(U0(Y), Y))+lambda0-a5*M^2*U0(Y)+la*R*a2*Theta0(Y)*(1+Qc*Theta0(Y)), (a1+g*l+6*delta1*(diff(U0(Y), Y))^2)*(diff(U1(Y), Y, Y))+12*delta1*(diff(U0(Y), Y))*(diff(U1(Y), Y))*(diff(U0(Y), Y, Y))-a2*R*(diff(U1(Y), Y))-(H^2*a2*l+M^2*a5)*U1(Y)+lambda1+la*R*a2*Theta1(Y)*(1+2*Qc*Theta0(Y)), (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y)))*(diff(Theta0(Y), Y, Y))/Pr-a3*R*(diff(Theta0(Y), Y))+a5*Ec*M^2*U0(Y)^2+a1*Ec*(diff(U0(Y), Y))^2+2*Ec*delta1*(diff(U0(Y), Y))^4+4*Rd*(Tp-1)*(diff(Theta0(Y), Y))^2*(1+(2*(Tp-1))*Theta0(Y)+(Tp-1)^2*Theta0(Y)^2)/Pr, (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y)))*(diff(Theta1(Y), Y, Y))/Pr-a3*R*(diff(Theta1(Y), Y))-l*a3*H^2*Theta1(Y)+2*a5*Ec*M^2*U0(Y)*U1(Y)+2*a1*Ec*(diff(U0(Y), Y))*(diff(U1(Y), Y))+8*Ec*delta1*(diff(U0(Y), Y))^3*(diff(U1(Y), Y))+4*Rd*(Tp-1)*(Theta1(Y)+(2*(Tp-1))*Theta0(Y)*Theta1(Y)+(Tp-1)^2*Theta0(Y)^2*Theta1(Y))*(diff(Theta0(Y), Y, Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))^2*((Tp-1)^2*Theta0(Y)*Theta1(Y)+(Tp-1)*Theta1(Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))*(diff(Theta1(Y), Y))*(1+(2*(Tp-1))*Theta0(Y)+(Tp-1)^2*Theta0(Y)^2)/Pr:  Cond := U0(0) = 0, U1(0) = 0, Theta0(0) = 0, Theta1(0) = 0, U0(1) = 0, U1(1) = 0, Theta0(1) = 1, Theta1(1) = 0:  RVals := [1, 2, 3]: Ans1[j] := {}: Ans[j] := {}: for j to numelems(RVals) do       Ans[j] := dsolve(eval([OdeSys, Cond], R = RVals[j]), numeric, output = listprocedure) end do: OdeSys1 := (a6+6*delta1*(diff(V0(Y), Y))^2)*(diff(V0(Y), Y, Y))-a7*R*(diff(V0(Y), Y))+lambda0-a10*M^2*V0(Y)+la*R*a7*Phi0(Y)*(1+Qc*Phi0(Y)), (a6+g*l+6*delta1*(diff(V0(Y), Y))^2)*(diff(V1(Y), Y, Y))+12*delta1*(diff(V0(Y), Y))*(diff(V1(Y), Y))*(diff(V0(Y), Y, Y))-a7*R*(diff(V1(Y), Y))-(H^2*a7*l+M^2*a10)*V1(Y)+lambda1+la*R*a7*Phi1(Y)*(1+2*Qc*Phi0(Y)), (a9+(4*Rd*(1/3))*(1+(Tp-1)^3*Phi0(Y)^3+3*(Tp-1)^2*Phi0(Y)^2+(3*(Tp-1))*Phi0(Y)))*(diff(Phi0(Y), Y, Y))/Pr-a8*R*(diff(Phi0(Y), Y))+a10*Ec*M^2*V0(Y)^2+a6*Ec*(diff(V0(Y), Y))^2+2*Ec*delta1*(diff(V0(Y), Y))^4+4*Rd*(Tp-1)*(diff(Phi0(Y), Y))^2*(1+(2*(Tp-1))*Phi0(Y)+(Tp-1)^2*Phi0(Y)^2)/Pr, (a9+(4*Rd*(1/3))*(1+(Tp-1)^3*Phi0(Y)^3+3*(Tp-1)^2*Phi0(Y)^2+(3*(Tp-1))*Phi0(Y)))*(diff(Phi1(Y), Y, Y))/Pr-a8*R*(diff(Phi1(Y), Y))-l*a8*H^2*Phi1(Y)+2*a10*Ec*M^2*V0(Y)*V1(Y)+2*a6*Ec*(diff(V0(Y), Y))*(diff(V1(Y), Y))+8*Ec*delta1*(diff(V0(Y), Y))^3*(diff(V1(Y), Y))+4*Rd*(Tp-1)*(Phi1(Y)+(2*(Tp-1))*Phi0(Y)*Phi1(Y)+(Tp-1)^2*Phi0(Y)^2*Phi1(Y))*(diff(Phi0(Y), Y, Y))/Pr+8*Rd*(Tp-1)*(diff(Phi0(Y), Y))^2*((Tp-1)^2*Phi0(Y)*Phi1(Y)+(Tp-1)*Phi1(Y))/Pr+8*Rd*(Tp-1)*(diff(Phi0(Y), Y))*(diff(Phi1(Y), Y))*(1+(2*(Tp-1))*Phi0(Y)+(Tp-1)^2*Phi0(Y)^2)/Pr:  Cond1 := V0(0) = 0, V1(0) = 0, Phi0(0) = 0, Phi1(0) = 0, V0(1) = 0, V1(1) = 0, Phi0(1) = 1, Phi1(1) = 0: for j to numelems(RVals) do     Ans1[j] := dsolve(eval([OdeSys1, Cond1], R = RVals[j]), numeric, output = listprocedure)  end do:
 > cols := [red, blue, black]:   plotA:= display           ( [ seq               ( odeplot                 ( Ans[k],                   [Y, U0(Y)],                   0..1,                   color=cols[k]                  ),                  k=1..numelems(RVals)                )              ],              'axes'= 'boxed',              labels=[Y,'U0']             );   plotB:=  display            ( [ seq                ( odeplot                  ( Ans[k],                    [Y, U1(Y)],                    0..1,                    color=cols[k]                  ),                  k=1..numelems(RVals)                )              ],              'axes'= 'boxed',              labels=[Y,'U1']            );   plotC:= display           ( [ seq               ( odeplot                 ( Ans1[k],                   [Y, V0(Y)],                   Y=0..1,                   color=cols[k]                 ),                 k=1..numelems(RVals)               )             ],            'axes'= 'boxed',            'linestyle' = 'dashdot',            labels=[Y,'U0']          );   plotD:=  display            ( [ seq                ( odeplot                  ( Ans1[k],                    [Y,V1(Y)],                    Y=0..1,                    color=cols[k]                  ),                  k=1..numelems(RVals)                )              ],              'axes'= 'boxed',              'linestyle' = 'dashdot',              labels=[Y,'U1']            );   plotE:= display           ( [ seq               ( odeplot                 ( Ans[k],                   [Y,Theta0(Y)],                   Y=0..1,                   color=cols[k]                 ),                 k=1..numelems(RVals)               )             ],             'axes'= 'boxed',             labels=[Y,'Theta0']           );   plotF:=  display            ( [ seq                ( odeplot                  ( Ans[k],                    [Y, Theta1(Y)],                    Y=0..1,                    color=cols[k]                  ),                  k=1..numelems(RVals)                )              ],              'axes'= 'boxed',               labels=[Y,'Theta1']            );   plotG:= display           ( [ seq               ( odeplot                 ( Ans1[k],                   [Y, Phi0(Y)],                   Y=0..1,                   color=cols[k]                 ),                 k=1..numelems(RVals)               )            ],            'axes'= 'boxed',            'linestyle' = 'dashdot',            labels=[Y,'Theta0']          );   plotH:=  display            ( [ seq                ( odeplot                  ( Ans1[k],                    [Y,Phi1(Y)],                    Y=0..1,                    color=cols[k]                  ),                  k=1..numelems(RVals)                )              ],              'axes'= 'boxed',              'linestyle' = 'dashdot',              labels=[Y,'Theta1']            );
 > plots:-display([plotE, plotG]); plots:-display([plotF, plotH]);
 >

op(4, s2)[][];

## The function...

'g2' does not contain the name 'q', and so is entirely independent of it.

The expression TC contains the variable names

{Cd, D, O1, O2, P, beta, h1, h2, i, q, r, t1, t2, theta1, theta2}.

However when you differentiate this expression with respect to t2, the result contains only the variables

{Cd, D, O2, beta, h2, i, t2, theta2},

so

{O1, h1, q, t1, theta1}

have been eliminated

## fsolve()...

is a numerical solver. Thus every quantity in your equation(s), other than the variable(s) you want fsolve() to compute has to have a numeric value.

I notice that you sayyou are trying to  "solve a system of differential equations" for which I would have expected you to be using dsolve(), with, or without the numeric option.

This would be so much easier to diagnose/fix if you uploaded your worksheet using the big green up-arrow in the Mapleprimes toolbar

## So why don't you...

just use the built-in help which comes with your Maple release?

As in

Help->What's new

from the menus in the standard GUI

## As far as I can tell...

Maple 2017 produces "correct" answer. All subsequent Maple releases produce weird answer.

Using lists rather than expressions sequences appears to make the problem go away in later releases (Maple 2018 - Maple2023)

## Can reproduce in Maple2023...

This is a bug, and it is a bad one.

If one uses sequences of names rather than integers, then one can see the difference in the calculations being performed depending on whether the Units() package is loaded or not. The second calculation in the attached makes sense: the first one makes no sense at all.

 > restart;
 > interface(version); with(Units):