An interesting problem!

It can be solved in different ways. In my view, the most short way is to use a double integral with the change of variables. The change  u=y/x, v=(a-x)/y  maps the original region on the rectangle.

restart;

solve({u=y/x, v=(a-x)/y}, {x,y}):

assign(%):

int(Student[MultivariateCalculus][Jacobian]([x,y], [u,v], output = determinant), [u=1/2..1, v=1..3]) assuming a>0; 

                                                            7/120*a^2

A := plots[implicitplot](max(2-r, r-5, 3*Pi*(1/4)-theta, theta-5*Pi*(1/4)) = 0, r = 0 .. 6, theta = 0 .. 2*Pi, coords = polar, axiscoordinates = polar, gridrefine = 3):

B := plottools[polygon]([[2*cos(3*Pi*(1/4)), 2*sin(3*Pi*(1/4))], [5*cos(3*Pi*(1/4)), 5*sin(3*Pi*(1/4))], seq([5*cos(3*Pi*(1/4)+(1/200)*Pi*i), 5*sin(3*Pi*(1/4)+(1/200)*Pi*i)], i = 1 .. 100), [2*cos(5*Pi*(1/4)), 2*sin(5*Pi*(1/4))], seq([2*cos(5*Pi*(1/4)-(1/100)*Pi*i), 2*sin(5*Pi*(1/4)-(1/100)*Pi*i)], i = 1 .. 49)], color = green):

plots[display](A, B);

To automate the plotting of complicated plane figures, and to calculate their areas and perimeters, you can see my work

http://www.maplesoft.com/applications/view.aspx?SID=146470

Of cause, the problem can easily be solved by brute force method:

N:=0:

for n from 1 to 2013 do

if (irem(n,3)=0 and irem(n+1,4)=0) or (irem(n,4)=0 and irem(n+1,3)=0) then

N:=N+1; fi;

od:

N; 

                     336

In fact, the problem can be solved analytically for any range. The decision is based on the following arguments:

1. If the number  n  is divisible by 3, and the number  n+1  is divisible by 4, then the general formula for all such numbers is obtained as the solution of Diophantine equation  3*k+1=4*m . Similarly, if  n  is divisible by 4, and  n+1  is divisible by 3.

2. If  a .. b  is any  real range (a<=b), the number of integer points in this range is  floor(b) - ceil(a)+1 .

P:=proc(N1, N2)  # N1 and N2 specify the range N1..N2

local sol1, sol2;

isolve(3*k+1=4*m);

sol1:=solve({3*rhs(%[1])>=N1, 3*rhs(%[1])<=N2-1});

isolve(4*k+1=3*m);

sol2:=solve({4*rhs(%[1])>=N1, 4*rhs(%[1])<=N2-1});

floor(rhs(sol1[2]))-ceil(lhs(sol1[1]))+floor(rhs(sol2[2]))-ceil(lhs(sol2[1]))+2;

end proc: 

Examples:

P(1, 2014);

P(10^10, 10^20);

                      336

      16666666665000000000

If you want to use a symbol  V , and that there was no contradiction with the already defined vector  V , you can write 

V:=Vector(5):

for i to 20 do 
...
V||i =...
...
od:

Example:

V:=Vector(5);

for i to 20 do

V||i:=i^2;

od:

V||9;

There are infinitely many solutions of the form  a=i^2,  b=2*i*j,  c=j^2  (i, j  are integers),  that is the consequence of the identity

i^2*4^n+2*i*j*6^n+j^2*9^n = (i*2^n+j*3^n)^2

It remains to prove that there are no other solutions.

