@Carl Love Maybe I'm misunderstanding your intention, but the localized version still knows about the parameter 'a' and that can still be set:
 

restart;
Sol:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):
#
#"Carl's Magic .m" method:
Localize:= proc(Sol, params::list(complexcons))
option remember;
   Sol('parameters'= params);
   sscanf(sprintf("%m", eval(Sol)), "%m")[]:
end proc:
##
res:=Localize(Sol,[1]);
res(parameters);  #        [a = 1.]
res(1); #                  [x = 1., y(x) = 2.71828133411964]
res(parameters=[2]); #     [a = 2.]
res(1); #                  [x = 1., y(x) = 7.38905384345779]          

I thought that you wanted a procedure like the one you get when using dsolve/numeric for a concrete value of 'a', here a = 1.
In fact Sol and res above are identical:

res:=Localize(Sol,[1]);
interface(verboseproc=3);
S:=convert(eval(Sol),string):
R:=convert(eval(res),string):
StringTools:-Compare(S,R); # true
StringTools:-Compare(R,S); # true

The procedure(s) are terrifyingly long.
I checked that R (and so S) actually has the information about the parameter value by setting the parameter to 7.23456789 instead of 1.
Then I did

StringTools:-Search( "7.23456789", R );
R[1279..1306];   #  "Array(1..2, [1.,7.23456789])"

###
Furthermore, if you use the "naive" approach, but add option remember it seems that you get good results for the plot3d problem:
 

restart;

Sol:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):

Q:= proc(Sol, params::list(complexcons))
option remember;
   Sol('parameters'= params);
   Sol
end proc:

Y:= (A,X)-> 
   `if`([A,X]::[numeric$2], eval(:-y(:-x), Q(Sol, [A])(X)), 'procname'(args)):
gc():
P2:= CodeTools:-Usage(plot3d(Y(A,X), A= -2..2, X= -2..2, grid= [99$2])):

@Carl Love You can save the procedure as an .m file.
 

restart;
ode:=diff(x(t),t)=x(t)+a;
res:=dsolve({ode,x(0)=1},numeric,parameters=[a],output=listprocedure);
X:=subs(res,x(t));
res(parameters=[5]);
save res, "F:/MapleDiverse/temp.m";
seq(X(i),i=1..10);
res(parameters=[7]);
seq(X(i),i=1..10); 
read "F:/MapleDiverse/temp.m";
X:=subs(res,x(t));
seq(X(i),i=1..10); 

I should have added these lines at the bottom:
 

res(parameters=[7]);
seq(X(i),i=1..10);

Is your statement:

"It turns out that setting the parameter's numeric value in the
dsolve numeric solution-procedure causes the loss of previous details
of the numeric solving, even if the parameter's value is the same."

documented somewhere?
Or how to see that it is actually happening?

@siddikmaple Try this:
 

CSTR:={p[b], p[s], p[0, b], p[0, c], p[0, p], p[0, s], p[1, b], p[1, c], p[1, p], p[1, s], p[2, b], p[2, c], p[2, p], p[2, s], p[3, b], p[3, c], p[3, p], p[3, s], r[0], r[1], r[2], r[3], t[0], t[1], t[2], t[3]}=~(0..1);
ans := DirectSearch:-SolveEquations([e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21, e22, e23, e24, e25, e26],CSTR, evaluationlimit = 1000000);

The answer doesn't seem to satisfy the equations very well though. Is there a solution at all?
 

If fsolve returns unevaluted then it wasn't able to find a solution. It could help to give it a starting point or ranges . See the help for fsolve.
To get better help upload a worksheet with the equations. Use the green arrow in the MaplePrimes editor to do that.

@Zeineb I checked your copy of my code in Maple 18.02. It works.
See the reply to you from acer about getting Maple 18.02.
##
I just tried the code in Maple 15. I see the problem. expand doesn't do what it does in 18.02 and 2018.
This modification (where expand is replaced by IntegrationTools:-Expand) works in Maple 15:
 

