I followed the link from the Wikipedia page to the paper that has the Matlab code. Then I meticulously translated that into efficient and compact Maple code. The totality of the Matlab code accounts for only the GRF:= ...; paragraph of my code! The result is a greyscale image. To add color, I put the high values on a blue scale and the low values on a red+green scale.

restart:

Digits:= 15:
GRF:= proc(n,r)
uses AT= ArrayTools, IT= ImageTools;
local
     GRF:= (IT:-FitIntensity@IT:-Convolution)(
          AT:-RandomArray(n, 'distribution'= 'normal'),
          Matrix(2*r+1$2, (i,j)-> `if`((i-r-1)^2 + (j-r-1)^2 <= r^2, 1, 0))        
     )
;     
     (IT:-Preview@IT:-CombineLayers)(
          IT:-FitIntensity(map[evalhf](x-> `if`(x >= .5, 0, x), GRF))$2,
          IT:-FitIntensity(map[evalhf](x-> `if`(x < .5, 0, x), GRF))
     )
end proc:

GRF(300, 13);

When the number of parameters of a BVP is greater than or equal to the number of conditions at one of the boundary points, it can be solved as an IVP with fsolve being used to determine the parameters. I like this idea because I think that the main IVP solver (rkf45) is more robust than any of the BVP solvers. The following worksheet produces a very high-accuracy solution in good time using little memory, although I haven't compared it timewise or memorywise to Preben's. Accuracywise, my value of omega computed at Digits = 10 has three more digits of accuracy than Preben's. (I don't mean to denigrate his solution, for he just used the default error control.)

Download BVP_as_IVP.mw

If I understand you correctly, you want to increase the voltage, V, by steps and do fsolve at the increased voltages, at each step using as initial values the solutions from the previous step. Is that correct? In the worksheet below, I've done all this using a step of 0.1. The keys to making this work is to leave V symbolic and to use unapply to turn f1 and f2 into procedures with parameter V. The other changes that I made are just to keep everything neat and symbolic until it needs to made numeric.

Download 2DOF_mod.mw

There are various uses for ~ in Maple. Your syntax is wrong. For elementwise operation, the ~ goes after the procedure name. So that should be int8~(round~(S)).

Simply take your first example, and change type~(r,anything) to evalb~(r).

@tomleslie 

I believe that the following graph shows what not satisfying the boundary conditions smoothly means. This picks up from the end of the worksheet that you gave.

plots:-display(seq(plots:-odeplot(sol[k], [eta, f(eta)], 0..1, color= col[k]), k= 1..2));

I think that the critics may be complaining about that bump in the red curve. If you simply change the method from bvp[middefer] to bvp[midrich], this artifact will disappear completely. I think that this answers the OP's Question "Is there any other way?"

Correcting the abserr problem pointed out by Preben will also correct the anomaly in the red curve.

You just need to include two t-> and remove t= from the domain specification:

plot(
     [
          t-> add(BS(bb[4], t, i, 3)*bb[1](i+1), i= 0..n),
          t-> add(BS(bb[4], t, i, 3)*bb[2](i+1), i= 0..n),
          0..1
     ],
     color= red, axes= normal, scaling= constrained, numpoints= 100, thickness= 2
);

Note that the last item in the square brackets is 0..1 rather than t= 0..1.

The above can be made more efficient by including option remember at the beginning of procedures BS and bb.

