Question: System taking too long to solve

I have a rather complicated procedure to solve. All the individual parts of the program are working fine, but when I try to put them together the simulation takes too long to solve. To my knowledge what I'm trying to do is not complicated in it of itself; each system is rather simple. But the whole does not seem to be working or is simply taking way longer than expected or is impractical for the task.

I'm supposed to be doing an integration over a range of 48 (48 hours if units are of interest) and so far the calculation seems to have taken 50 minutes with no answer. Just doing it over 1 hour range takes about 15 minutes. I've tried to check if there is a units problem (thus some factor making the calculations bizarre) but so far I've found none.

I would like it to be much quicker if possible since many of the "constants" involved must be tested to optimize.

Since the procedure is rather lengthy I'll try to go along and descibe what is going on before presenting the actual code.

First there is a system of differential equations which should give me back Fo2(z), Fco2(z). However the value to be used on the next step is the average value of a function of these functions Po2(z)=Fo2(z)/(Fo2(z)+Fco2(z)+Fn2)*P(z)-----> we obtain the average value of the function Po2m and Pco2m in a close range (0,zf)

However, the value of these functions are dependent on the variables of the next set of differential equations. In  this sense when we obtained Po2m we were looking for Po2m(a,b,c) where a b and c are different variables for the next sytem.

Finaly we wish to resolve this last system of a (o2) b (co2) as well as several other variables (the system of equations is rather large)

As mentioned previously each different step works just fine. I have no problem getting Po2m(a,b,c) and it gives me reasonable values (I'ved tried to reduce as much as possible the integration time)

and the last system, when we omit the Po2 and leave in a constant works just fine and gives me plausible values as well

The real system is as follows (there are more parts to it but are irrelevant to the problem, however I might have defined a variable there which is not present in the rest, if you need it please let me know)

Before presenting the code I would like to thank anyone who is willing to take it a look.

# Constants related to the equations, 
> vmax:=.0108: Km:=10:
> teo:=0.0475: Kte:=0.063:
> Ki:=25: Tc:=25: E:=200:  zf:=4.0:

> Q:=24: 
> Fn2:=.79*Q/(8.314e-5*(25+273.15)): # Used
> Fo2o:=.21*Q/(8.314e-5*(25+273.15)): # Used as initial condition for the  first diff system following system
> eh:=.1:
> kla:=650:
> Ho2:=769.23: #L*atm/mol
> Hco2:=29.41: # L*atm/mol
> Po:=1: #considering 1 atm at top of the reactor
> Ar:=10/zf:

#Ecuations describing the evolution in the gases
> facc:=1+1/(.01/c^3+1):
>P:=z->Po+0.986(zf-z):
> #Factor 1000 is to chance concentration to mol/l to mol/m3
> nsol3:=dsolve({diff(Fo2(z),z)=-(1-eh)*Ar*kla*facc*1000*(Fo2(z)*P(z)/(Ho2*(Fo2(z)+Fco2(z)+Fn2))-a),diff(Fco2(z),z)=(1-eh)*Ar*kla*1000*(b-Fco2(z)*P(z)/(Hco2*(Fo2(z)+Fco2(z)+Fn2))), Fo2(0)=Fo2o, Fco2(0)=0},numeric,'range'=0..zf,parameters=[a,b,c],'output'=listprocedure):
> Fco2s:=eval(Fco2(z),nsol3):
> Fo2s:=eval(Fo2(z),nsol3):

#Procedures in order to obtain Po2m(a,b,c) and Pco2m(a,b,c)
> Po2m:=proc(aValue,bValue,cValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)) then
> return 'procname'(args);
> end if;
> Fo2s('parameters'= [a = aValue, b=bValue,c=cValue]);
> Fco2s('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[6](1/zf*int(Fo2s(z)*P(z)/(Fo2s(z)+Fco2s(z)+Fn2),z=0.001..zf,'method'=_d01ajc));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:
>
> Pco2m:=proc(aValue,bValue,cValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)) then
> return 'procname'(args);
> end if;
> Fo2s('parameters'= [a = aValue, b=bValue,c=cValue]);
> Fco2s('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[6](1/zf*int(Fco2s(z)*P(z)/(Fo2s(z)+Fco2s(z)+Fn2),z=0.001..zf,'method'=_d01ajc));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:

> #More Constants for the following system
> umax:=.6:
> Kgx:=1.5:
> tu:=.0355:
> Ktx:=.0876:
> Kgm:=.75:
> txo:=1.5:
> Ysx:=.55:
> m:=.07:
> Mrx:=26.5: #g/cmol
> Mrg:=30: #g/Cmol

#Ecuaciones
> mu:=umax*G(t)/(Kgx+G(t))*tu*(exp(Ktx*Tc)-txo):
> qg:=mu/Ysx+m*G(t)/(Kgm+G(t)):
> ro2:=(4.2*mu/Mrx-4*qg/Mrg)*X(t)/(-2.): #please check sign, should be in mol/Lh
> rco2:=(qg/Mrg-mu/Mrx)*X(t): #should be in mol/lh
> rA:=vmax*E*A(t)/(Km+A(t))*teo*exp(Kte*Tc)*Ki/(Ki+G(t)):
> facc2:=1+1/(.01/(ro2)^3+1):
> DX:=diff(X(t),t)=mu*X(t):
> DA:=diff(A(t),t)=-rA:
> DG:=diff(G(t),t)=rA-qg*X(t):
> DCo2:=diff(Co2(t),t)=rco2-kla*facc2*(Co2(t)-.001/Hco2):
> Do2:=diff(O2(t),t)=-ro2+kla*facc2*(Po2m(O2(t),Co2(t),ro2)/Ho2-O2(t)):

#Resolution
> nsol4:=dsolve({DX,DA,DG,DCo2,Do2,X(0.00001)=.5,A(0.00001)=Af,G(0.00001)=Gf,Co2(0.00001)=0,O2(0.00001)=0.000273},numeric,maxfun=0,'minstep'=0.001, range=0..48,'output'=listprocedure,known=[Pco2m,Po2m]):
> Ar3:=eval(A(t),nsol4):
> Gr3:=eval(G(t),nsol4):
> Co2r:=eval(Co2(t),nsol4):
> O2r:=eval(O2(t),nsol4):
> Xr:=eval(X(t),nsol4):

I tried to restrain the calculation as much as possible (maxfun had to be equal to 0 since it would stop beforehand if not) by making the step bigger. 

NOTE: in the initial conditions for the final system, you can use Af=105 and Gf=15