restart;
A := Matrix(2, 2, [0, -1, 1, -1]);
X := t-> <x(t), y(t)>;
g := X(t)/(t+1)^2;
alpha := -1/4;
B := exp(-2*alpha*(u-t))*Matrix(2, 2, [0, 1, 1, 0]);
eq := diff~(X(t), t) = A.X(t)+int~(B.X(u), u = 0 .. t)+g;
eqs:=Equate(lhs(eq),rhs(eq));
## Now using IntegrationTools:-Expand instead of expand:
eqs:=IntegrationTools:-Expand(eqs);
S:=[y1(t)=int(exp(u/2)*y(u), u = 0 .. t),x1(t)=int(exp(u/2)*x(u), u = 0 .. t)];
sys1:=subs(rhs~(S)=~lhs~(S),eqs);
sys2:=diff~(S,t);
SYS:={sys1[],sys2[]};
## x(0)=1 and y(0)=1 are given; The last two follow from S.
ics:={x(0)=1,y(0)=1,x1(0)=0,y1(0)=0};
res:=dsolve(SYS union ics,numeric);
plots:-odeplot(res,[x(t),y(t)],-0.9..6);

From the image I see you are not using numerical solution. You need to do this here. For that you also need initial conditions.
Upload a worksheet using the green arrow in the MaplePreimes editor.

@Zeineb Could you upload a worksheet in which you have tried to implement my code?
Then I can have a look at it.

@Zeineb As Tom Leslie points out the code I gave works in Maple 18.02 as well. Elementwise operations were introduced in Maple 13. So the code should work in all versions from 13 and up.
Since you report the error message
" Error, (in dsolve/numeric/process_input) input system must be an ODE system, got independent variables {t, u} "

clearly you didn't just copy and paste the code into a fresh worksheet. The 'u' in "independent variables {t, u}"  must come from the integration variable used in your integrals.
The integrals should have disappeared totally in SYS.

restart;
dsolve(diff(y(x),x)=x);

Result:

  Warning, this version of the "Physics Updates" is intended

   for Maple release 2018.1, while you are using Maple 2018.2

  . Please check at http://www.maplesoft.com/products/maple/feat\

  ures/physicsresearch.aspx for the availability of a version

   of the "Physics Updates" for the copy of Maple you are using.
                              1  2      
                       y(x) = - x  + _C1
                              2         

The warning would be confusing to many. What has physics got to do with this ode one might ask.
I know the answer,  but I still find the warning unfortunate.

A consolation is that the warning only appears once after a restart. The next (or same) dsolve command doesn't result in a warning.

Followup:
I just downloaded Physics Updates from the Cloud (yes I know a link was given in the first warning).
The warning this time was:
"`Warning, this version of the "Physics Updates" is intended for Maple release ``2018.2`` and _libraryversion <= `1356656`, while you are using Maple release ``2018.2`` with _libraryversion = `1362973`. `

So maybe the earlier warning should have included a warning not to use the Cloud :-)

Note: I have now tried the link and as I thought downloaded and installed the Physics update:
Result: The same warning as from the Cloud.

@Lottie The first line in the procedure Q serves the purpose of making Q(A) return unevaluated if A is not an actual real constant (or doesn't evaluate to one).

if not A::realcons then return 'procname(_passed)' end if;

You could have written this instead:
 

if not type(A,realcons) then return 'Q(A)' end if;

If the name ttxx is unassigned then Q(ttxx) will return (literally) Q(ttxx).
The reason I wrote the line the way I did is simply that it's the way I usually do it.
It allows for inclusion of all arguments actually passed, but that is irrelevant here.

@Lottie That you could do to assure yourself of your results, but remember that your Lyapunov function V is just one that works, i.e. proves local asymptotic stability of (0,0).
The region in the plane where dV is negative doesn't necessarily correspond to the basin of attraction. In this case it certainly doesn't.

While I can't explain what I see from the lack of context, I hazard a guess that it may have to do with the fact that series remembers.
 

@ERA The OP wrote:
" I am working on an assignment and have some doubts about my solution to plotting the position of an anharmonic oscilliator with the Runge-Kutta method. "
Thus he might find it comforting that dsolve/numeric seems to agree with his results.

@ERA Surely you are right, but the OP said:
" I am working on an assignment and have some doubts about my solution to plotting the position of an anharmonic oscilliator with the Runge-Kutta method. "
Thus I thought it might be a comfort to him that dsolve/numeric appears to agree.