If you want to create a list of indeterminate length and the value of each list element can be computed independently of the others (including that the stopping criterion doesn't depend on what is computed), then you should probably use a seq command. It is faster than a for loop. But there are a great many situations where the computations are not independent or the stopping criterion must be re-evaluated at each step. In those cases, you must use a for loop (which includes do loops and while loops). The best way to create a list in a for loop is this:

MyList:= Vector();
for n from 1 to 10 do
     MyList(n):= n+1;
end do;
MyList:= convert(MyList, list);

Note that this creates a list of size 10, not 9, which is what you said.

To reiterate what Tom Leslie said, a "list in {} form" is a set. Lists and sets have many similarities, but they're not the same thing. In particular a set isn't a type of list. To create a set in a for loop, do as above, but change the last line to

MyList:= convert(MyList, set);

If you have a group of expressions e1, ..., e9, they can be made into a set of equations equated to 0 like this:

MyEqns:= {e||(1..9)} =~ 0;

The reason for your print problems is a feature of the old-style matrices and arrays which isn't worth explaining here because you shouldn't be using them anyway. Your code can be vastly simplified and modernized to the new Matrices and Arrays like this:

restart:

Der:=proc(A,n)
local i,j,k,l,C;
     eval(
          Matrix((n,n), symbol= C),
          solve({
               seq(seq(seq(add(
                    A[i,j,k]*C[k,l] = A[k,j,l]*C[i,k]+A[i,k,l]*C[j,k],
                    k= 1..n), l= 1..n), j= i+1..n), i= 1..n-1)
          })
     )
end proc:

A1:= Array((1..2)$3, {(1,1,2)= 1}, storage= sparse):
Der(A1,2);

The above code for procedure Der is functionally equivalent to Tom Leslie's. If the Array A is highly sparse, as in the above example, there'll several occurences of the escaped local variable C. This may be a nuisance, and you may wish to replace it with a global variable. In that case, remove C from the local variables, and change the procedure's other four occurences of C to _C.

That is done by converting the first ODE's solution function for f(x) into something that'll return unevaluated when passed a symbolic argument and return numeric when passed a numeric argument and then declaring that as known in the dsolve command for the second ODE. Like this:

Sol_f:= dsolve({diff(f(x),x)=f(x), f(0)=1}, numeric):
F:= x-> `if`(x::complexcons, eval(f(':-x'), Sol_f(x)), 'procname'(x)):
Sol_g:= dsolve({diff(g(x),x)=F(x)*g(x), g(0)=1}, numeric, known= [F]);

Kitonum's Answer works also, but there are times when you want to solve the two as separate systems. And the above (using option known) will work for any numeric procedure, not just one that comes from a prior dsolve.

By the way, your title is perfect. I don't know what you see wrong with it.

That is surprising that it can't do it, because it can do it immediately if you replace the floats with symbols:

int(x*(1+a*x^4)^e*(C1-C2*x^4/(1+a*x^4)), x= 0..1);

     (1/8)*(4*C2*a^(5/4)*GAMMA(3/2+e)*GAMMA(1-e)*a^(-3/4+e)*hypergeom([1-e, -1/2-e], [1/2-e],                      -1/a)+4*C1*a^(1/4+e)*hypergeom([-e, -1/2-e], [1/2-e], -1/a) * a^(5/4)*GAMMA(3/2+e)*GAMMA(-e) * e-(4*        (1/2+e))*Pi^(3/2)*sec(Pi*e)*(C1*a*e+(1/2)*C2)) / (e*a^(3/2)*GAMMA(3/2+e)*(2*e+1)*GAMMA(-e))

evalf(eval(%, [C1= 2.906663106, a= 1.38, e= -0.4311594203, C2= 3.458929096]));

      1.00000000009749

This is the same (except for a little round-off error in the last digit) as what you'll get if you apply numeric integration (evalf(Int(...))) to your original integral.

If you want to compute an integral of a dsolve(..., numeric) solution, and one of the limits of integration is the initial point (or a boundary point) of the system, then the best way usually---the way that gives the best control of truncation error---is to build the integral directly into the system of ODEs, like this: Add to the system the ODE

diff(IntW(y), y) = W(y)

and the initial condition

IntW(ymin) = 0.

Then int(W(y), y= ymin..k) (which can't actually be directly computed like that) is equivalent to

eval(IntW(y), sol(k)); 

Edit: This idea comes from Joe Riel.

In the 2D input mode, when a variable is immediately followed by a left parenthesis, it is assumed that the variable is a function symbol being applied to the arguments inside the parentheses. In order to imply multiplication, you need to put a space between the variable and the parentheses. In other words, you need to change

p(4x^2+6) 

to

p (4x^2+6)

Oh, and with(student) serves no purpose in your worksheet. So why include it?

Use nested for loops, like this:

S:= 0:
for m from 1 to 5 do
     for j from 0 to m-1 do
          k:= m-1-j;
          S:= S + 2^m*(D@@j)(f)*(D@@k)(f)
     end do
end do;

where f(j,k) is the summand.

Or, do you really need k? This works:

add(add(2^m*(D@@j)(f)*(D@@(m-1-j))(f), j= 0..m-1), m= 1..5);

Yes, you are essentially right about the problem being due to extra brackets. The cause of the brackets is that XACT and YGPA are created as two-dimensional Arrays with a degenerate second dimension. Essentially, they're 120 x 1 matrices whereas you're probably thinking of them as vectors. XACT[9], for example, extracts the ninth row, not the ninth element. Those are essentially the same thing---the rows only have one element each---but to Maple your L1 looks like a list of pairs of 1x1 matrices rather than a list of pairs of numbers. We could convert XACT and YGPA to Vectors, and if I was going to work with them at length, I'd probably do that. But for this job that's not necessary. And they don't need to be converted to lists either. They can be plotted as they are, like this:

F:= Statistics:-Fit(a+b*x, XACT, YGPA, x);

Scatter:= plot(<XACT|YGPA>, style= point, symbol= circle, color= blue, symbolsize= 9):
Line:= plot(F, x= (min..max)(XACT), thickness= 2):
plots:-display(
     [Scatter, Line],
     title= "ACT vs GPA",
     labels= ["ACT", "GPA"],
     ytickmarks= [1,2,3,4],
     view= [14..35, 0.5..4.0]
);

The command to combine plots is plots:-display.

By the way, just out of curiosity, are those GPAs measured at high-school graduation, college graduation, or some other time?