a := [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:

b := convert([ 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011, 0.011], fraction):

c := convert([ -0.88, -8.87, -0.86, -0.82, -0.77, -0.71, -0.66, -0.62, -0.57, -0.54, -0.89, -0.88, -0.85, -0.81, -0.76, -0.71, -0.66, -0.61, -0.57, -0.53], fraction):

P1:=CurveFitting[PolynomialInterpolation](a[1..10], c[1..10], x):

P2:=CurveFitting[PolynomialInterpolation](a[11..20], c[11..20], x):

F:=unapply(expand(P1*(y-b[11])/(b[1]-b[11])+P2*(y-b[1])/(b[11]-b[1])), x, y); 

The polynomial  F  is exact on your data:

is([seq(F(a[i],b[i]), i=1..20)] = c);

                          true

Apply a function  tan  to both sides of the equation and then after of simplifications  use   isolve  command. 

The only answer   {x = 0, y = 4}

From conditions of the problem all the angles are easy to find, and from the law of sines  all sides of the triangle can be expressed by  a :

is(cos(Pi/7)^2+cos(2*Pi/7)^2+cos(4*Pi/7)^2=5/4);

is(sin(Pi/7)/sin(2*Pi/7)/a+sin(Pi/7)/sin(4*Pi/7)/a=1/a);

                                               true

                                               true

You can find out the nature of the roots of a cubic equation, without solving the equation itself, through its coefficients in terms of the discriminant and the resultant (Delta  and  Delta[0]). See  http://en.wikipedia.org/wiki/Cubic_function 

This is the bug in Maple 16.00 . Use other versions or build individually in any version as follows (Carl Love's example):

A := plots[polarplot](cos(t), t = -(1/2)*Pi .. (1/2)*Pi, color = red):

B := plots[polarplot](sin(t), t = 0 .. Pi, color = blue):

plots[display](A, B, thickness = 2);

solve({sin(2*x) = 1/2, x <= 2*Pi, -2*Pi <= x}, AllSolutions, explicit);

I do not know how to use  events  option to your problem, but another way is possible. Yourselves set the parameters  h  and  epsilon . Then decide your equation by  dsolve({de, x(0)=1}, numeric)  (rkf45 method will by default)  and then in a loop by changing the argument with step =  h  interrupt the loop when abs( x(t)- x(t-h) )< epsilon .

Example:

restart;

de := diff(x(t), t) = -x(t):

ds := dsolve({de, x(0) = 1}, numeric):

h := 10^(-3):  epsilon := 10^(-6):  a := 0.:  b := 1.:

for i do

a := a+h: c := rhs(ds(a)[2]):

if abs(c-b) < epsilon then break end if:

b := c:

end do:

op(ds(a));

plots[odeplot](ds, 0 .. a);

I have not found a command in Maple that extracts the coefficient of the specified monomial  in a multivariate  polynomial and so I wrote a simple procedure that solves this problem. The procedure called  coefff returns the coefficient of the monomial  t  in the polynomial  P (T - the set of variables of the polynomial).

coefff:=proc(P, T, t)
local L, H, i, k:

L:=[coeffs(P, T, 'h')]:  H:=[h]:  k:=0:

for i from 1 to nops(H) do

if H[i]=t then k:=L[i] fi:

od:

k;

end proc:

Example:

P:=x^2-y^2+x-100:

coefff(P,{x,y},x*y),  coefff(P,{x,y},y^2),  coefff(P,{x,y},1);

                                      0, -1, -100

For specified values of the all parameters solutions can be found without any problems:

a1, b1, a2, b2, a3, b3 := 1.0 $ 6:

dC1 := diff(C1(t), t) = -a1*C1(t)+b1*C2(t):

dC2 := diff(C2(t), t) = a1*C1(t)-(b1+a2)*C2(t)+b2*C3(t):

dC3 := diff(C3(t), t) = a2*C2(t)-(b2+a3)*C3(t)+b3*O1(t):

dO1 := diff(O1(t), t) = a3*C3(t)-b3*O1(t):

syst := dC1, dC2, dC3, dO1:

ics := C1(0) = 1.0, C2(0) = 0., C3(0) = 0., C1(t)+C2(t)+C3(t)+O1(t) = 1.0: funct := [C1(t), C2(t), C3(t), O1(t)]:

sol := dsolve([syst, ics], funct);

plot(eval(funct, sol), t = 0 .. 5, color = [red, blue, green, yellow], legend = [C1(t), C2(t), C3(t), O1(t)]);

All works if some corrections have been done:

C:=proc(n,N)

binomial(N, n);

end proc:

theo:=proc(a1::symbol, a2::symbol, N1::integer, N2::integer)

local i;

add(add(d!*C(i,N1)*C(d-i,N2)/a1^i/a2^(d-i), i=0..min(N1,d)), d=0..N1+N2);

end proc:

Example:

theo(a,b,3,4);

The reason for the error is not clear.