# Question:unable to plot 3d and 2d plots

## Question:unable to plot 3d and 2d plots

Maple 2018

 > restart:
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >

 >
 > ODE:=[(a2+K)*(diff(U0(eta), eta, eta))/a1-Ra*(diff(U0(eta), eta))+lambda0/a1-a5*M1^2*U0(eta)/a1+K*(diff(N0(eta), eta))/a1+la*Ra*Theta0(eta)*(1+Qc*Theta0(eta)), (a2+K)*(diff(U1(eta), eta, eta))/a1-H^2*l1*U1(eta)-Ra*(diff(U1(eta), eta))+lambda1/a1-a5*M1^2*U1(eta)/a1+K*(diff(N1(eta), eta))/a1+la*Ra*(Theta1(eta))(1+2*Qc*Theta0(eta)), diff(N0(eta), eta, eta)-Ra*a1*Pj*(diff(N0(eta), eta))-2*n1*N0(eta)-n1*(diff(U0(eta), eta)), diff(N1(eta), eta, eta)-Ra*a1*Pj*(diff(N1(eta), eta))-2*n1*N1(eta)-n1*(diff(U1(eta), eta))-H^2*a1*Pj*l1*N1(eta), (a4/(a3*Pr)-delta*Ra^2/H^2+4*Rd*(1+(Tp-1)^3*Theta0(eta)^3+3*(Tp-1)^2*Theta0(eta)^2+(3*(Tp-1))*Theta0(eta))/(3*a3*Pr))*(diff(Theta0(eta), eta, eta))-Ra*(diff(Theta0(eta), eta))+a5*Ec*M1^2*U0(eta)^2/a3+(a2+K)*Ec*(diff(U0(eta), eta))^2/a1+Q*Theta0(eta)/a3+4*(diff(Theta0(eta), eta))^2*Rd*(3*(Tp-1)+6*(Tp-1)^2*Theta0(eta)+3*(Tp-1)^3*Theta0(eta)^2)/(3*a3*Pr), (a4/(a3*Pr)-delta*Ra^2/H^2+4*Rd*(1+(Tp-1)^3*Theta0(eta)^3+3*(Tp-1)^2*Theta0(eta)^2+(3*(Tp-1))*Theta0(eta))/(3*a3*Pr))*(diff(Theta1(eta), eta, eta))-(H^2*l1+2*Ra*delta*l1+Ra)*(diff(Theta1(eta), eta))+(Q/a3-delta*H^2*l1^2)*Theta1(eta)+2*(a2+K)*Ec*(diff(U0(eta), eta))*(diff(U1(eta), eta))/a1+2*a5*Ec*M^2*U0(eta)*U1(eta)/a3+4*(diff(Theta0(eta), eta, eta))*Theta1(eta)*Rd*(3*(Tp-1)+6*(Tp-1)^2*Theta0(eta)+3*(Tp-1)^3*Theta0(eta)^2)/(3*a3*Pr)+4*Rd*(diff(Theta0(eta), eta))^2*(6*(Tp-1)^2*Theta1(eta)+6*(Tp-1)^3*Theta0(eta)*Theta1(eta))/(3*a3*Pr)+4*Rd*(diff(Theta1(eta), eta))*(diff(Theta0(eta), eta))*(6*(Tp-1)+6*(Tp-1)^3*Theta0(eta)^2+12*(Tp-1)^2*Theta0(eta))/(3*a3*Pr)]:
 > (LB,UB):= (0,1): BCs:= [      U0(0) = 0, U1(0) = 0, N0(0) = 0, N1(0) = 0, Theta0(0) = 0, Theta1(0) = 0, U0(1) = 0, U1(1) = 0, N0(1) = 0, N1(1) = 0, Theta0(1) = 1, Theta1(1) = 0 ]:
 >
 > Params:= Record(        M1=  1.2, Rd=0.8,la=0.8,n1=1.2,Q=0.2,Pj=0.001,Ra=0.8,Ec=1,    Pr= 21,   delta= 0.2,    t1= (1/4)*Pi, lambda0=2,lambda1=3,   Qc= 0.1,    l1= 1,K=0.4,H=3 ,deltat=0.05  ):
 > NBVs:= [      a1**D(U0)(0) = `C*__f` , # Skin friction coefficient  (a4+(4*Rd*(1/3))*(1+(Tp-1)*(Theta0(0)+0.1e-2*exp(l1*t1)*Theta1(0)))^3)*((D(Theta0))(0)+0.1e-2*exp(l1*t1)*(D(Theta1))(0)) = `Nu*`    # Nusselt number      ]: Nu:= `Nu*`: Cf:= `C*__f`:
 > Solve:= module() local    nbvs_rhs:= rhs~(:-NBVs), #just the names    Sol, #numeric dsolve BVP solution of any 'output' form    ModuleApply:= subs(       _Sys= {:-ODEs[], :-BCs[], :-NBVs[]},       proc({           M1::realcons:=  Params:-M1,          Pr::realcons:= Params:-Pr,          Rd::realcons:= Params:-Rd,          la::realcons:= Params:-la,          Tp::realcons:= Params:-Tp,          n1::realcons:= Params:-n1,          Q::realcons:= Params:-Q,          Pj::realcons:= Params:-Pj,          Ra::realcons:= Params:-Ra,          Ec::realcons:= Params:-Ec,          t1::realcons:=  Params:-t1,          delta::realcons:= Params:-delta,          lambda0::realcons:= Params:-lambda0,          lambda1::realcons:= Params:-lambda1,          Qc::realcons:= Params:-Qc,          K::realcons:= Params:-K,          l1::realcons:= Params:-l1,          H::realcons:= Params:-H       })          Sol:= dsolve(_Sys, _rest, numeric);          AccumData(Sol, {_options});          Sol       end proc    ),    AccumData:= proc(       Sol::{Matrix, procedure, list({name, function}= procedure)},       params::set(name= realcons)    )    local n, nbvs;       if Sol::Matrix then          nbvs:= seq(n = Sol[2,1][1,Pos(n)], n= nbvs_rhs)       else          nbvs:= (nbvs_rhs =~ eval(nbvs_rhs, Sol(:-LB)))[]       fi;       SavedData[params]:= Record[packed](params[], nbvs)    end proc,    ModuleLoad:= eval(Init); export    SavedData, #table of Records    Pos, #Matrix column indices of nbvs    Init:= proc()       Pos:= proc(n::name) option remember; local p; member(n, Sol[1,1], 'p'); p end proc;       SavedData:= table();       return    end proc ;    ModuleLoad() end module:
 >
 >
 > #procedure that generates 3-D plots (dropped-shadow contour + surface) of an expression ParamPlot3d:= proc(    Z::{procedure, `module`}, #procedure that extracts z-value from Solve's dsolve solution    X::name= range(realcons), #x-axis-parameter range    Y::name= range(realcons), #y-axis-parameter range    FP::list(name= realcons), #fixed values of other parameters    {       #fraction of empty space above and below plot (larger "below"       #value improves view of dropped-shadow contourplot):       zmargin::[realcons,realcons]:= [.05,0.15],       eta::realcons:= :-LB, #independent variable value       dsolveopts::list({name, name= anything}):= [],       contouropts::list({name, name= anything}):= [],       surfaceopts::list({name, name= anything}):=[]        } ) local    LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),    Zremember:= proc(x,y)    option remember; #Used because 'grid' should be the same for both plots.       Z(          Solve(             LX= x, LY= y, FP[],             #Default dsolve options can be changed by setting 'dsolveopts':             'abserr'= 0.5e-7, 'interpolant'= false, 'output'= Array([eta]),               dsolveopts[]          )       )    end proc,    plotspec:= (Zremember, RX, RY),    C:= plots:-contourplot(       plotspec,       #These default plot options can be changed by setting 'contouropts':       'grid'= [25,25], 'contours'= 5, 'filled',       'coloring'= ['yellow', 'orange'], 'color'= 'green',       contouropts[]    ),    P:= plot3d(       plotspec,       #These default plot options can be changed by setting 'surfaceopts':       'grid'= [25,25], 'style'= 'surfacecontour', 'contours'= 6,       surfaceopts[]    ),    U, L #z-axis endpoints after margin adjustment ;    #Stretch z-axis to include margins:    (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(       zmargin[],       (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range    );    plots:-display(       [          plots:-spacecurve(             {                [[lhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),L]], #yz backwall                [[rhs(RX),rhs(RY),U],[rhs(RX),lhs(RY),U],[rhs(RX),lhs(RY),L]]  #xz backwall             },             'color'= 'grey', 'thickness'= 0          ),          plottools:-transform((x,y)-> [x,y,L])(C), #dropped-shadow contours          P       ],       #These default plot options can be changed simply by putting the option in the       #ParamPlot3d call:       'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 'axes'= 'frame',       'labels'= [lhs(X), lhs(Y), Z], 'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],       'caption'= nprintf(cat("%a = %4.2f, "\$nops(FP)-1, "%a = %4.2f"), (lhs,rhs)~(FP)[]),       'captionfont'= ['TIMES', 14],       'projection'= 2/3,          _rest    ) end proc:
 >
 >
 >
 > ParamPlot3d(    GetNu,Q= 0..5, Rd= 0..5, [        Pr= 21   ],    labels= [Q, gamma, Nu] );