> 
ODES := diff(f(eta), eta)*diff(f(eta), eta$2)+X*diff(f(eta), eta)+Gr*(theta(eta))+Gc*(phi(eta))M*(S*g(eta)f(eta))2*(S)*diff(f(eta), eta) = 0,
diff(theta(eta), eta$2)2*Pr*(S)*(diff(theta(eta), eta))+Pr*Ec*(diff(f(eta), eta))^2+Pr*Df*(diff(phi(eta), eta$2)) = 0,
diff(phi(eta), eta$2)2*Sc*(S)*(diff(phi(eta), eta))K*phi(eta)+Sr*diff(theta(eta), eta$2)=0, R^2*diff(g(eta), eta$2)+P*diff(f(eta), eta)2*R*P*(S)*diff(g(eta), eta)=0;


(1) 
> 
bcs:= f(0) = 0,
f(1) = 1,
g(0) = 0,
g(1) = 1,
theta(0) = 0,
theta(1) = 1,
phi(0) = 0,
phi(1) = 1


(2) 
> 
params:=[ Sc = 0.22, Pr = 0.71,S=1,K=0.5,Sr=2,Df=2,X=2,
Gr=10,Gc=10,M=0.4,R=6,P=3,Ec=0.22
];


(3) 
> 
phiVals:=[0.01, 0.1, 0.2, 1]:
Mvals:= [3, 5, 7, 9]:
ans:=Matrix( numelems(Mvals)*numelems(phiVals)+1, 5):
ans[1,..]:= < 'M'  'phi'  expr1  expr2 expr3>:
for k from 1 by 1 to 4 do
mv:= Mvals[k]:
for j from 1 by 1 to 4 do
pVal:=phiVals[j]:
sol:=dsolve( eval([ODES, bcs], params), numeric,method=bvp[middefer]);
ans[3*(k1)+j+1,..]:= < mv 
pVal 
R__e^(0.5)*sh= eval( diff(phi(eta), eta), [sol[], params[]])(0) 
R__e^(0.5)*NU= eval( diff(theta(eta), eta),[sol[], params[]])(0) 
R__e^(0.5)*C[f]=eval( diff(f(eta), eta,eta), [sol[], params[]])(0)
>;
od:
od:
ans;


(4) 
> 
for k from 1 by 1 to 3 do
plot( [ seq( [ seq( [ans[j,1], rhs(ans[j,2+k]) ], j=i..10,3 ) ], i=2..4 ) ],
color=[red, green, blue],
labels=[typeset(M), typeset( lhs(ans[2,2+k]) )],
labelfont=[times, bold, 20],
legend=[typeset(phi=0.01),typeset(phi=0.1),typeset(phi=0.2)],
legendstyle=[font=[times, bold, 20]],
title=typeset( ans[1,2+k], " versus ", M, " parameterized by ", phi),
titlefont=[times, bold, 24]
)
od;

> 
sol:=dsolve( eval([ ODES, bcs], params), numeric);

> 
plots:odeplot( sol,
[ [eta, theta(eta)],
[eta, phi(eta)],
[eta, f(eta)]
],
eta=0..10,
color=[red, green, blue]
);

