I don't understand what you want to do.
Your worksheet contains the assignment
E := sum((V[i]-x[i]*((1+b)/(b*x[i]^4+1))^((1+b)/(4*b)))^2, i = 1 .. 78);
and nothing else.
What do you want to do with it?

Apparently the singularity is a problem.
Try
plots[contourplot](1/(x^2+y^2), x=-1..1, y=-1..1,contours=[seq(1..50,5)],grid=[100,100]);
plots[contourplot3d](1/(x^2+y^2), x=-1..1, y=-1..1,view=0..100);

You may try
?2DMathshortcutkeys

@Carl Love It appears that emmantop wants to (has to) write his own algorithm:

http://www.mapleprimes.com/questions/202416-How-Do-I-Write-Maple-Programme-For-Nonlinear-Ode

In that link the tag 'homework' appears. So I guess it is 'has to'.

Have you already got a program for the linear case? There the set of equations for the discrete version can be solved by LinearSolve immediately.
That is not so in the nonlinear case, where you will need to use Newton's method, but LinearSolve can still be used in each iteration.
The initial work in making a procedure handling the nonlinear case is basically the same as for the linear case. So making a procedure handling the linear case first makes sense to me.
I would make the procedure so it takes a system of first order as input.
 

@J4James In your test case this could be done like this:
restart;
eq1:=diff(f(y), y$4)-(diff(f(y), y$2));
bcs:=f(h1) = (1/2), f(h2) = -1/2, D(f)(h1) = -1, D(f)(h2) = -1:
h1:= 1+cos(x):
h2:=-1-cos(x+g):
db:=eq1,bcs:
d1 := subs(g=1,[db]);
F:=proc(xx,yy) local fp; global x,y;
   fp:=subs(dsolve(eval(d1,x=xx), numeric,output=listprocedure),f(y));
   fp(yy)
end proc;
   
F(0,.3); #Test of f(x,y)
plot3d(F,0..1,-2..2,grid=[20,20]);

##################
While the method above works it is unnecessarily slow: In doing the plot there will be n^2 calls to dsolve if the grid is [n,n]. The number of calls can be reduced to n if we change to:
Fx:=proc(xx) option remember; global y,x;
   if not type(xx,numeric) then return 'procname(_passed)' end if;
   subs(dsolve(eval(d1,x=xx), numeric,output=listprocedure),f(y))
end proc;
Fx(0)(.3); #Test
plot3d(Fx(x)(y),y=eval(h2..h1,g=1),x=-2..2);

 

@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)});