Here's a picture of it.

You seem to be using ln sqrt(3) instead of ln(sqrt(3)).  The logarithm function is a function and needs parentheses.  Otherwise Maple thinks you're multiplying.

In  that case it may be better to use the transformation method, which does not require derivatives of f.  The change of variables

t = x^(1/kappa)

takes

int(f(x)*x^(kappa-1),x=0..1)

to

1/kappa*int(f(t^(1/kappa)),t=0..1)

I think all Maple procedures that do file input/output, other than those involving Maple libraries, archives or helpo databases (which would use libname or savelibname) should use currentdir.

I imagine "Home directory" would mean the same thing to Maple that it means to the operating system: in Windows this would be from the environment variables HOMEDRIVE and HOMEPATH.

There are three standard tricks:

(1) break up the interval of integration into pieces, one near the singularity (x=0 in this case), one near infinity, and one where everything behaves nicely.

(2) transform an improper integral to a nice smooth one on a finite interval

(3) subtract off the singular part and integrate that exactly.

For example, to integrate x^(kappa-1)*f(x) for x from 0 to 1, where f(x) is an analytic function in a neighbourhood of this interval and 0 < kappa < 1, you can use the Taylor series f(x) = sum(c[i]*x^i, i=0..4) + g(x), where g(x) = O(x^5) as x -> 0.  Thus x^(kappa-1)*g(x) should be sufficiently nice (with 4 continuous derivatives) to behave well with Simpson's Rule, while int(x^(kappa-1+i), x=0..1) = 1/(kappa+i).
 

If what you want is (the first few terms of) a power series solution, use the option series in dsolve, after determining initial values at x=0.

The differential equation can be used to express higher derivatives in terms of p(x) and q(x).  You should be able to use rifsimp in the DEtools package for this.  See www.mapleprimes.com/forum/analytic-derivatives-of-numerical-solutions-of-an-ode-system

I think it goes like this.

  currentdir is the command used to check or change the current working directory, the one that procedures such as read and fprintf should use by default.  What this is initially will depend on how Maple is started.  For example, if you launch Maple by clicking on a worksheet file, currentdir() will return the directory containing that file.

  kernelopts(homedir)  is  the user's "home directory".  It is not
possible to change this within Maple.

  UserDirectory in maplesys.ini (on a Windows system) points to the directory where Maple would keep the file maple.ini (and .ini files for various users on a multi-user system).

 

This is a strange and rather alarming bug in dsolve(..., numeric).  The strangest part is that it hasn't been found and corrected long ago.  Here's a slightly simpler example.

> dsolve({diff(x(t),t)=-x(t)/t, x(1)=2}, numeric)(0.9);

The answer should be essentially

[t = .9, x(t) = 2.22222222222222232]

and it is, up to Maple 7.  But from Maple 8 up it is

[t = 0.9, x(t) = 1.90299400449141576]

I would guess that Maple is trying to make some sort of transformation, and ends up solving the wrong differential equation.

Although I don't see it documented in the help, it seems to me that legends in plots are for CURVES or POINTS objects and not for POLYGONS.  To fix it, you could make an outline of the polygon at the same time as the polygon itself.

> outlinedpolygon := proc(L::list)
      uses plots;
      display(polygonplot(args),
        pointplot([op(L),L[1]], args[2..-1], style=line))
end proc;

> outlinedpolygon([[0,0],[1,1],[1,0],[0,0]],
     colour=blue,legend="stuff");

Another way, perhaps slightly more elegant:

g := zip(()->args, e, f);

Or you can use

g := ListTools:-Interleave(e,f);

Tell Maple to do the integration:

> GeneralSolution := dsolve(SYS1,useint);

{C[10](z) = sin((-BR[1]-BR[2])^(1/2)*z)*_C4+cos((-BR[1]-BR[2])^(1/2)*z)*_C3+((BR[4]+BR[5])*C[50]+BR[3])/(BR[1]+BR[2]), C[51](z) = -1/2*(-2*BR[1]*_C2-2*BR[1]*_C1*z-BR[1]*z^2*Ped*u[0]*C[50*x]-BR[1]*z^2*BR[3]*C[50]-2*_C2*BR[2]+BR[5]*BR[3]*z^2+2*BR[5]*_C3*cos((-BR[1]-BR[2])^(1/2)*z)-z^2*Ped*u[0]*C[50*x]*BR[2]-2*_C1*z*BR[2]+2*BR[5]*_C4*sin((-BR[1]-BR[2])^(1/2)*z)-z^2*BR[3]*C[50]*BR[2]+C[50]*BR[5]^2*z^2+BR[5]*C[50]*BR[4]*z^2)/(BR[1]+BR[2])}

Express the solutions as functions:

> GS1 := map(t -> op(0,lhs(t)) = unapply(rhs(t), z), GeneralSolution);

GS1 := {C[10] = (z -> sin((-BR[1]-BR[2])^(1/2)*z)*_C4+cos((-BR[1]-BR[2])^(1/2)*z)*_C3+((BR[4]+BR[5])*C[50]+BR[3])/(BR[1]+BR[2])), C[51] = (z -> -1/2*(-2*BR[1]*_C2-2*BR[1]*_C1*z-BR[1]*z^2*Ped*u[0]*C[50*x]-BR[1]*z^2*BR[3]*C[50]-2*_C2*BR[2]+BR[5]*BR[3]*z^2+2*BR[5]*_C3*cos((-BR[1]-BR[2])^(1/2)*z)-z^2*Ped*u[0]*C[50*x]*BR[2]-2*_C1*z*BR[2]+2*BR[5]*_C4*sin((-BR[1]-BR[2])^(1/2)*z)-z^2*BR[3]*C[50]*BR[2]+C[50]*BR[5]^2*z^2+BR[5]*C[50]*BR[4]*z^2)/(BR[1]+BR[2]))}

