Surely you want vl and vg to be positive. I suspect you also want them > bd, as otherwise you might have trouble with the 1-(1/v)*bd in the denominator, and vg > vl. You can use int rather than Int for your integral, if you tell Maple to make those assumptions:

> assume(vl > bd, vg > vl);
  constraint_1:= value(constraint_1);

constraint_1 := -.1000000000 ``(20.*vg+9.*vg*ln(2.*vl-1.)*vl-20.*vl-9.*vl*ln(2.*vg-1.)*vg)/vl/vg `` = (vl-vg)*(.9/vg/(1-.5/vg)-2/vg^2) Some experimenting with plots shows there is a solution near vl = 0.77, vg = 1.76:

> plots[implicitplot]([constraint_1,constraint_2], vl=0 .. 2, 
    vg = 0 .. 2,colour=[red,blue]);

And then you can tell fsolve to look for a solution near there:

> fsolve({constraint_1,constraint_2}, {vl = 0.7 .. 0.9, 
      vg = 1.7 .. 1.9});

{vg = 1.775178294, vl = .7447957481}

The output of C (and other CodeGeneration routines) is normally printed in the worksheet as a side-effect; the command itself returns NULL. The CodeGeneration Options help page says you can use the output option to send this output to a file or to a string. Thus:

C(bb, optimize, output="trygenC.cpp");

ought to work. Well, I don't know if it will work under Matlab, but in normal Maple it would work.

The -T option (and all others) are described in the help page ?maple. If you're starting Maple by clicking on an icon (e.g. in Windows XP), you can change the command line by right-clicking the icon, selecting Properties, and editing the Target line. But I doubt that changing this will help: you're more likely to be running into a hardware limit than a software limit, trying to allocate more memory than your system can handle. What might help is not running other applications concurrently, using Classic or command-line Maple rather than Standard, perhaps changing some of the settings in your Control Panel to free up more memory for programs. But the first thing you should ask yourself is whether the calculation you're trying to do is one that it's reasonable to expect Maple to do. It's all too easy to ask Maple to do something that would require more memory than is available in the whole world.

There is an unfortunate ambiguity in terminology here, and it is not always clear which version of covariance, variance, or standard deviation is meant. But it would be a good idea for Covariance, Variance and StandardDeviation to be consistent here. In particular, Covariance(A,A) should be the same as Variance(A). But this is not the case in Maple for data. As far as I can tell, for a data set (say a list), Variance and StandardDeviation use the definition with N-1, while Covariance uses N. For example:

> with(Statistics): A:= [1,2,3]:
  StandardDeviation(A),Variance(A),Covariance(A,A);

1., 1., .6666666667

I suspect the < problem struck again. I wish this could be fixed in some way: at least have an option to disable html tags (and that option should probably be the default). Anyway, to produce < you need to type < ; (without the space).

What the [1..7,1..8] does is take the submatrix consisting of the top 7 rows of A - tI (where t is left as a variable). Call that matrix B(t). Since B(t) has rank 7, its null space is spanned by a single vector v(t) (with entries rational functions of t) for which B(t)*v(t) = 0. My claim is that "generically", if lambda is an eigenvalue of A, v(lambda) should be an eigenvector of A for eigenvalue lambda. The reasoning is as follows: 1) The denominators of the entries of v(t) may be 0 for some values of t, but "generically" these will not be eigenvalues of A, because they don't depend on the last row of A while the eigenvalues do. So v(lambda) is well-defined and B(lambda)*v(lambda) = 0. 2) Similarly, the numerators may be 0 for some values of t, but "generically" not all will be 0 for eigenvalues of A, so v(lambda) is not the 0 vector. 3) Unless A has a left eigenvector whose last entry is 0 (which "generically" won't happen), the first 7 rows of B(lambda) will be linearly independent. Thus the null space of B(lambda) has dimension 1, and is spanned by v(lambda). But an eigenvector of A for eigenvalue lambda will also be in the null space of B(lambda), so v(lambda) must be a scalar multiple of it (and thus an eigenvector itself).

The reason the posted code didn't come out is that < is interpreted as the start of an HTML tag. You have to use < ; (without the space). BTW it would be better in general to use

  if type(x,numeric)

rather than

  if is(x,numeric)

but that doesn't make any difference in this case. This should work:

> evalf(Limit(k(x),x=1,left));

Unfortunately it doesn't, because of a bug: evalf(Limit(...)) doesn't work for something that is constant as the limiting point is approached. For example:

> evalf(Limit(1, x=0));

Limit(1,x = 0)

