Question: How to plot this equation

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})