@k20057 5
1. The Maple package Statistics has a procedure called NonlinearFit. I'm using the long name Statistics:-NonlinearFit for NonlinearFit. An alternative is to issue the command with(Statistics); after which you can use the name NonlinearFit instead.
2. NonlinearFit handles models like f(E)=B/(A*exp(E/kT)+1) containing unknown parameters, here B,A,kT. E is the variable. M is a matrix containing the experimental data with the E's in the first column and the values of the corresponding experimental f's in the second, here I shall refer to those as fexp(E).
3. NonlinearFit tries to find values for A,B,kT which minimizes the sum of squares of f(E)-fexp(E) over the given E's, i.e. add( (f(E[i])-fexp(E[i]) )^2, i=1..N) is attempted minimized.
4. The input to NonlinearFit is the sequence: Model (f(E) ), data (M), variable (E).

Edited: My remark about the value of A was due to a stupid typo.

Looking in a book (Kittel) on Solid State Physics I found a version of the distribution written like this:
f(E)=1/(exp(E-E_F)/(k*T)+1).
Written like that A would be given by A=exp(-E_F/(k*T)), thus positive.
So I wonder about the images you posted.

Maybe a normalization is needed (reflected in the parameter B):
f:=E->B/(A*exp(E/kT)+1);
Now you could either write down 3 equations using 3 points and solve for A,B,kT or (better) use all the points thereby getting an f that is best possible in some sense.
f:=E->B/(A*exp(E/kT)+1);
eq1:=f(0)=1.8;
eq2:=f(1)=1.6;
eq3:=f(2)=1.2;
solve({eq1,eq2,eq3},{A,kT,B});
plot(eval(f(E),%),E=0..9);
#The latter approach:
M:=Matrix([[0,1.8],[1,1.6],[2,1.2],[3,0.8],[4,0.4],[5,0.2],[6,0],[7,0],[8,0],[9,0]]);
res:=Statistics:-NonlinearFit(f(E),M,E);
plot([res,M],E=0..9,style=[line,point],symbolsize=25);

@blink23 Maybe ModuleLoad is what you need.
?ModuleLoad
So remove the protect command, make ModuleLoad a local and add these lines:
ModuleLoad:=proc() protect(i,j,k,epsilon); printf("Protecting i,j,k,epsilon\n") end proc;
ModuleLoad();

Will that be OK?

Just a comment: The boundary value problem as stated is solved without problem by dsolve.

@J4James Have you tried?
Take e.g. x=0.1. P2[0.1] is a difference between 2 procedures:
P2[0.1];
returns
proc(y) ... end proc-proc(y) ... end proc
Try plotting
plot(P2[0.1],-1-cos(0.1+1)..1+cos(0.1),view=-2..0);
You will see that P2[0.1] is constant as I said above, but it is a procedure in y.
You need the value of the constant. You could use any y in the interval, but the interval is specific to each x. y=0, however, is in the interval for all x.
When plotting a matrix using plot you need a matrix whose columns have numbers as elements, not procedures.
######################
Comment.
I have assumed that your equation is a simple example of something more complicated. If not you could easily solve the problem symbolically.
After having defined d1 and P1 do

res:=dsolve(d1);
res2:=eval(P1,res);
plot(res2,x=0..1);
Notice that the plot looks quite like the plot from the numeric approach.

@Alejandro Jakubi In ?solve,details it states under the heading 'Description':
"The solve command solves one or more equations or inequalities for the specified unknowns. The unknowns may be names, including indexed names (though for efficiency reasons, indexed names should be avoided when possible), or functions. Both indexed names and functions are considered to be independent of each other and of all other unknowns."

So functions ought to be OK as unknowns.

The following two versions do work, though:
p:=(x-3)/(9*x^5-48*x^4+73*x^3-6*x^2-36*x-8);
plot(p, x = -3 .. 3, y = -2 .. 2, discont = [usefdiscont=true]);
plot(unapply(p,x), -3 .. 3, y = -2 .. 2, discont = true);
and the command
discont(p,x);
works fine.
Also doing:
debug(discont);
plot(p, x = -3 .. 3, y = -2 .. 2, discont=true);
shows no sign of a problem although the plot doesn't conform to the findings of discont.

@Kitonum Having defined G already, you can just do:
H:=G@G;

@Carl Love Actually, the square brackets are unnecessary due to the special evaluation rules for evalf.
Notice the difference between the following two calls to evalf:
{x = 4/3+(1/3)*sqrt(7)}, {x = 4/3-(1/3)*sqrt(7)};
evalf(%);
evalf({x = 4/3+(1/3)*sqrt(7)}, {x = 4/3-(1/3)*sqrt(7)});

@mahmood180 If you want to keep the intervals, then don't simplify. You use fnormal and identify so you end up with a piecewise expression, which exactly simplies to the one in which only 0 and 1 are dividing points (only in the simple case in your second worksheet).
To see that clearly try this:
f:=eval(g1(t), sol);
nops(f);
[seq(`+`(op(i..2+i,f)),i=1..7,3)];#List with 3 elements
simplify~(%); #Simplifying each of the 3.
`+`(op(%)); #From list to sum
simplify(%);

