It may then be more suitable to use a piecewise approximation to the solution - this allows it to be used more flexbly later on. So for your problem, generate the piecewise solution for a given set of parameters, and z=0..2. As a simple example: > sol := op([2,2],dsolve({diff(y(x),x)=-y(x),y(0)=1},numeric, > range=0..2,output=piecewise)); sol := { { undefined , x < 0. 0.999621690653975 - 0.972619698891573 x 2 + 0.486310338843125 (x - 0.0277621231571284) 3 - 0.162115462190059 (x - 0.0277621231571284) 4 + 0.0402113832410770 (x - 0.0277621231571284) , ... Is this suitable?

It may then be more suitable to use a piecewise approximation to the solution - this allows it to be used more flexbly later on. So for your problem, generate the piecewise solution for a given set of parameters, and z=0..2. As a simple example: > sol := op([2,2],dsolve({diff(y(x),x)=-y(x),y(0)=1},numeric, > range=0..2,output=piecewise)); sol := { { undefined , x < 0. 0.999621690653975 - 0.972619698891573 x 2 + 0.486310338843125 (x - 0.0277621231571284) 3 - 0.162115462190059 (x - 0.0277621231571284) 4 + 0.0402113832410770 (x - 0.0277621231571284) , ... Is this suitable?

To answer each of the questions in turn: (1) why cannot I evaluate the value of f(t) at t=5 by using f(5)? You just need to do this a little differently. Consider: > A := 100; v := 50; rho := 0.996; r := 0.9; theta := 0.993; A := 100 v := 50 rho := 0.996 r := 0.9 theta := 0.993 > x :=5000; c := 200; m := 100000; l := 5000; x := 5000 c := 200 m := 100000 l := 5000 > solution := dsolve([ > diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, > diff(f(t),t)=((c - s(t)) / l) * m + x, > diff(h(t),t) = x, > f(0) = 0, s(0) = 0, h(0) = 0], > numeric, output=operator); bytes used=8307108, alloc=917336, time=0.41 solution := [t = (proc(t) ... end proc), f = (proc(t) ... end proc), h = (proc(t) ... end proc), s = (proc(t) ... end proc)] > f := eval(f,solution); f := proc(t) ... end proc So 'f' is now a procedure that can be used to evaluate the 'f' values at different points: > f(100); 7 -0.399982377673172718 10 > f(1000); 9 -0.489993379968400420 10 Please be warned though that it is not such a good idea to assign 'f' then use it again. Generally it is better to use something else like 'fsol' for 'f'. Note that if you try to apply dsolve again your original system of DEs is no longer a system of DEs, but is a system with a 'f' function, and attempting to specify the initial condition for 'f' gives an identity: > f(0) = 0; 0. = 0 (2) For each fixed x, there are curves s(t), f(t), and h(t). Given s(tau) = c, I want to find f(tau)=?. How can I do that? Do I need to find tau first, then find f(tau)? How? Yes, continuing on from prior, we have: > s := eval(s,solution); s := proc(t) ... end proc > fsolve(s(tau)=100); -0.6237217782 > val := fsolve(s(tau)=100,tau=0..100); val := 3.000783529 > f(val); 24992.3520950843448 So use 'fsolve' to compute the tau value at which s(tau)=100, then evaluate with 'f'. (3) My ultimate goal is for x in a range, say x=500..5000, find the minimum of f(tau) (since for each x there is a f(tau)). How can I do that? The easiest way to do this is to delay the evaluation of the DE solution, leaving 'x' a variable, then create an interface procedure to evaluate it at a specific time value, or as follows: > A := 100; v := 50; rho := 0.996; r := 0.9; theta := 0.993; A := 100 v := 50 rho := 0.996 r := 0.9 theta := 0.993 > c := 200; m := 100000; l := 5000; c := 200 m := 100000 l := 5000 > solution := dsolve([ > diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, > diff(f(t),t)=((c - s(t)) / l) * m + x, > diff(h(t),t) = x, > f(0) = 0, s(0) = 0, h(0) = 0], > numeric, output=operator); solution := [t = (proc(t) ... end proc), f = (proc(t) ... end proc), h = (proc(t) ... end proc), s = (proc(t) ... end proc)] > f := eval(f,solution); f := proc(t) ... end proc Now the 'f' solution cannot be evaluated until 'x' is assigned: > f(1); Error, (in f) global 'x' must be assigned to a numeric value before obtaining a solution Now the interface procedure, say you want values at t=100, then do: > tval := 100: > getval := proc(xval) global x; > if type(xval,'numeric') then > x := xval; f(tval); > else > 'getval'(xval); > end if; > end proc: Now 'getval' is a function that computes f(t=100;x), and it can be plotted: > plot(getval,500..5000); * -2.4e+06 *A -2.6e+06 +A -2.8e+06 + A + A -3e+06 + A -3.2e+06 + AA -3.4e+06 + AA -3.6e+06AA + AA -3.8e+06 AA + AA + AAA -4e+06 AAAAAA AAAAAAAAAAAAAAAAAAAAA +-+--+--+--+---+--+-**********************************+---+--+--+--+---+--+ 1000 2000 3000 4000 5000 Or solved: > x := 'x': > Optimization[Minimize](getval(x),x=500..5000); 7 [-0.413681014780686982 10 , [2480.90973324573588 = 2480.90973324573588]] - Allan