Question: Improve Algorithm Dsolve

Hi,

I have the following code which solves a system od ode (3 ode) with constraint parameter. actually, this algorithm is time consuming and is not effective as I see that some people solve this odes with matlab in an efficient (qucik) way. So, I ask you if it is possible help me to improve the efficiency of the following code

thanks

restart:
eq1 := diff(u(eta), eta, eta)+diff(u(eta), eta)/ (eta+zet)+Gr*Pr*T(eta)-Nr*phi(eta)+sigma1;
eq2 := diff(T(eta), eta, eta)+diff(T(eta), eta)/(eta+zet)+Nb*diff(T(eta), eta)*diff(phi(eta), eta)+Nt*diff(T(eta), eta)^2-u(eta);
eq3 := diff(phi(eta), eta,eta)+diff(phi(eta), eta,eta)/(eta+zet)+ Nt/Nb*(diff(T(eta), eta, eta)+diff(T(eta), eta)/(eta+zet))-Le*u(eta);
zet:=0.1:
Nt:=0.05;
Nb:=0.05;
Le:=10;
Pr:=1;
Gr:=5;
Nr:=0.1;
Ree:=5;
#sigma:=2;



p:=proc(pp2) option remember; global res,F0,F1,F2;
if not type([pp2],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({subs(sigma1=pp2,eq1)=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0,phi(1-zet)=0,T(0)=0,T(1-zet)=0}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
evalf((2/(1-zet^2))*Int(F0(eta)*(eta+zet),eta=0..1-zet))-Ree*Pr;
end proc;

fsolve(p(pp2)=0,pp2=(-80)..(80));   # I dont know what exactly this parameter has
#res := dsolve({eq1=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0,phi(1-zet)=0,T(0)=0,T(1-zet)=0}, numeric);
 

sigma1:=%;
                             
res := dsolve({eq1=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0,phi(1-zet)=0,T(0)=0,T(1-zet)=0}, numeric);
with(plots):
odeplot(res,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1-zet);

Amir