The boundary conditions should be expressed using D rather than diff:

> BC1 := {(C[10])(h) = 1, D(C[10])(s) = 0, D(C[51])(h) = 0,D(C[51])(s) = Da[1]*(C[50]-J3)};

> eval(BC1,GS1);

{-1/2*(-2*BR[1]*h*BR[3]*C[50]-2*BR[1]*h*Ped*u[0]*C[50*x]-2*BR[1]*_C1-2*BR[5]*_C3*sin((-BR[1]-BR[2])^(1/2)*h)*(-BR[1]-BR[2])^(1/2)+2*BR[5]*C[50]*BR[4]*h+2*BR[5]*BR[3]*h-2*h*BR[3]*C[50]*BR[2]+2*BR[5]*_C4*cos((-BR[1]-BR[2])^(1/2)*h)*(-BR[1]-BR[2])^(1/2)+2*C[50]*BR[5]^2*h-2*h*Ped*u[0]*C[50*x]*BR[2]-2*BR[2]*_C1)/(BR[1]+BR[2]) = 0, -1/2*(-2*BR[1]*s*BR[3]*C[50]-2*BR[1]*s*Ped*u[0]*C[50*x]-2*BR[1]*_C1-2*BR[5]*_C3*sin((-BR[1]-BR[2])^(1/2)*s)*(-BR[1]-BR[2])^(1/2)+2*BR[5]*C[50]*BR[4]*s+2*BR[5]*BR[3]*s-2*s*BR[3]*C[50]*BR[2]+2*BR[5]*_C4*cos((-BR[1]-BR[2])^(1/2)*s)*(-BR[1]-BR[2])^(1/2)+2*C[50]*BR[5]^2*s-2*s*Ped*u[0]*C[50*x]*BR[2]-2*BR[2]*_C1)/(BR[1]+BR[2]) = Da[1]*(C[50]-J3), cos((-BR[1]-BR[2])^(1/2)*s)*(-BR[1]-BR[2])^(1/2)*_C4-sin((-BR[1]-BR[2])^(1/2)*s)*(-BR[1]-BR[2])^(1/2)*_C3 = 0, sin((-BR[1]-BR[2])^(1/2)*h)*_C4+cos((-BR[1]-BR[2])^(1/2)*h)*_C3+((BR[4]+BR[5])*C[50]+BR[3])/(BR[1]+BR[2]) = 1}

Now the next step should be to solve this system for the arbitrary constants _C1, _C2, _C3, _C4.  Unfortunately, this doesn't work, for a mathematical reason: your boundary conditions are not well posed.  One of the solutions of the homogeneous system is constant: C[10](z) = 0, C[51](z) = 1.  In other words, you can add a constant to C[51] in any solution of your differential equations and get another solution.  But your boundary conditions only involve the derivative of C[51](z), so they are unaffected by adding a constant to C[51](z).  Thus the boundary conditions can't determine the arbitrary constants.

Do you want to do this exactly (symbolically) or numerically?  For numerical optimization of a function of one variable, Optimization[Maximize] with the option method=branchandbound will attempt to find the global maximum in a given finite interval.  For example:

> with(Optimization):
  Maximize(sin(x^2/2) + sin(5*sqrt(Pi)*x/2), x = 0 .. 10, method=branchandbound);

[2., [x = 1.77245385093359986]]

A polynomial whose roots all have negative real parts is said to be stable. You can use DynamicSystems[RouthTable] to determine when that is the case: the polynomial can contain symbolic parameters, which are assumed real.  With the option stablecondition = true, it returns a Boolean expression which must be true in order for the polynomial to be stable.

> DynamicSystems[RouthTable](x^3 -c*x  + 16, x, stablecondition = true);

         false

On the other hand:

> DynamicSystems[RouthTable](x^3  + x^2 - c*x + b, x, stablecondition = true);

0 < b and 0 < -c-b

Some hints:

1) Your statements

f := x -> x;
a := a;  b := b; 

don't make sense.  f, a and b are formal parameters of your procedure.  That is, when somebody calls the procedure, e.g.

Approx(2, 3, x -> x^2, 5);

the statements in your procedure will be executed with a, b, f and N replaced by the values 2, 3, x -> x^2, and 5 that were specified.

2) All your procedure Approx has to do is use the formula for a right Riemann sum.
After your h := (b-a)/N, you just need one more statement.

3) You'll probably find sum useful.

There are two problems.  The first is premature evaluation.  Your code calls pLayer(x) where x is a symbolic variable.  Since (when x is symbolic) `mod`(16*x, 63) evaluates to 16*x, the result is piecewise(1 <= x and x <= 62, 16*x, `in`(x, {0, 63}, x), and that is what gets plotted.  One way to fix this problem would be to plot the function pLayer rather than the expression pLayer(x).  Thus

> plot(pLayer, 0 .. 63);

Unfortunately that won't work either.  The reason here is that `mod` doesn't work on floats, just on integers.  I'm not quite sure what you're really trying to plot, i.e. what you expect pLayer(x) to be when x is not an integer.  If you want to plot it for integer values of x, you could try

> plot([seq([x, pLayer(x)], x = 0 .. 63)], style=point);

The premature evaluation issue doesn't arise here because seq has special evaluation rules.