@MagnusIC3 Your assumptions on what is small and what is not seem to be quite ad hoc or arbitrary, and leave me feeling a bit unfomfortable.  Perhaps a more systematic approach can help.

Generally "smallness" makes sense in the context of a problem which is expressed in terms of a parameter, let's say h, and where we are interested in the behavior of the solutions as h approaches zero.  In fluid mechanics, for instance, h can be the Reynolds number, and the limiting case of h going to zero would lead from the nonlinear Navier-Stokes equations to the linear Stokes equations of motion.

The general method of linearizing a nonliner equation involving such a paramater h is to express the solution as a power series in h, substitute in the differential equation, and keep the lowest order terms.  So in your case, let's say

u(x,t) = 0*u[0](x,t) + h*u[1](x,t) + h^2*u[2](x,t) + ...
b(x,t) = 0*b[0](x,t) + h*b[1](x,t) + h^2*b[2](x,t) + ...

where I have written u and b for your uu0 and beta to simplify the notation.  Note that I have set the coefficients of u[0] and b[0] in the expansion above to zero because you said both u(x,t) and b(x,t) are small.

We substitute these in your PDEs, collect the like powers of h, and pick the coefficients of the lowest powers of h:

expansion :=
  u(x,t) = 0*u[0](x,t) + h*u[1](x,t) + h^2*u[2](x,t),
  b(x,t) = 0*b[0](x,t) + h*b[1](x,t) + h^2*b[2](x,t);

eval(rhs(KFBC)-lhs(KFBC), [expansion]):
expand(%):
collect(%, h):
coeff(%, h, 1) = 0;

eval(rhs(Eulerx[1])-lhs(Eulerx[1]), [expansion]):
expand(%):
collect(%, h):
coeff(%, h, 1) = 0;

This is a linear system of PDEs for the unknows b[1](x,t) and u[1](x,t). 

In a clarification to your question, you wrote:

   1D Boundary condition is A[0]=A[1] and A[M]=A[M-1] if I consider the grid points step is h=1/M.

Thus, your space-grid is indexed 0, 1, ..., M-1, M.  That is not good since the notation becomes messy.  The natural choice is to define M so that the space-grid is 0, 1, ..., M, M+1.  That way the indices 0 and M+1 will correspond to the boundary points, and the indices 1,2,...,M will be the interior points.  Note that after this change you will have h = 1/(M+1).

The next step is to express your equation in matrix form.  You should verify that in the case M=5 it takes the form:

Generalization to arbitrary M should be obvious.  Note the rightmost vectors on the first and the last lines.  They contain the values of A[] and a[] at the boundaries, that is, A[0], A[6], a[0], a[6].  Now, according to your boundary conditions, these may be replaced with A[1], A[5], a[1], a[5], respectively. The result is an equation that transforms the array a[] to the array A[].  Only the interior points, that is those of index 1,2,...,M, come into play.

The left hand side of this matrix equation is linear in A, therefore it may be written as P.A, where P is an MxM matrix.  The equation takes the form PA=Q, where Q is the mess on the right hand side, every bit of which is known.  So we compute A through A = P^(-1).Q.  This takes us from a[] to A[].  Then you may itereate, as you have noted, beginning with an initial condition.

The implementation in Maple is quite straightforward if you know a few things about the LinearAlgebra package.  Ask if you have questions about it.

Nota Bene: What I have sketched here is correct, however, it is not computationally efficient -- solving PA=Q by inverting P is quite wasteful.  Note that P is a tridiagonal matrix.  You will find efficient algorithms for solving a linear system with a triadiagonal coefficient matrix in most numerical analysis books.

First, fix the typo in the equation, as noted in Carl Love's response.  Then do the substitution as you have indicated.  You will find that your differential equation takes the form:

(...)*e + (...)*e^2 + ... + (...)*e^11 = 0

where I have written e for epsilon to simplify my typing.

In the small parameter method, one sets the parenthesized expressions to zero.  You will find that the first parenthesized expression is

Solving this we get:

The second parenthesized expression is

We substitute the previously computed x[1](t) here, then solve the resulting equation for x[2](t).  Then we combine the results by setting

This solution involves four integration constants _C1, _C2, _C3, _C4.  These are determined from the initial conditions as follows. Let's say the initial conditions are

where the coefficients a[1], a[2]. etc., are given.  Evaluatate x(t) and D(x)(t) at t=0 and match the coefficients of like powers of epsilon against those in the initial conditions to deternine the constants _C1, _C2, _C3, _C4.

We may obtain a pretty good solution through Euler's explicit method.  The following picuture shows the result.  The blue and red graphs correspond to step sizes of 0.1 and 0.01, respectively.  Decreasing the step size any further has no effect, therefore we can be confident that the red graph is the correct solution.

Here is the worksheet that produced that:  delay-equation.mw.

Note added later:

The solution shown above is nonsense.  In my worksheet I had a multiplication instead of a division by mistake.  After fixing the error we find that the solution exists on the interval [0, 0.84].  Beyond that point the "delay" becomes negative, and therefore the algorithm attempts to look up a not-yet-determined future value of y.  The solution ceases to exist beyond that point.  This is consistent with what was observed in Preben Alsholm's approach.

The solution on the interval of existence looks like this:

