The Maple share library was replaced (as of Maple 6) by the
on-line Maple Application Center. However, I can't find
IntSolve there. In any case, I think IntSolve could
handle a single integral equation (of certain types) but
not a system of integral equations.
If you give us an example of the type of system you're
trying to solve, we might have some ideas of how to solve it.

You mean like this?
```
> imagefile := cat(kernelopts(datadir), "/help/ImageTools/fjords.jpg");
with(ImageTools):
img := Read(imagefile):
Preview(img);
Preview(Rotate(img,90));
```

You might try Import in the ExcelTools package (in Maple 11, Standard Worksheet). You can simply say e.g.
> A:= ExcelTools[Import]("c:/myfolder/myfile.xls");

I don't know why your h needs abs twice: since
abs(Omega*beta) >= 0, abs(Omega*beta)^(3/2) >= 0
as well.
Anyway, you can try something like this:
implicitplot3d(piecewise(h(beta,Omega) >0,
F(beta,Omega,delta),undefined),beta=-1..1,Omega=-20 .. 20,delta=-6..6,axes=box);

Try something like this:

sys:= {seq(D(x[i])(t)=x[i](t) - 2*x[5-i](t),i=1..4),
seq(x[i](0)=i,i=1..4)}:
sol:= dsolve(sys,[seq(x[i](t),i=1..4)],numeric,
output=listprocedure);
S:= subs(sol,add(x[i](t),i=1..4));
e:= 1.5:
plot(S,1 .. e);
fsolve(S(t) = t);

An error seems to be occurring somewhere in the statement
YP := ... in procedure p1. Actually you shouldn't
assign a value to YP at all, just to the entries
of YP. YP comes in to the procedure p1
as an array(1 .. N), and you're throwing that out and
replacing it by an array(1 .. N-1). And then
you try to assign a value to YP[N], which would not
exist if YP was an array(1 .. N-1). I don't understand
what you're trying to do in p1, but I suspect that
what you should be doing is something like
for i from 1 to N-1 do
YP[i] := ...
end do;
YP[N] := ...

Do you mean something like this?
plot([piecewise(t <=1,t,2-t^2), piecewise(t <=1,t^2,2-t),t=0..2]);

**roots** returns the roots of a polynomial in the rationals or an algebraic number field. In the case of x^2 + x + 1, you need your field to contain a square root of -3, so roots(x^2 + x + 1,I*sqrt(3)) would work. If you
don't know or care what field the roots should be in, you should probably stick with solve and fsolve.

As Jacques noted, if dsolve can't find a closed-form
solution it probably doesn't exist. You might be
able to use a numeric solution or a series solution.
See the help pages ?dsolve,numeric and ?dsolve,series.

Can you give us an example where this occurs?

Yes, this should be correct. Of course to actually
get "true" or "false" rather than an equation, you need to apply evalb to it.

To get explicit solutions for quartics, set
_EnvExplicit := true
before using solve. Alternatively, you can
try convert(..., radical) on the RootOf expressions.
Warning: these explicit results tend to be very
complicated. Personally, I find the RootOf notation
to be much nicer. There's really not much you can
do with the radical solution that you can't do with
the RootOf.

OK, I see the problem. This version should work in any version
of Maple from Maple 8 on:
PlotFeasibleRegion:= proc(S::{set(`<=`), list(`<=`)},
xr:: name = range, yr:: name = range)
# S is a set or list of non-strict inequalities
# xr and yr of the form x = a..b, y = c..d
# where x and y are the axis variables
# The usual plot options, such as colour=..., can be specified
# after these
local SS, pts, i, j, n, x, y, a, b, c, d, R;
x:= lhs(xr); y:= lhs(yr);
a:= op(1,rhs(xr)); b:= op(2,rhs(xr));
c:= op(1,rhs(yr)); d:= op(2,rhs(yr));
SS:= convert(S,set) union {x >= a, x <= b, y >= c, y <= d};
n:= nops(SS);
pts:= {seq(seq(solve({convert(SS[i], equality),
convert(SS[j], equality)}),
j=i+1 .. n), i=1..n-1)};
pts:= select(type,pts,set(name=realcons));
pts:= select(p -> andmap(is,subs(p,SS)), pts);
pts:= map(subs, pts, [x,y]);
R:=simplex[convexhull](pts);
if R = {} then return NULL end if;
plots[polygonplot](R,args[4..-1])
end proc;
Cheers,
Robert

For example, with some made-up data:
```
> X:= [$0..9];
Z:= [1.85, 1.54, 1.35, 1.23, 1.17, 1.16, 1.17, 1.22, 1.29, 1.38];
S:=dsolve((D@@2)(y)(x)=a*D(y)(x)+b*y(x)+c);
res:= [seq(eval(rhs(S),x=X[j])-Z[j],j=1..10)];
Optimization[LSSolve](res);
```

[.26471056267204e-4, [a = -.226946305957795996, b = .542887865054803298e-1, c = -.290795758320473961e-1, _C1 = 1.09634198339014442, _C2 = .217404308130228680]]
... and to see how well it fits:
```
> yfit:= subs(%[2],rhs(S));
P1:= plot(yfit, x = 0 .. 9);
P2:= plot([seq([X[i],Z[i]],i=1..10)],
style=point,symbol=circle, colour=black):
plots[display](P1,P2);
```

This is a case where the roots a/2 (+/-) sqrt(a^2+4b)/2
happen to be real. I don't think you'd be able to use the
same form in a case where they are complex: LSSolve
wouldn't like running into complex numbers. In that case
you'd want to use the other form of the solution, involving
sin and cos. In either case, you may also need to restrict
the ranges of a and b in the LSSolve command to avoid
complex numbers.
Cheers,
Robert

An alternative to Jacques's solution is to use **print** in your procedure to produce
any plots before the end of the procedure. For example
```
Energy:=proc(n)
.
.
.
print(plots[listdensityplot](...));
.
.
.
plot(...)
end proc;
```