In x axis (0..0.5,0.05) and y axis (0..1)

then how to plot this values 

[.486935382154125, .485087274176440, .483255914856304, .481441076124814, .479642533998987, .477860068520125, .476093463645627, .474342507167362, .472606990609939, .470886709216958, .469181461771770]

restart;

PDEtools[declare]((f, g)(x), prime = x);

de1 := diff(f(x), x, x, x, x)-(H*H)*(diff(f(x), x, x))-R*(diff(f(x), x, x))*(diff(f(x), x))+R*(diff(f(x), x, x, x))*f(x);

de2 := diff(g(x), x, x)-(H*H)*g(x)-R*(diff(f(x), x))*g(x)+R*(diff(g(x), x))*f(x); R := 1; H := 1;

dd1 := {de1 = 0, de2 = 0, f(0) = 0, f(1) = 1, g(0) = 1, g(1) = 1, (D(f))(0) = 0, (D(f))(1) = 0};

r1 := dsolve(dd1, numeric, output = Array([0., 0.5e-1, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.00]));

restart;

PDEtools[declare]((f, g)(x), prime = x);

de1 := diff(f(x), x, x, x, x)-(H*H)*(diff(f(x), x, x))-R*(diff(f(x), x, x))*(diff(f(x), x))+R*(diff(f(x), x, x, x))*f(x);

de2 := diff(g(x), x, x)-(H*H)*g(x)-R*(diff(f(x), x))*g(x)+R*(diff(g(x), x))*f(x); R := 5; H := 5;

dd1 := {de1 = 0, de2 = 0, f(0) = 0, f(1) = 1, g(0) = 1, g(1) = 1, (D(f))(0) = 0, (D(f))(1) = 0};

r2 := dsolve(dd1, numeric, output = Array([0., 0.5e-1, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.00]));

odeplot({r1, r2}, [x, f(x)]);
print(`output redirected...`); # input placeholder
                f(x) will now be displayed as f
 

restart; N := 4; de1 := A*(diff(f(eta), eta, eta, eta))+n*(-(diff(f(eta), eta, eta)))^(n-1)*(diff(f(eta), eta, eta, eta))-S*(diff(f(eta), eta))+(2-n)*eta*(diff(f(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(f(eta), eta, eta))/(1+n)-(diff(f(eta), eta))^2-g(eta)*(diff(f(eta), eta, eta))+(M*M)*(diff(f(eta), eta)) = 0, A*(diff(g(eta), eta, eta, eta))+(-(diff(f(eta), eta, eta)))^(n-1)*(diff(g(eta), eta, eta, eta))-(n-1)*(diff(g(eta), eta, eta))*(diff(f(eta), eta, eta, eta))*(-(diff(f(eta), eta, eta)))^(n-2)-S*(diff(g(eta), eta))+(2-n)*eta*(diff(g(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(g(eta), eta, eta))/(1+n)-(diff(g(eta), eta))^2+g(eta)*(diff(g(eta), eta, eta))-(M*M)*(diff(g(eta), eta)) = 0, (1+E*j(eta))*(diff(j(eta), eta, eta))+E*(diff(j(eta), eta))^2+2*Pr*n*f(eta)*g(eta)*(diff(j(eta), eta))/(1+n)-Pr*S*(2-n)*eta*(diff(j(eta), eta))/(1+n)+Pr*(Nb*(diff(j(eta), eta))*(diff(h(eta), eta))+Nt*(diff(j(eta), eta))^2)+Pr*lambda*j(eta) = 0, diff(h(eta), eta, eta)+2*Le*Pr*n*f(eta)*g(eta)*(diff(h(eta), eta))/(1+n)-Le*Pr*S*(2-n)*eta*(diff(h(eta), eta))/(1+n)+Nt*(diff(j(eta), eta, eta))/Nb = 0, f(0) = 0, (D(f))(0) = 1, g(0) = 0, (D(g))(0) = alpha, (D(j))(0) = -b*(1-j(0))/(1+E*j(0)), (D(h))(0) = -d*(1-h(0)), (D(f))(N) = 0, (D(g))(N) = 0, j(N) = 0, h(N) = 0; d1 := subs(alpha = .2, M = .4, A = 1, S = .1, n = .5, Pr = 4, E = 1.5, Nb = .5, Nt = .2, Le = 1, lambda = .2, b = 1.2, d = .5, [de1]); da1 := dsolve(d1, numeric, output = operator, maxmesh = 2048, method = bvp[midrich], abserr = 10); with(plots); restart; N := 4; de2 := A*(diff(f(eta), eta, eta, eta))+n*(-(diff(f(eta), eta, eta)))^(n-1)*(diff(f(eta), eta, eta, eta))-S*(diff(f(eta), eta))+(2-n)*eta*(diff(f(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(f(eta), eta, eta))/(1+n)-(diff(f(eta), eta))^2-g(eta)*(diff(f(eta), eta, eta))+(M*M)*(diff(f(eta), eta)) = 0, A*(diff(g(eta), eta, eta, eta))+(-(diff(f(eta), eta, eta)))^(n-1)*(diff(g(eta), eta, eta, eta))-(n-1)*(diff(g(eta), eta, eta))*(diff(f(eta), eta, eta, eta))*(-(diff(f(eta), eta, eta)))^(n-2)-S*(diff(g(eta), eta))+(2-n)*eta*(diff(g(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(g(eta), eta, eta))/(1+n)-(diff(g(eta), eta))^2+g(eta)*(diff(g(eta), eta, eta))-(M*M)*(diff(g(eta), eta)) = 0, (1+E*j(eta))*(diff(j(eta), eta, eta))+E*(diff(j(eta), eta))^2+2*Pr*n*f(eta)*g(eta)*(diff(j(eta), eta))/(1+n)-Pr*S*(2-n)*eta*(diff(j(eta), eta))/(1+n)+Pr*(Nb*(diff(j(eta), eta))*(diff(h(eta), eta))+Nt*(diff(j(eta), eta))^2)+Pr*lambda*j(eta) = 0, diff(h(eta), eta, eta)+2*Le*Pr*n*f(eta)*g(eta)*(diff(h(eta), eta))/(1+n)-Le*Pr*S*(2-n)*eta*(diff(h(eta), eta))/(1+n)+Nt*(diff(j(eta), eta, eta))/Nb = 0, f(0) = 0, (D(f))(0) = 1, g(0) = 0, (D(g))(0) = alpha, (D(j))(0) = -b*(1-j(0))/(1+E*j(0)), (D(h))(0) = -d*(1-h(0)), (D(f))(N) = 0, (D(g))(N) = 0, j(N) = 0, h(N) = 0; d2 := subs(alpha = .2, M = .4, A = 1, S = .1, n = .5, Pr = 5, E = 1.5, Nb = .5, Nt = .2, Le = 1, lambda = .2, b = 1.2, d = .5, [de2]); da2 := dsolve(d2, numeric, output = operator, maxmesh = 2048, method = bvp[midrich], abserr = 10); with(plots); p1 := odeplot([da1, da2], [eta, diff(f(eta), eta), linestyle = 1, color = "Red", thickness = 2], labels = ["η", "f' g' "], labeldirections = [HORIZONTAL, VERTICAL]); p4 := odeplot(da2, [[eta, f(eta)], [eta, g(eta)], [eta, h(eta)], [eta, j(eta)]], color = [red, green, blue, black]); plots[display]({p1})

Hi, There is a problem in solving ODE using dsolve/numeric code.

Ba.mw

How to get a function value for this HPM code. in the equations p-value 0 to N also some problem.