@acer 

Please ignore my last message. I spoke it over with my tutor and it isn't really applicable to the system. The problems that I seemed to be having were tied to some of the physical values.

Even if I was never able to make both Po2m and Pco2m work at the same time in the simulation I was able to use them to make a corrleation to the parameters which I later used for the simulation.

Once again thank you so much for all your help!

@acer 

Please ignore my last message. I spoke it over with my tutor and it isn't really applicable to the system. The problems that I seemed to be having were tied to some of the physical values.

Even if I was never able to make both Po2m and Pco2m work at the same time in the simulation I was able to use them to make a corrleation to the parameters which I later used for the simulation.

Once again thank you so much for all your help!

@acer 

Once again I want to thank you for taking your time on the problem... you've been immensely helpful. However I've made a huge mistake on the function I want out of the first differential equation (incorrect physical interpretation) and I think now it should be able to work. 

However, I'm getting trouble with the "allowing the parameters to be function"  part. 

In essence what I would want to have now is:

For O2:

DPo2m=(Fo2s(z)*P(z)/(Ho2*(Fo2s(z)+Fco2s(z)+Fn2))-aValue)

so I require the logarithmic mean value of this function : (do2(0)-do2(zf))/ln(do2(0)/abs(do2(zf))           note: the abs is not normally present but I need it because there tends to be a slight invertion of the function which gives negative values.

I'm trying to do this calculation with the same procedure that I had for evaulating the integral but it doesnt seem to work.

Here is the latest worksheet: tried to put as much enfasis where I'm having the problem

Final_Mod.mw (note that the range of 48 is rather arbitrary, could be more or less)

Again and again, I cannot thank you enough for your help

@acer 

Once again I want to thank you for taking your time on the problem... you've been immensely helpful. However I've made a huge mistake on the function I want out of the first differential equation (incorrect physical interpretation) and I think now it should be able to work. 

However, I'm getting trouble with the "allowing the parameters to be function"  part. 

In essence what I would want to have now is:

For O2:

DPo2m=(Fo2s(z)*P(z)/(Ho2*(Fo2s(z)+Fco2s(z)+Fn2))-aValue)

so I require the logarithmic mean value of this function : (do2(0)-do2(zf))/ln(do2(0)/abs(do2(zf))           note: the abs is not normally present but I need it because there tends to be a slight invertion of the function which gives negative values.

I'm trying to do this calculation with the same procedure that I had for evaulating the integral but it doesnt seem to work.

Here is the latest worksheet: tried to put as much enfasis where I'm having the problem

Final_Mod.mw (note that the range of 48 is rather arbitrary, could be more or less)

Again and again, I cannot thank you enough for your help

First I want to thank you for your time and the proposed modifications. I would be lying if I said I understood everything but I do get most of it. 

Thankfuly the system at least it converges in a reasonable time now (again thank you very much). I'm a bit worried that I'm getting negative values of O2, is there anyway of setting up a boundry? I've tried another "physical" method that seems to be working but i just wanted to check.

Also I forgot to put in Pco2m in my system when I posted the question equation. I was trying to "edge up" to the full problem. DCo2 should read:

DCo2:=diff(Co2(t),t)=rco2-kla*(Co2(t)-Pco2m(O2(t),Co2(t),ro2)/Hco2);

I've tried using the same procedure por Pco2m as you did for Po2m but it won't converge. 

with regards to some of your comments. I was using in the integral method _d01ajc because it seemed to give me the fastest response time for the integral (which I thought was my problem initially) so maybe there is no need to work around the method, just simply omit it.

with regards to the results the graphs for O2 are a bit werid but it is also heavy on initial hypothesis in the deduction of the equations. So it could just be the aproximation.

Here is the code I'm using (you can ignore the first part since its only how to get Af and Gf)

Final_Modelisation_4.mw

Again I cannot stress this enough, thank you very much for your VERY quick and usefull answer. Even if I don't resolve the Pco2m issue, you've been of great help

First I want to thank you for your time and the proposed modifications. I would be lying if I said I understood everything but I do get most of it. 

Thankfuly the system at least it converges in a reasonable time now (again thank you very much). I'm a bit worried that I'm getting negative values of O2, is there anyway of setting up a boundry? I've tried another "physical" method that seems to be working but i just wanted to check.

Also I forgot to put in Pco2m in my system when I posted the question equation. I was trying to "edge up" to the full problem. DCo2 should read:

DCo2:=diff(Co2(t),t)=rco2-kla*(Co2(t)-Pco2m(O2(t),Co2(t),ro2)/Hco2);

I've tried using the same procedure por Pco2m as you did for Po2m but it won't converge. 

with regards to some of your comments. I was using in the integral method _d01ajc because it seemed to give me the fastest response time for the integral (which I thought was my problem initially) so maybe there is no need to work around the method, just simply omit it.

with regards to the results the graphs for O2 are a bit werid but it is also heavy on initial hypothesis in the deduction of the equations. So it could just be the aproximation.

Here is the code I'm using (you can ignore the first part since its only how to get Af and Gf)

Final_Modelisation_4.mw

Again I cannot stress this enough, thank you very much for your VERY quick and usefull answer. Even if I don't resolve the Pco2m issue, you've been of great help

@acer 

the first sytem a parametrization of an elipse x=a*cos(c*z), y=b*sin(c*z)

diff(x,z)=a*c*-sin(ct)=-a*c(y/b)

diff(y,z)=b*c (x/a)

the medium values I am not having a problem with at all. 

xm=a*sin(c*deltaz)/c*deltaz, which i decided to call α*a since sin(c*deltaz)/c*deltaz is constant

     so xm=α*a 

ym=b*(1-cos(c*deltaz))/c*deltaz

    similarly ym=β*b 

However when y try to use these solutions (Xm and Ym) into yet again a "known answer" problem this time the answer doesn't seem very precise at all

So this is my test:

I'll call

a=β*cos(t)----> da/dt=-β*sin(t)--->dx/dt=-(β/α)*b=Xm/α

b=α*sin(t)---> db/dt=α*cos(t)---> dx/dt=(α/β)*a=Ym/β

To make things very easy I say α=β  therefore I use (arbitrarily) the values c=1/2 and deltaz=Pi

thus α=β=2/Pi

Therefore the system is

da/dt=2/Pi*Xm(a(t),b(t),1/2.,Pi)  with a(0)=.6366198 or a(0.01)=.6365879
db/dt=2/Pi*Ym(a(t),b(t),1/2,Pi)   with b(0)=0  or b(0.01)=.--63661

So my code is as follows (the expresions for Xm and Ym are left unchanged since they give very accurate results)

>nsol2:= dsolve({a(0.01)=.6365879, b(0.01)=.0063661, diff(a(t),t)=-Pi/2*Ym(a(t),b(t),1/2.,Pi), diff(b(t),t)=Pi/2Xm(a(t),b(t),1/2.,Pi)},numeric, relerr=1e-9, 'output'=listprocedure, known=[Xm,Ym]):
> as:=eval(a(t),nsol2):
> bs:=eval(b(t),nsol2):

I used the 0.01 since I was still getting problems at a(0) or b(0)

I tested a couple of points afterwards and I'm getting good results now (I made a post with wrong answeres, double checked the equations and found a couple of errors in my code)

I'm extremely greatful for all the help you guys have provided. 

@acer 

the first sytem a parametrization of an elipse x=a*cos(c*z), y=b*sin(c*z)

diff(x,z)=a*c*-sin(ct)=-a*c(y/b)

diff(y,z)=b*c (x/a)

the medium values I am not having a problem with at all. 

xm=a*sin(c*deltaz)/c*deltaz, which i decided to call α*a since sin(c*deltaz)/c*deltaz is constant

     so xm=α*a 

ym=b*(1-cos(c*deltaz))/c*deltaz

    similarly ym=β*b 

However when y try to use these solutions (Xm and Ym) into yet again a "known answer" problem this time the answer doesn't seem very precise at all

So this is my test:

I'll call

a=β*cos(t)----> da/dt=-β*sin(t)--->dx/dt=-(β/α)*b=Xm/α

b=α*sin(t)---> db/dt=α*cos(t)---> dx/dt=(α/β)*a=Ym/β

To make things very easy I say α=β  therefore I use (arbitrarily) the values c=1/2 and deltaz=Pi

thus α=β=2/Pi

Therefore the system is

da/dt=2/Pi*Xm(a(t),b(t),1/2.,Pi)  with a(0)=.6366198 or a(0.01)=.6365879
db/dt=2/Pi*Ym(a(t),b(t),1/2,Pi)   with b(0)=0  or b(0.01)=.--63661

So my code is as follows (the expresions for Xm and Ym are left unchanged since they give very accurate results)

>nsol2:= dsolve({a(0.01)=.6365879, b(0.01)=.0063661, diff(a(t),t)=-Pi/2*Ym(a(t),b(t),1/2.,Pi), diff(b(t),t)=Pi/2Xm(a(t),b(t),1/2.,Pi)},numeric, relerr=1e-9, 'output'=listprocedure, known=[Xm,Ym]):
> as:=eval(a(t),nsol2):
> bs:=eval(b(t),nsol2):

I used the 0.01 since I was still getting problems at a(0) or b(0)

I tested a couple of points afterwards and I'm getting good results now (I made a post with wrong answeres, double checked the equations and found a couple of errors in my code)

I'm extremely greatful for all the help you guys have provided. 

@pagan4515

Thank you very much for the solution. It's exactly what I needed!

I'm creeping iup the complexity up to the simulation I have to do, but when I try to apply the same concept but to a system of differential equations. During the "second system"  I get the following error message 

Error, (in f) unable to store '.33333333333333*Int(xs,0. .. 3.)' when datatype=float[8]

This is the code I'm using:

restart;
> sys:= [x(0)=a,y(0)=0, diff(x(t),t)=-a*c/(b)*y(t), diff(y(t),t)=b*c/(a)*x(t)];
> nsol:= dsolve(sys,numeric, parameters=[a,b,c], 'output'=listprocedure):
> xs:=eval(x(t),nsol):
> ys:=eval(y(t),nsol):

Xm:=proc(aValue,bValue,cValue,zValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)and type(zValue,numeric) ) then
> return 'procname'(args);
> end if;
> xs('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[14](1/zValue*Int(xs,0..zValue));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:
>
> Ym:=proc(aValue,bValue,cValue,zValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)and type(zValue,numeric) ) then
> return 'procname'(args);
> end if;
> ys('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[14](1/zValue*Int(ys,0..zValue));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:

## Up to this point everying seems to be working just fine (tested for some random values of a, b, and c and got back the right values

nsol2:= dsolve({a(0)=0, diff(a(t),t)=Xm(a(t),1,2.,3.)},numeric, relerr=1e-7, 'output'=listprocedure, known=[Xm]):
> as:=eval(a(t),

nsol2):Error, (in f) unable to store '.33333333333333*Int(xs,0. .. 3.)' when datatype=float[8]

Any ideas

(sorry as I said previously I have practically 0 training in programming)

Thank you very much in advance

@pagan4515

Thank you very much for the solution. It's exactly what I needed!

I'm creeping iup the complexity up to the simulation I have to do, but when I try to apply the same concept but to a system of differential equations. During the "second system"  I get the following error message 

Error, (in f) unable to store '.33333333333333*Int(xs,0. .. 3.)' when datatype=float[8]

This is the code I'm using:

restart;
> sys:= [x(0)=a,y(0)=0, diff(x(t),t)=-a*c/(b)*y(t), diff(y(t),t)=b*c/(a)*x(t)];
> nsol:= dsolve(sys,numeric, parameters=[a,b,c], 'output'=listprocedure):
> xs:=eval(x(t),nsol):
> ys:=eval(y(t),nsol):

Xm:=proc(aValue,bValue,cValue,zValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)and type(zValue,numeric) ) then
> return 'procname'(args);
> end if;
> xs('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[14](1/zValue*Int(xs,0..zValue));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:
>
> Ym:=proc(aValue,bValue,cValue,zValue)
> option remember;
> if not (type(aValue,numeric) and type(bValue,numeric) and type(cValue,numeric)and type(zValue,numeric) ) then
> return 'procname'(args);
> end if;
> ys('parameters'= [a = aValue, b=bValue,c=cValue]);
> try
> forget('evalf/int'); forget(evalf);
> evalf[14](1/zValue*Int(ys,0..zValue));
> catch "unable to store":
> catch "not possible":
> end try;
> end proc:

## Up to this point everying seems to be working just fine (tested for some random values of a, b, and c and got back the right values

nsol2:= dsolve({a(0)=0, diff(a(t),t)=Xm(a(t),1,2.,3.)},numeric, relerr=1e-7, 'output'=listprocedure, known=[Xm]):
> as:=eval(a(t),

nsol2):Error, (in f) unable to store '.33333333333333*Int(xs,0. .. 3.)' when datatype=float[8]

Any ideas

(sorry as I said previously I have practically 0 training in programming)

Thank you very much in advance

@Markiyan Hirnyk 

So I came up with a simpler system than the one I have but should emulate it well enough. On the plus side it has an easy analytical solution to double check our procedure

The system is as follows:

First we solve:

dx/dz=a      ----> x=a*z+b 
For my current problem I need the average value of the function on the interval (0,Zo) so we end up with

xm= aZo/2+b

The second differential equation is:

da/dt=xm  which is solved very easily by hand

say we use as initial condition a(0)=0

My first attempt was as follows:

 > eq1:=diff(x(z),z)=a;
 > sol1:=dsolve({eq1,x(0)=2},numeric,parameters=[a],'output'=listprocedure):
 > xs:=rhs(sol1[2]):   #this last function lets me "obtain" the x function out of the procedure output

Next a procedure to generate the average value of the function:

>Xm:=proc(aValue,zValue) global xs, a, x, z; xs('parameters'= [a = aValue]); 

try return eval(1/zValue*int(xs(z),z=0..zValue)); catch: return "not possible"; 

end try; end proc:

Up to this point I've checked every step with real values and have gotten the right aswer,

so the "functions" (I'm not accustomed to programming language so don't know if it's functions

or procedures) work properly

 Finally I try inserting it into my last differential equation:

     >eq2:=diff(a(t),t)=Xm(a(t),2);
     > sol2:=dsolve({eq2,a(0)=0},numeric,'output'=listprocedure):
     > as:=rhs(sol2[2]);

and i get the following responses:

     Error, (in dsolve/numeric/process_parameters) parameter values must evaluate to 

numeric, got a = a(t)

     Error, (in dsolve/numeric/process_input) system must be entered as a set/list 

of expressions/equations

     Error, invalid input: rhs received sol2[2], which is not valid for its 1st argument, 

expr

So up to this point this was my original problem that I posted. In retrospect I should have 

come up with it before

Now I try the "second thought" you suggested

 Again I try:

> sol3:=b->dsolve({diff(x(z),z)=b,x(0)=2},numeric,'output'=listprocedure):
> xs3:=rhs(sol3[2](b))

However I get the following response:

Warning, The use of global variables in numerical ODE problems is deprecated, and will be 

removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

Error, invalid input: rhs received [z = proc (z) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i;

 option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error 

"invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15;

 if _EnvInFsolve = true then _xout := evalf[_EnvDSNumer..., which is not valid for its 1st argument, expr

I've tried different combinations to "recuperate" my x solution and I get similar answers 

(again sorry if the language isn't precise enough) 

Hope my problem has become clearer and thank you for the follow up

@Markiyan Hirnyk 

So I came up with a simpler system than the one I have but should emulate it well enough. On the plus side it has an easy analytical solution to double check our procedure

The system is as follows:

First we solve:

dx/dz=a      ----> x=a*z+b 
For my current problem I need the average value of the function on the interval (0,Zo) so we end up with

xm= aZo/2+b

The second differential equation is:

da/dt=xm  which is solved very easily by hand

say we use as initial condition a(0)=0

My first attempt was as follows:

 > eq1:=diff(x(z),z)=a;
 > sol1:=dsolve({eq1,x(0)=2},numeric,parameters=[a],'output'=listprocedure):
 > xs:=rhs(sol1[2]):   #this last function lets me "obtain" the x function out of the procedure output

Next a procedure to generate the average value of the function:

>Xm:=proc(aValue,zValue) global xs, a, x, z; xs('parameters'= [a = aValue]); 

try return eval(1/zValue*int(xs(z),z=0..zValue)); catch: return "not possible"; 

end try; end proc:

Up to this point I've checked every step with real values and have gotten the right aswer,

so the "functions" (I'm not accustomed to programming language so don't know if it's functions

or procedures) work properly

 Finally I try inserting it into my last differential equation:

     >eq2:=diff(a(t),t)=Xm(a(t),2);
     > sol2:=dsolve({eq2,a(0)=0},numeric,'output'=listprocedure):
     > as:=rhs(sol2[2]);

and i get the following responses:

     Error, (in dsolve/numeric/process_parameters) parameter values must evaluate to 

numeric, got a = a(t)

     Error, (in dsolve/numeric/process_input) system must be entered as a set/list 

of expressions/equations

     Error, invalid input: rhs received sol2[2], which is not valid for its 1st argument, 

expr

So up to this point this was my original problem that I posted. In retrospect I should have 

come up with it before

Now I try the "second thought" you suggested

 Again I try:

> sol3:=b->dsolve({diff(x(z),z)=b,x(0)=2},numeric,'output'=listprocedure):
> xs3:=rhs(sol3[2](b))

However I get the following response:

Warning, The use of global variables in numerical ODE problems is deprecated, and will be 

removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

Error, invalid input: rhs received [z = proc (z) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i;

 option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error 

"invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15;

 if _EnvInFsolve = true then _xout := evalf[_EnvDSNumer..., which is not valid for its 1st argument, expr

I've tried different combinations to "recuperate" my x solution and I get similar answers 

(again sorry if the language isn't precise enough) 

Hope my problem has become clearer and thank you for the follow up

Thank you for your answer Markiyan. However  I think i've tried it and more recent versions of Maple "dislike"this approach and reccomend using the parameters option.

Just to make sure I'll try it out again

Thank you for your answer Markiyan. However  I think i've tried it and more recent versions of Maple "dislike"this approach and reccomend using the parameters option.

Just to make sure I'll try it out again