@acer 

Oh wow! You did exactly what I wanted to do: using a parameter to define the dsolve and then loop through the parameter! Not only the memory is reduced from 550MB to 360MB, more importantly the computation time has been reduced from 100sec to 40sec, which is not what I expected but even better! 

Previously when I used a parameter it will give me an error, that roughly says: using parameter in dsolve/numeric is deprecated, and will be removed in a future release. So apparently it's the compile mode that did the magic. I can't find it under Maple's Help page on dsolve, can you please explain what it does exactly?

Thank you so much this is the best solution!

@acer 

The parameter option can't be used for dsolve/numeric, at least that's my understanding. I'm not sure how to use events, how are they different from piecewises?

The DE system still has 2D Math? Which part? I thought I got rid of the 2D Math after Carl suggested so a few days ago.

@Carl Love 

Yes I'm using Maple 12, so apparently Maple 18 is better for the memory control.

@acer 

Yes the problem occurs while I have suppressed all outputs. Benzene-meta-disulfur-AC-current-loop-iteration.mw

Attached is the file. The huge beast block of input of sys1 is my DE system. sol1 is the (dsolve,numeric) which I suspect is the major culprit for the memory consumption. I have shortened the dsolve evaluation and managed to reduce total memory consumption to around 550MB.

Thank you guys for the suggestions, but doing 'a' or NULL doesn't seem to work, the memory just goes up. The discussion Alejandro suggested is interesting to read, and now I'm quite convinced that Maple doesn't handle memory well. It turns out that my loop uses memory up to 1.5G, barely acceptable by my PC, so I think I'll just brute force it. If in the future I need to run bigger loops I'll just break it to pieces, since each iteration in my current calculation is independent from one another. 

@acer 

Thank you for the further advice. I kind of get the feeling here, so functions and procedures seem to be different things in Maple. Functions are the elementary ones with analytic expressions, e.g. sin(x), exp(-t); yet procedures are more general than that. I looked up D on Maple's Help which says it can differentiate procedures, so that's why I couldn't differentiate my solution before, because I was using diff. So I guess I need to learn more about procedures then, because due to the nature of my study/work I almost always use non-elementary numeric procedures.

Thank you very much!

@Carl Love 

I see! So basically the good expression is only with the function's name, but without the (x) or (t) attached to it. No wonder everytime I tried to use:

source(t):=sin(0.1* t)  it will pop up a window asking what I'm trying to do with this, and I always select define a function...

The reason why originally I used f1(t):= eval(a1(t), sol1) and still able to plot might be due to the fact that the a1(t) here is a numeric solution of dsolve. I have noticed that the numeric solution of dsolve extracted by listprocedure is not a standard function -- it's not differentiable, and not reliably integrable (sometimes yes sometimes no depending on the upper bound of the definite integral, apparently having something to do with the maximum iterations allowed for the dsolve,numeric step). In my own case I had to use finite difference to appoximate the derivative and convert it to a Vector before I can Mean(Vector) it after using with(Statistics), hence got around the integration step.

If you have more tips on how to use individual solutions of dsolve, numeric better I'd greatly appreciate it!

If not I have already learned a lot from you anyways! Thank you so much!

@Carl Love 

Oh wow! It worked! Thank you! Could you please briefly explain why this is correct, if there is a reasoning at all? I asked because I have never defined functions this way, I always do as the following example:

source(t):=sin(0.1* t)

In fact it's wrong to define a function like the following:

source:=sin(0.1* t)

I never quite figured out when to have (t) after the function name and when not.

Another related question is that, now that I can use:

f1:= eval(a1(t), sol1)  and f1(100)

to extract numerical values, I did one step further to calculate the finite difference to approximate the derivative of a1(t):

d1(t):=f1(t)-f1(t-1)

d1(100)

And it seemed to work properly, but this time I had to use d1(t), not d1. However, when subsequently trying to plot it:

plot(d1(t), 1..2000)  returns Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

Turns out if I do:

plot(d1, 1..2000)  it works just fine

So this (t) is really confusing me now!

Thank you very much!

Thank you! This is very helpful, I think an event with a halt would do it nicely.

Thank you! This is very helpful, I think an event with a halt would do it nicely.

I got how to write events now, thank you! My original problem which was not well formulated couldn't be solved with this though, but this is very still enlightening. The thing is I may need a better conditional like "if" to regulate how the problem behaves. I have re-formulated in the following:  

 Main Goal: to solve a system of ODE that involves two unknowns dx(t)/dt=f(x,y)+sin(t) and dy(t)/dt=g(x,y) with the condition that x(t)>=0 and y(t)>=0 (initially x(t)=y(t)=0)

     Sub-objectives:

                                a. solve for x(t) and y(t) normally until one or both x(t) and y(t) are <0

                                b. when x(t) goes below zero solve a new system x(t)=0 and dy(t)/dt=g(x,y), while keeping track of f(x,y)+sin(t), f(x,y)+sin(t) will eventually go positive then we go back solving the original system

                                c. when y(t) goes below zero solve a new system dx(t)/dt=f(x,y)+sin(t) and y(t)=0, while keeping track of g(x,y), g(x,y) will eventually go positive then we go back solving the original system

                                d. if both x(t) and y(t) goes below zero we keep track of both f(x,y)+sin(t) and g(x,y), they will eventually go positive, then we go back solving the original system

I got how to write events now, thank you! My original problem which was not well formulated couldn't be solved with this though, but this is very still enlightening. The thing is I may need a better conditional like "if" to regulate how the problem behaves. I have re-formulated in the following:  

 Main Goal: to solve a system of ODE that involves two unknowns dx(t)/dt=f(x,y)+sin(t) and dy(t)/dt=g(x,y) with the condition that x(t)>=0 and y(t)>=0 (initially x(t)=y(t)=0)

     Sub-objectives:

                                a. solve for x(t) and y(t) normally until one or both x(t) and y(t) are <0

                                b. when x(t) goes below zero solve a new system x(t)=0 and dy(t)/dt=g(x,y), while keeping track of f(x,y)+sin(t), f(x,y)+sin(t) will eventually go positive then we go back solving the original system

                                c. when y(t) goes below zero solve a new system dx(t)/dt=f(x,y)+sin(t) and y(t)=0, while keeping track of g(x,y), g(x,y) will eventually go positive then we go back solving the original system

                                d. if both x(t) and y(t) goes below zero we keep track of both f(x,y)+sin(t) and g(x,y), they will eventually go positive, then we go back solving the original system

Thank you for the suggestion and I think it's very helpful. I read about events and couldn't figure out how to correctly write the event trigger. Naively I think it should be like 

events=[[x(t)>0, do_something],[y(t)>0, do_something]] 

but the program returns an error message: event1 trigger is of incorrect form: x(t)>0

Also I couldn't understand the way you wrote it either, you wrote 

events=[[x(t)-1,halt],[y(t)-2,halt]]

what does x(t)-1 mean for the event trigger?

Thank you for the suggestion and I think it's very helpful. I read about events and couldn't figure out how to correctly write the event trigger. Naively I think it should be like 

events=[[x(t)>0, do_something],[y(t)>0, do_something]] 

but the program returns an error message: event1 trigger is of incorrect form: x(t)>0

Also I couldn't understand the way you wrote it either, you wrote 

events=[[x(t)-1,halt],[y(t)-2,halt]]

what does x(t)-1 mean for the event trigger?