@J4James I corrected an error in my reply on singularity. For f to have a fourth derivative (also) at r=-k you need f'''(-k)=0 (in addition to f'(-k)=1), not the equality of f'''(-k) and f''(-k).
This change doesn't really alter anything. The conclusion is the same.
You could try
dsolve(ode-1,f(r),formal_series); #You need the homogeneous eq. But f(r) is a particular solution to ode.
Notice that the two solutions are second degree polynomials and that they are not really different because of the arbitrary constants.
Such a solution will not solve your entire problem because it only has 2 arbitrary constants, so you won't be able to satisfy your 4 boundary conditions.
Another comment:
You can piece together solutions from the intervals -h..-k and -k..h like this:
ode as before.
bcs:=f(-h) = 1/2, D(f)(-h) = 1, f(h) = -1/2, D(f)(h) = 1:
h0:=evalf(1+0.5*sin(1)):
res:=dsolve(eval({ode,f(h) = -1/2, D(f)(h) = 1,f(-k)=h-k+1/2,D(f)(-k)=1},{k=0.5,h=h0}),numeric,method=bvp[midrich]); #Interval -k..h
plots:-odeplot(res,[seq([r,diff(f(r),[r$i])],i=0..3)],-0.5..h0,color=[red,blue,green,black],thickness=[1,1,2,2]); p1:=%:
#And choosing the straight line solution on -h..-k:
p2:=plot([r+h0+0.5,1,0,0],r=-h0..-0.5,legend=["f","f'","f''","f'''"],color=[red,blue,green,black],thickness=[1,1,2,2]):
plots:-display(p1,p2);
You notice that f and f' are continuous,but that f'' is not.
You may experiment to look for other solutions also satisfying the continuity conditions for f and f'.
I haven't found any other.

@mahmood180 I don't understand what you mean.
You could try
u0:=simplify((g2(t)-g1(t))^2); #Giving you an expression piecewise in t
u:=int(u0,t=0..2); #Giving you an expression in the c's
Now minimizing u gives the minimum of u as 0 and the values of the c's as 1.
That is easily checked independently:
eval(u,indets(u,name)=~1); #Returning zero
Clearly sqrt(u) then also has minimal value 0 at those values for the c's.

Try minimizing u:=int((g2(t)-g1(t))^2, t=0.. 2). Clearly that one is minimal where sqrt(u) is.

@J4James Using the simplification suggested by Thomas Richard:
Eq1:=simplify(Eq1,size);
#and then further
ode:=select(has,Eq1,diff);
##
We see that the value of B is irrelevant, but now the singularity at r=-k becomes more visible.
You are trying to solve a bvp on an interval -h..h which has this singularity as an interior point: In your case k=0.5 and h >= 1 (assuming that you always have your epsilon>=0).
No wonder numerical solution is difficult!
In my answer I had k=h=1, thus the singularity was at the left end point forcing me to use a midpoint method.

###
Notice that your equation doesn't involve f, but only the derivatives of f.
ode2:=subs(diff(f(r),r)=g(r),ode);
Boundary conditions g(-h)=1,g(h)=1. But you have the singularity at r=-k of course.
By doing
collect(ode2/(r+k)^3,diff);
you see that if g is going to be 3 times differentiable on -h..h then you need g(-k)=1.
(Edit begin) Also you need
u:=(k+r-1)*(diff(g(r), r, r))+(diff(g(r), r))+(-g(r)+g(-k))/(r+k); #Notice 1 replaced by g(-k)
to have the limit 0 for r->-k. This simply means that the second derivative of g must be zero at r=-k since
limit(u,r=-k);
returns -(D@@2)(g)(-k).
(Edit end)
The bvp for g on -h..-k has the constant solution g(r)=1. On the interval -k..h you have the same solution g(r)=1. Thus g(r)=1 on -h..h is a solution to ode2 so that f(r)=r+C would be a solution to ode. That solution will not solve both of the conditions f(-h)=-1/2 and f(h)=1/2. 
The conclusion is no solution f to the bvp can be 4 times differentiable at r=-k. So you have to find a way to deal with the singularity at r=-k.

@Aakanksha You get the same error in Maple 18. I really don't know anything about MeijerG, but could the parameter lists be inconsistent?
There is a procedure that checks the indices:

showstat(`MeijerG/check_indices`);

From looking at the link
http://en.wikipedia.org/wiki/Meijer_G-function

it appears that the requirements for the indices are that a[k]-b[j] must not be a positive integer for k =1..n, j=1..m.
In your case n=1 and a[1]=-15/2. m=4. So in your case a[1]-b[j] is not a positive integer for any j=1..m since
the b's are -5/2, 5/2, -15/2, -15/2.
Are there other requirements?

In Maple 18 I tried MeijerG on your parameter lists:
param:= [[-15/2], [-13/2]], [[-5/2, 5/2, -15/2, -15/2], []];
MeijerG(param,10.0); #returned .4656987417+0.*I
MeijerG(param,10); #gave the error you reported
So I'm suspecting a bug in the exact evaluation of MeijerG.
I submitted an SCR.