If that works at all, then u is not what Maple calls a vector. It must be a table (that's what you get if you just assigned values to some u[i] without first defining u). In that case a better way would be

> N:= max(op(map(op,[indices(u)])));

But just in case you want this to work on things that could also be vectors or Vectors or lists, try this:

> maxindex:= proc(u)
  if type(u,list) then nops(u)
  elif type(u,vector) then linalg[vectdim](u)
  elif type(u,table) then max(op(map(op,[indices(u)])))
  elif type(u,Vector) then LinearAlgebra[Dimension](u)
  else error "input is a %1",whattype(u)
  end if
end proc;
> maxindex(u); 

It does look complex, but turns out to be real (and also correct). Try this:

> S:= sum(binomial(2*k,k),k=1..n);
  seq(add(binomial(2*k,k),k=1..n) = simplify(S),n=0..5);

0 = 0, 2 = 2, 8 = 8, 28 = 28, 98 = 98, 350 = 350 The point is that hypergeom([1,n+3/2],[n+2],x) is being evaluated outside the radius of convergence of its series, using an analytic continuation which is complex at x=4, and that cancels the imaginary part.

In general, suppose you have a set S of linear inequalities in several variables, and you want to know whether there is a point that satisfies them all. If they were not strict inequalities, you could use simplex[feasible]. For strict inequalities, what you can try is this. Add another variable eps, with eps >= 0 and eps <= 1, replace each constraint A < B by A <= B - eps, and then try maximizing eps subject to these constraints. The answer is yes if and only if there is an optimal solution with objective value > 0. The following procedure should do it:

> myfeas:= proc(s::set({`<`,`<=`,`=`}))
   local C, strict, nonstrict, eps, R;
   strict,nonstrict:= selectremove(type,s,`<`);
   C:= map(t -> lhs(t) <= rhs(t)-eps, strict) union   
      nonstrict union {eps >= 0, eps <= 1};
   R:= simplex[maximize](eps,C);
   if R = {} then false
   else evalb(subs(R,eps) > 0)
   end if
 end proc;

In your example:

> myfeas({mU > 1-u, mU <1-u, mU > 0, u >0, u<1, mU <1, 
    mU >3/4*(1-u), mU < u,mU <1/2*u, mU < 1/4*u});

false

If you have the eigenvalues, you might do something like this. Suppose your 8 x 8 Matrix is A, and t is a variable that isn't one of your parameters.

> with(LinearAlgebra):
> NullSpace((A-t*IdentityMatrix(8))[1..7,1..8]);

"Generically" the result should contain a single vector v such that the first 7 entries of A.v - t*v are 0, and when t is an eigenvalue of A, v should be an eigenvector.

I suspect there's something wrong with what you're trying to plot. What does your "eq" look like? Perhaps you could upload the worksheet and data file.

If this is related to your previous question about an 8 by 8 jacobian matrix... I don't think there's any way to do this in parallel, but it may be possible to avoid completely. Presumably this must be a symbolic calculation, as a numerical one would be blindingly fast. The eigenvalues will be the roots of the characteristic polynomial, which has degree 8. I suspect you may have symbolic parameters, which could make that polynomial hellishly complicated (quite easy to have 8! different terms). The eigenvectors will then involve complicated rational functions of those eigenvalues. Perhaps the question to ask now is, what do you want to do with those eigenvectors? Do you really want them in explicit form?

The y values for this function corresponding to your x values in fs[1..200,1] could be calculated by

> yvals:= map(t -> eval(eq,x=t), fs[1..200,1]);

You could plot both your data and the curve from Fit by

> with(plots):
  xmin, xmax:= (min,max)(seq(fs[i,1],i=1..200));
  display([pointplot(fs), plot(eq,x=xmin..xmax)]);

Or perhaps better would be to plot the residuals, i.e. the pairs [x, y - eq].

> pointplot([seq([fs[i,1],fs[i,2]-yvals[i]],i=1..200)]);

A slightly better idea might have been to use the eigenvectors of the Jacobian matrix to get a starting point on the tangent plane to the stable manifold. I say "slightly better" because by the time you get out to k = 55/6, numerical errors may have destroyed any improvement that this makes. But anyway:

> F := [2*c^2+21/50*c-1/100*h^2-1/5*c*h-9/2000*h+49/4000,
         h*(c-1/5*h+1/5), 2/25*k-h-1/5*k*h];
  J:= eval(VectorCalculus[Jacobian](F,[c,h,k]),
         {c=1/50, h=11/10, k=-55/7});

J := RTABLE(157427492,MATRIX([[7/25, -61/2000, 0], [11/10, -11/50, 0], [0, 4/7, -7/50]]),Matrix)

> with(LinearAlgebra):
  E,V:= Eigenvectors(J);

E, V := RTABLE(169535840,MATRIX([[-7/50], [3/100+1/200*1158^(1/2)], [3/100-1/200*1158^(1/2)]]),Vector[column]), RTABLE(169625816,MATRIX([[0, -427/16000*(17/2+1/4*1158^(1/2))/(-25/4+1/8*1158^(1/2)), -427/16000*(17/2-1/4*1158^(1/2))/(-25/4-1/8*1158^(1/2))], [0, 119/400+7/800*1158^(1/2), 119/400-7/800*1158^(1/2)], [1, 1, 1]]),Matrix) Note that in this case the first and third eigenvalues are negative. So we want a linear combination of the first and third columns of V.

> p0:= <1/50,11/10,-55/7> + 1/10*(Column(V,3)-Column(V,1));
  sys:= {diff(c(t),t) = 2*(c(t))^2 + 21/50*c(t) -  
     1/100*(h(t))^2 - 1/5*c(t)*h(t) - (9/2000)*h(t) + 49/4000,
   diff(h(t),t) = h(t)*(c(t) - 1/5*(h(t)-1)),
   diff(k(t),t) = 8/100*k(t) - h(t) - 1/5*k(t)*h(t)};  
  Sol2:= dsolve(sys union {c(0)=p0[1],h(0)=p0[2],k(0)=p0[3]},
    {c(t),k(t),h(t)}, numeric,output=listprocedure);
  plots[odeplot](Sol2,[t,k(t)],t=-100..0,
     view=[-100..0,-55/7 .. 55/6]);

It does appear to hit k = 55/6 somewhere near t=-80.

> ksol:= subs(Sol2,k(t)): 
  t2:= fsolve(ksol(t)=55/6, t=-80 .. -60);

t2 := -76.88920001

> ic:= {c(0) = subs(Sol2,c(t))(t2),
   h(0) = subs(Sol2,h(t))(t2),
   k(0) = 55/6};

ic := {h(0) = .501371567495818016, c(0) = -.187351197287361840e-1, k(0) = 55/6}