@wenny You should try executing the worksheet again. I did not get that error message.
If you get that error once again then try this:

indets({g1,g2,g3},name);

you should be seeing the unassigned names and hopefully what you see are just the variables.

@ I'm apalled to learn that some math professor out there cares a hoot about how you indicate the end of proof as long as you do it!

@jonlg Maple simply cannot do the integral even when T1, T2, U0, L0 are given concrete values.

I dropped the case t<=0 because I thought (perhaps wrongly) that you had no interest in t<0.

piecewise works like this: if the first condition is satisfied then the following expression is returned. If not then the next condition is evaluated, etc. If none of these is satisfied then 0 is returned. That 'otherwise' case is specified explicitly in my example as the last argument; when given explicitly it can be anything.

Could you give us a (simple) example of what you have in mind, thereby answering questions like 'Does the parameter have to be determined also or can it chosen arbitrarily' ?

@Konstantin@ A common (and recommended) method is to include the integrals in the system.

Suppose as a simple example that you want to find the integral int(z(x)*q(x),x=0..a), where q is a known function.
Then you could include the ode diff(w(x),x)=z(x)*q(x) with w(0)=0 in the system:

sys:={diff(z(x), x$2)+z(x) = -2, diff(w(x),x)=z(x)*q(x), z(0)=0, z(1)=1, w(0)=0}.

Then the matrix returned with the option=Array(...) will include w-values, and of course w(a) will be your integral.

Did you actually try this:

dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series',order=10);

It works for me in Maple 18, Maple 15, and Maple 12. I'm pretty confident that it works in all versions from 12 and up. Maybe it also works even in some earlier releases. I don't have access to those on this machine.

As for series:

series(cos(x),x=0,10);

@Markiyan Hirnyk The system DE1, DE2 is clearly NOT autonomous, since rhs(DE2) depends on t explicitly through cos(t).
When you tried using DEtools[autonomous] you had a wrong independent variable at the end.

DEtools[autonomous]({DE1, DE2}, {y, z}, t); #t not x

@Konstantin@ I should have written 1 and 3, not 1 and 2. Notice that I go on quoting those points.

The puzzling part about 1 is that read for me works on .m files in the sense that something is read. 'read' doesn't produce an output, it just reads stuff. In this case this means that S1 is assigned to the list you see, when you write S1; after having executed the read statement.
I thought that your original problem was that indeed something is read, but that something (S1) doesn't work as expected. It is difficult to read your original question any other way.

@wenny You should use the numerical solver, fsolve instead of solve. The integrals are inert integrals (Int instead of int). When fsolve gets to work it will evaluate the integrals numerically.
You will have more of a chance of success if you provide starting values or ranges.

@Konstantin@ I use Maple 18.02 on Windows 7. Which Maple version are you using?

I just had occasion to repeat the recipe I gave and with the same result as mentioned earlier.

I'm puzzled by your points 1 and 2:

that you cannot read the file you saved as an .m file
and
that when you save with any other extension everything works. (!!?)

I tried with 2 different extensions: .mpl and .txt. It didn't make the read-in value of S1 any better. But this time the reason was another, I got the error you report in your point 3. So do you really mean that everything works for you if you save with a different extension than .m ?

I tried your ode with initial values instead of boundary values. There a similar remedy I attempted with the procedure `dsolve/numeric/SC/IVPdcopy` didn't succeed. But see the remark about the range option at the bottom.

If you try

F1(0.12345);

you get the error:

Error, (in F1) argument #2, must be the target rtable

That may be because in F1 there is an Array named _dat whose second element contains pointers to places in the memory. These pointers are unique to the session. (I hope this is a somewhat accurate description).
Thus when saved in one session and read into another they will not point to anything (relevant).

showstat(F1);

The first element of _dat (out of 4) is a (large) procedure. The remaining were in my session:
2 = (Array(1..3, [4466537074,4466537218,4466537314])), 3 = [x, z(x), diff(z(x),x)], 4 = []

pointto(4466537074); #Error

If you try the same in the session where you save S1 you will see that the similar pointto command outputs a procedure. The 3 numbers refer to the procedures for x, z(x), diff(f(x),x) respectively.

If you take another look at the boundary value problem you will see that in F1 there is an Array called data quite like _dat. Its second element is an Array of three pointers (addresses).

That the copying device described above works and the pointers are apparently not needed after that is interesting.
#### The range argument in the ivp-problem
Trying with the range argument in the ivp problem showed that indeed a similar device worked here. Executing and copying in a fresh worksheet the procedure `dsolve/numeric/SC/IVPval` and assigning the pasted version to `dsolve/numeric/SC/IVPval` in our "read"-session was successful.

CONCLUSION: To avoid the trouble use output=Array(...) if you need to save data between sessions.

