Hi,

I get the error in the following code

restart:

gama1:=0.01:

zet:=0;

#phi0:=0.00789:

Phiavg:=0.02;

lambda:=0.01;

Ha:=1;

0

0.02

0.01

1

rhocu:=2/(1-zet^2)*int((1-eta)*rho(eta)*c(eta)*u(eta),eta=0..1-zet):

eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(1-Ha^2*u(eta))+((1/(eta)+1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));

eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*(-2/(1-zet^2)*rho(eta)*c(eta)*u(eta)/(p2*10000)+( (a[k1]+2*b[k1]*phi(eta))/(1+a[k1]*phi1[w]+b[k1]*phi1[w]^2)*diff(phi(eta),eta)+k(eta)/k1[w]/(eta)*diff(T(eta),eta) ));

eq3:=diff(phi(eta),eta)+phi(eta)/(N[bt]*(1+gama1*T(eta))^2)*diff(T(eta),eta);

/ d / d \\ mu1[w] (1 - u(eta))

|----- |----- u(eta)|| + -------------------

\ deta \ deta // mu(eta)

/ / d \\

| mu_phi |----- phi(eta)||

| 1 \ deta /| / d \

+ |--- + -----------------------| |----- u(eta)|

\eta mu(eta) / \ deta /

/ /

| |

/ d / d \\ 1 | | rho(eta) c(eta) u(eta)

|----- |----- T(eta)|| + ------ |k1[w] |- ----------------------

\ deta \ deta // k(eta) | | 5000 p2

\ \

/ d \

(a[k1] + 2 b[k1] phi(eta)) |----- phi(eta)|

\ deta /

+ -------------------------------------------

2

1 + a[k1] phi1[w] + b[k1] phi1[w]

/ d \\\

k(eta) |----- T(eta)|||

\ deta /||

+ ---------------------||

k1[w] eta ||

//

/ d \

phi(eta) |----- T(eta)|

/ d \ \ deta /

|----- phi(eta)| + ------------------------

\ deta / 2

N[bt] (1 + 0.01 T(eta))

mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):

k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):

rhop:=3880:

rhobf:=998.2:

cp:=773:

cbf:=4182:

rho:=unapply( phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):

c:=unapply( (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):

mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):

a[mu1]:=39.11:

b[mu1]:=533.9:

mu1[bf]:=9.93/10000:

a[k1]:=7.47:

b[k1]:=0:

k1[bf]:=0.597:

zet:=0.5:

#phi(0):=1:

#u(0):=0:

phi1[w]:=phi0:

N[bt]:=0.2:

mu1[w]:=mu(0):

k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1):

eq2:=subs(phi(0)=phi0,eq2):

eq3:=subs(phi(0)=phi0,eq3):

#A somewhat speedier version uses the fact that you really need only compute 2 integrals not 3, since one of the integrals can be written as a linear combination of the other 2:

Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10,B;

global Q1,Q2;

print(pp2,fi0);

if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:

res := dsolve(subs(p2=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=lambda/(phi(1)*rhop/rhobf+(1-phi(1)))*D(u)(1),D(u)(0)=0,phi(1)=phi0,T(1)=0,D(T)(1)=1}), numeric,output=listprocedure):

F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):

INT0:=evalf(Int((1-eta)*F0(eta),eta=0..1-zet));

INT10:=evalf(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet));

B:=(-cbf*rhobf+cp*rhop)*INT10+ rhobf*cbf*INT0;

a[1]:=2/(1-zet^2)*B-10000*pp2;

a[2]:=INT10/INT0-Phiavg;

Q1(_passed):=a[1];

Q2(_passed):=a[2];

if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if

end proc;

#The result agrees very well with the fsolve result.

#Now I did use a better initial point. But if I start with the same as in fsolve I get the same result in just about 2 minutes, i.e. more than 20 times as fast as fsolve:

Q1:=proc(pp2,fi0) Q[1](_passed) end proc;

Q2:=proc(pp2,fi0) Q[2](_passed) end proc;

Optimization:-LSSolve([Q1,Q2],initialpoint=[6.5,exp(-1/N[bt])]);

proc(pp2, fi0) ... end;

proc(pp2, fi0) ... end;

proc(pp2, fi0) ... end;

HFloat(6.5), HFloat(0.006737946999)

the error is :

Error, (in Optimization:-LSSolve) system is singular at left endpoint, use midpoint method instead

how can I fix it.

Thanks

**Amir**