@Mariusz Iwaniuk I'm not convinced. The method used is ftoc (try infolevel[int]:=5), but the antiderivative is discontinuous.
No indication that that is taken into account.

u:=convert(cos(2*x)/(1+2*sin(3*x)^2), exp);
U:=int(u,x);
plot(Re(U),x=-1..Pi+1);

@vv If we use the integral in the form int(-cos(2*x)/(-2+cos(6*x)), x ) instead, it appears by plotting (!) that that antiderivative is continuous. As pointed out in my answer Maple doesn't find the definite integral with the integrand written like that though.

@Carl Love Thanks, Carl, for your detailed explanation. I think I got it now. I'm slow.

@ecterrab The other day I deleted the Physics Update folder (named 2018) and after that the warning disappeared.
The point of that exercise was to find out if the warning would go away. Since it did, the warning will only come up for people actually downloading the Physics update. Thus they would presumably know that they have such a thing.
Thus I'm not too worried anymore.

Clearly, the hack above can be turned into a procedure.
Input the dsolve solution (res) and the number of odes of the corresponding first order system (ndiff).

restart;
FreezeParams:=proc(res,ndiff::posint,$) local R,params,ind,i,j,S;
  params:=rhs~(res('parameters'));
  ind:=[ListTools:-SearchAll('undefined',params)];
  R:=convert(eval(res),string);
  S:=NULL;
  j:=0;
  for i from 1 to nops(params) do
    if not member(i,ind) then 
      j:=j+1;
      S:=S,cat("Y",[ndiff+j])=convert(params[i],string)
    end if;
  end do;
  R:=StringTools:-Subs({S},R);
  parse(R)
end proc:
###  
## Example.
res:= dsolve({diff(y(x),x,x)=a*y(x)+b, y(0)=c,D(y)(0)=d}, numeric, parameters= [a,b,c,d]):
res0:=FreezeParams(res,2); 
res(parameters);
res0(parameters);
res(parameters=[a=1]);
res1:=FreezeParams(res,2); 
res0(parameters);
res1(parameters);
res(parameters=[b=2,c=3]); # Two more set
res2:=FreezeParams(res,2);
res2(parameters);
res2(parameters=[7,8,9,10]); #Attempting to set parameters
res2(parameters); # Not set
res0(parameters); #Unchanged
res1(parameters); # Unchanged
res(parameters=[d=4]); # All set.
res3:=FreezeParams(res,2);
result:=res(1); 
res3(1);
res(parameters=[0,0,0,0]);
res(1);
res3(1); # same as 'result'.

@Carl Love It seems that the solution procedure in dsolve/numeric uses Y[2] for the parameter in the case of one parameter in a system of one equation:
 

restart;

res:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):
res(parameters=[2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs("Y[2]"="2.",R):
res2:=parse(R2):
res2(1); # As result above
res2(parameters);  # Query
res2(parameters=[7]); # Seemingly set
res2(1); # But it isn't

Guessing the (obvious) generalisation:

restart;
### First order, two parameters:
res:= dsolve({diff(y(x),x)=a*y(x)+b, y(0)=1}, numeric, parameters= [a,b]):
res(parameters=[1,2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs({"Y[2]"="1.","Y[3]"="2."},R):
res2:=parse(R2):
res2(1);
res2(parameters); 
res2(parameters=[7,9]);
res2(1);
######################################
restart;
### Second order, two parameters:
res:= dsolve({diff(y(x),x,x)=a*y(x)+b, y(0)=1,D(y)(0)=0}, numeric, parameters= [a,b]):
res(parameters=[1,2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs({"Y[3]"="1.","Y[4]"="2."},R):
res2:=parse(R2):
res2(1);
res2(parameters); 
res2(parameters=[7,9]);
res2(1);

@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.