Here is the corrected worksheet: delay-equation2.mw

Let's look at a simple case first.  Suppose we want to expand the function f(c*x) about x, knowing that c is close to 1.  We have:

f(x+h) = f(x) + f'(x) * h + 1/2*f''(x) * h^2 + ...

We are interested in x + h = c*x, that is, h = c*x - x.  Therefore we get

f(c*x) = f(x) + f'(x) * (c*x - x) + 1/2*f''(x)*(c*x - x)^2 + ...

OK, enough practice. Turning now to your actual problem, we let:

xi := sqrt(m/(m+1))*(x+1/(2*sqrt(m))) - x;

eta := sqrt(m/(m+1))*(y-1/(2*sqrt(m))) - y;

mtaylor(f(x+a,y+b), [a,b], 2);

subs(a=xi, b=eta, %);

You don't need Maple to evaluate that integral.  You can tell the answer is 1/2 just but looking at the following picture:

The red curve is the graph of the function y = (1 - x^5)^(1/5).   The graph of the function y = x is shown in black.  The blue lines correspond to x = 2^(-1/5) and y = 2^(-1/5).  These divide the 1x1 square into six pieces of areas a, b, c, as indicated. Looking at the geometric meaning of your integral, we see that

• Your integral on the interval 0 < x < 2^(-1/5) calculates the area below the red and above the black curve, therefore it equals a + b.
• Your integral on the interval 2^(-1/5) < x < 1 calculates the area below the black and above the red curve, therefore it equals c.

Consequently, the entire integral equals a + b + c, which is clearly 1/2 of the area of the 1x1 square.

Remark 1: Note that the exponent 1/5 is immaterial to the argument.  Therefore, in general, for any positive p we have:

Remark 2: Come to think of it, the radical is immaterial as well.  The same argument shows that if the function f := [0,1] -> [0,1] is

1. continuous and monotone;
2. f(0)=1, f(1)=0;
3. the graph of y = f(x) is symmetric about the line y = x,
   

then

Addendum:

When I posted this solution yesterday, I had attempted to fill the various parts of the diagram with colors but I didn't know how.  After some search, I found a way of doing it in a post by tkolokol in MaplePrimes from 2009.  I made the following diagram with tkolokol's technique:

The diagram shows the unit square [0,1]x[0,1] in the xy plane as divided into four regions, the diagram being symmetric about the line y=x.  Think of the curve that connects the square's northwest and southeast corners as the graph of a function y = f(x).

Since the green and blue regions are congruent, then

    area of "yellow U green" = area of "yellow U blue" = 1/2.

But the area of "yellow U green" may be expressed as an integral. It follows that

The precise requirements on the function f are listed in Remark 2 above.

sin(-(1/6)*Pi+(1/2)*arccos(1/3));

expand(%^2):
combine(%, sincos):
sqrt(%);

In the code that you have shown, the intent is unclear.  On the one hand you say f is exp(x-t), therefore it is a scalar.  On the other hand, you use f as an operand in ScalarMultiply which is designed to multiply two vectors.  You need to resolve that inconsistency first.

You have:

divE := diff(e*x/(4*Pi*r^3), x)+diff(e*x/(4*Pi*r^3), y)+diff(e*x/(4*Pi*r^3),z);

It's likely that you meant:

divE := diff(e*x/(4*Pi*r^3), x)+diff(e*y/(4*Pi*r^3), y)+diff(e*z/(4*Pi*r^3),z);

Solve the first equation for theta(y) without any bounday conditions.  You will get the general solution in terms of Bessel functions and two arbirary constants C1 and C2.  Applying the condition theta(0)=1 eliminates C2 (because otherwise the solution will be undefined at y=0), and forces the value of C1 to be 1.  Applying the condition

then introduces a relationship between the various constants of your problem.  You will have to decide what to do with that relationship.  Neither Maple nor anyone else can help you with that because only you know what the equations are about.

Your worksheet needs a few changes before you can do what you have asked.

1. Remove all evalm().  They don't serve any purpose, and they complicate things.
2. Remove all asterisked variables.  Use b[i] instead of a[i]^`*`.

Then, after ATild has been defined, do:

ATildf := unapply(ATild, [a[0],a[1],a[2],b[1],b[2]]);

Now you may say:

ATildf(P,Q,R,S,T);

Here is the desired volume:

See the attached worksheet, volume-of-a-spherical-wedge-ver2.mw, for the derivation.

I have replaced the H₁ and H₂ with a and b to simplify the presentation.

You have

Eq:= (G-0.043)(Lc/2)+(-G)(Ls/2)+G*Lcon-0.043*Lc=14;

That should be

Eq:= (G-0.043)*(Lc/2)+(-G)*(Ls/2)+G*Lcon-0.043*Lc=14;

That will fix the problem and you will get a solution from fsolve().

However, your equation has two solutions.  There are at least two ways to get the second solution:

1.  You can tell fsolve() the interval to search for the second solution.  You can get a rough idea of where that solution lies by plotting the graph of Eq - 14.
   
2. Since your equation is particularly simple, you may find the two solutions at once with solve(), as in solve(Eq), instead of fsolve(Eq).

A side comment: The subject line of your inquiry refers to simultaneous equations.  But Eq is a single equation, there is nothing simultaneous about it.