@wenny Certainly it is possible that there is a solution in some other range than {gam0=0..230, gam1=0.1..10}.

So do you have any guesses as to where a solution (if any exists) should lie?
Should gam0 and/or gam1 be positive?

The reason I tried this:

evalf(eval([g1,g2],{gam0=0.01,gam1=3.62}));

was simply to check that actual numbers came out.
If I couldn't find some gam0 and gam1 for which numbers were produced there would be no reason to try fsolve.

Does the problem come from some application (in physics or engineering e.g.)?

@tira I looked at your latest worksheet Exact_for_u.mw.

I corrected a few things, and did manage to get a graph in the end.
In two places a[0], which in your case is negative, appears under a square root sign. The one place is in u2, where it appears as (exp(y*sqrt(Pr*a[0])), the other place is in u7, where it appears as sqrt(-alpha[0]*a[0]).
In u2 the combined result is real (it appears). In u7, however, the result is that u7 is purely imaginary.
Thus I tried changing the sign under the square root as you will see below. Then u7 becomes real and thus also U since all the other are real.
However, there are two obvious problems: The values of u6,u7, and u8 are huge, and so is U. Also U is negative.

restart;
alpha[0] := 1/R; a[0] := (Pr-1)/(Pr*R); Gr0 := Gr/(Pr*R);
u1 := -erfc(y*sqrt(Pr)/(2*sqrt(.4)))-a[0]*((.4+(1/2)*y^2*Pr)*erfc(y*sqrt(Pr)/(2*sqrt(.4)))-y*sqrt(.4*Pr)*exp(-y^2*Pr/(4*.4))/sqrt(Pi));
u2 := (1/2)*exp(.4*a[0])*(exp(y*sqrt(Pr*a[0]))*erfc(y*sqrt(Pr)/(2*sqrt(.4))+sqrt(.4*a[0]))+exp(-y*sqrt(Pr*a[0]))*erfc(y*sqrt(Pr)/(2*sqrt(.4))-sqrt(.4*a[0])));
u3 := y*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u4 := .4*a[0]*y*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u5 := y*exp(.4*a[0])*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u6 := y*sqrt(alpha[0])*(Int(Int(exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
#u7 := y*sqrt(alpha[0]*a[0])*(Int(Int((t-s)*exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
u7 := y*sqrt(-alpha[0]*a[0])*(Int(Int((t-s)*exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
u8 := y*sqrt(alpha[0])*exp(.4*a[0])*(Int(Int(exp(-y^2/(4*R*u)+u+alpha[0]*s-a[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
Pr := .71; R := .2; Gr := .5;
evalf(eval([u1,u2,u3,u4,u5,u6],y=1.12345));
evalf(eval([u7,u8],{y=1.12345,t=1}));
U:=Gr0*(u1+u2+u3+u4-u5+u6+u7-u8)/a[0]^2;
evalf(eval(U,{y=.1,t=.1234}));
evalf(eval(U,{y=1,t=1.234}));
evalf(eval(U,{t=0,y=1}));
plot(eval(U,t=1),y=0..10);

@tira I took a look at the worksheet Exact.mw.
There are some syntax problems:

1. Many missing multiplication signs.
2. e^x is used a couple of times instead of exp(x), which is the correct syntax.
3. List brackets are used in the definition of u(y,t). They should be ordinary parentheses ().

A misprint: At the start is written alpha[o], later alpha[0].

Suggestions:

Just use the assignments  u1:= ... etc. No need for u1(y,t).
Use capital U in U:=Gr[0]*(u1+u2+u3+u4-u5+u6+u7-u8)/a[0]^2.

Most important of all: Use 1D input instead of 2D-input. It is so much easier to catch syntax errors when using 1D-math input (also known as Maple Input): In Maple go to Tools/Options/Display/Input display and choose Maple Notation. Then click the button Apply Globally (or Apply to Session).
This change applies to new worksheets. It won't change what you have written already.

Finally, there is no explanation of what is going on in the worksheet Exact.mw.. Is the claim that U solves the pde system and the boundary and initial conditions mentioned in one of your first worksheets?

@tira This is just a comment about your plots. It applies to either one of the worksheets.
I don't see the point in using seq(..., j=4), so this would be less opaque:
p[i] := pds:-plot(u, t = 4*/10, y = 0 .. 6, legend = [i]);

Another comment. You can make an animation with Gr as animation parameter like this:

Gr:='Gr': #Clearing Gr from the value assigned in the loop
pp:=proc(gr) local pds;
  pds:=pdsolve(eval(PDE,Gr=gr), IBC, numeric, spacestep = 1/100);
  pds:-plot(u, t = 4*/10, y = 0 .. 6,title=typeset("Gr = ",gr))
end proc;
plots:-animate(pp,[gr],gr=0..2,paraminfo=false,trace=24);

#You may want to leave out 'trace = 24'.