With your equations as written you can define d and a0(t) with piecewise:

ics:={a[1](0)=0.29, D(x)(0)=.1362, x(0)=0, theta(0)=10, D(theta)(0)=0.366, a[2](0)=0};
d:=piecewise(x(t)+a[1](t)>0,x(t)+a[1](t),0);
a0(t):=piecewise(x(t)+a[1](t)>0,-x(t),0)+piecewise(x(t)+a[1](t)>0,0,1)*a[1](t);
#Notice the somewhat complicated way a0[t] has been defined. That makes it possible for convertsys to solve for diff(a[1](t),t).
#The first attempt
#a0(t):=piecewise(x(t)+a[1](t)>0,-x(t),a[1](t));
#doesn't work
#Now you can continue with
res:=dsolve({eq1,eq2,eq3,eq4} union ics,numeric);
plots:-odeplot(res,[[t,theta(t)],[t,x(t)],[t,a[1](t)],[t,a[2](t)]],0..1,view=-10..15,color=[red,green,blue,black]);

If I understand you correctly this should be a solution:
eq:=f(n)*c(n+1) + g(n)*c(n) + h(n)*c(n-1)=0;
N:=2:
eqs:=eval([seq(eq,n=1..N)],c(N+1)=0);
LinearAlgebra:-GenerateMatrix(eqs,[seq(c(i),i=0..N)]);

From the help page:
"The is function returns FAIL if it cannot determine whether the property is always satisfied. This is a result of insufficient information or an inability to compute the logical derivation."

You are assuming that c<>0. Since f depends on a/c, b/c, and d/c only (besides x) you will either have no solution or infinitely many.
Try this.

f:=x->(a*x + b)/(c*x + d); #First version
(numer@(lhs-rhs))~([f(3)=1,f(4)=8,f(8)=-4,f(10)=-5]);
M,Z:=LinearAlgebra:-GenerateMatrix(%,[a,b,c,d]);
LinearAlgebra:-Rank(M);
#So you only got the zero solution which, however, is unacceptable.
#A more reasonable approach takes c = 1:
f:=x->(a*x + b)/(x + d);
solve([f(3)=1,f(4)=8,f(8)=-4,f(10)=-5],[a,b,d]); # No solution
(numer@(lhs-rhs))~([f(3)=1,f(4)=8,f(8)=-4,f(10)=-5]);
M,Z:=LinearAlgebra:-GenerateMatrix(%,[a,b,d]);
LinearAlgebra:-Rank(M); #Obviously not 4
LinearAlgebra:-Rank(<M|Z>); # Is 4

But clearly you must answer the question posed by Carl first. It is not clear how your solve-example is related to L.

Maybe you mean two triangles and the 4 points are the non-origin points?

A scalar version of the question was asked a few days ago and involved using the package DEtools.
http://www.mapleprimes.com/questions/152326-Differential-Operator#comment152362

Here is that method extended to the matrix case:
restart;
with(DEtools):
_Envdiffopdomain:=[Dx,x]:
A:=Matrix([[ mult(u(x),Dx)+v(x), a(x)], [b(x), mult(v(x),Dx)-u(x)]]);
 B:=< f(x),g(x)>;
C:=Vector(2,i->add(diffop2de(A[i,j],B[j]),j=1..2));

I used another method which was not nearly as elegant.
However, I shall repeat it here together with an extension to the matrix case:

M(p) will be the operator of multiplication by p (a function by default):

M:=proc(p,{constant::truefalse:=false}) local x;
   if constant then
      f->unapply(p*f(x),x)
   else
      f->unapply(p(x)*f(x),x)
   end if
end proc;
A:=Matrix([[ M(u)@D+M(v), M(a)], [M(b), M(v)@D-M(u)]]);
 B:=< f,g>;   
#`&p` will work as the noncommutative matrix multiplication but replaces multiplication by composition
`&p`:=proc(d,v) local a,b;
   a:=freeze~(d).freeze~(v);
   b:=thaw~(evalindets(a,`*`,apply@op));
   eval(b)
end proc;
C:=A &p B;
map(apply,C,y);
apply~(C,y);

I made a few but important changes:

modalemassaigvvoetgangers:=
unapply(piecewise(qverschilprocentueel>=5/100,modalemassametvoetgangers,
modalemassazondervoetgangers),hoogte);
qigvvoetgangers:=unapply(piecewise(qverschilprocentueel>=5/100,qtotaal,qaanwezig),hoogte);
eigenfrequentie:=unapply(22.4*sqrt(buigstijfheid*10^3/(qigvvoetgangers(hoogte)*10^2*lengte^4))/(2*Pi),hoogte);

reductiefactorvoetgangershoogte := proc (hoogte) local egf:=eigenfrequentie(hoogte);
  piecewise(egf <= 1.25,0,egf< 1.7,(egf-1.25)/(1.7-1.25),egf<= 2.1,1,egf < 2.3,(2.3-egf)/(2.3-2.1),egf <= 2.5,0,egf < 3.4,.25*(egf-2.5)/(3.4-2.5),egf <= 4.2,.25,egf < 4.6, .25*(4.6-egf)/(4.6-4.2), 0)
end proc;

Also reductiefactorvoetgangerseigenfrequentie could be rewritten using piecewise, but it works as it is.

Since
'2/0';
results in an error (by automatic simplification) I see no way of handling the problems as written.

However, 2/sin(Pi) evaluates to 2/0, but '2/sin(Pi)' is not automatically simplified, so there you could use try ... end try:

for i in [1, 2/'sin(Pi)', -3, 1/4, 5] do
  try
   if i::posint then print(i)  end if;
  catch:
  end try
end do;

try select(x->is(x::posint), {1, 2/'sin(Pi)', -3, 1/4, 5}) catch: end try;

If P is an Array (not array) then you could use ArrayTools:-FlipDimension:

P:=Array([[1,4],[2,3],[3,2],[4,1],[6,5],[6,1]]);
ArrayTools:-FlipDimension(P,1);

There is a missing multiplication sign between the two parentheses in the denominator:
The denominator should be
(x^2+y^2-4*x)*(x^2+y^2+4*x)
it is
(x^2+y^2-4*x)(x^2+y^2+4*x)
With this change you don't need to do anything else.
But you can when needed use the optional argument  view=0..16.

The operator of multiplication by a constant p in a function space could be defined by pop:
pop:=f->(x->p*f(x));
pop(g)(8);
(pop@@2)(g)(z);
((D-pop)@@2)(f)(x);
expand(%);

Since pop refers to p you could make a procedure M taking p as an argument and returning the operator of multiplication by p. You could e.g. do it like this:

M:=proc(p) local x; f->unapply(p*f(x),x) end proc;
((D-M(q))@@2)(f)(x);
expand(%);
#########
To handle also non-constant p you could do:
M:=proc(p,{constant::truefalse:=false}) local x;
   if constant then
      f->unapply(p*f(x),x)
   else
      f->unapply(p(x)*f(x),x)
   end if
end proc;

((D-M(q))@@2)(f)(x);
expand(%);

((D-M(q,constant))@@2)(f)(x);
expand(%);

If you insert a space after 'if' in  'if(v<=0)'  then it parses. 

Assuming that by poincaré coordinates you mean coordinates p1,p2,p3 related to x,z,w by
p1=x/(1+d), p2=z/(1+d), p3=w/(1+d), where d = sqrt(x^2+z^2+w^2) then you could do something like this:

f := (x, z, w) ->x^2+x*z-2*z+3;
h := (x, z, w) ->-3*x-3*z+z^2+x*z;
k := (x, z, w) -> w*x+w*z-3*w;
sys := diff(x(t), t) = f(x(t), z(t), w(t)), diff(z(t), t) = h(x(t), z(t), w(t)), diff(w(t), t) = k(x(t), z(t), w(t));
Equilibrium points in x,z,w:
pts:=solve({f(x, z, w) = 0, h(x, z, w) = 0, k(x, z, w) = 0}, explicit);
Poincare coordinates p1, p2, p3
pvar:=[seq(p[i](t),i=1..3)];
P:=sqrt(add(p[i](t)^2,i=1..3));
varchange:=[x(t),z(t),w(t)]=~(pvar/~(1-P)); #Variable change
PDEtools:-dchange(convert(varchange,set),{sys},[seq(p[i](t),i=1..3)]);
sysP:=solve(%,{seq(diff(p[i](t),t),i=1..3)}):
DEtools[DEplot3d](sysP,pvar,t=-10..10,[[p[1](0)=.1,p[2](0)=.2,p[3](0)=.3]],stepsize=.1);
#The right hand sides of sysP:
RHS:=evalindets(rhs~(sysP),function,f->op(0,f));
#The inverse change:
d:=sqrt(x^2+z^2+w^2);
invvarchange:=[p[1]=x/(1+d),p[2]=z/(1+d),p[3]=w/(1+d)]; #As given in my introduction
#Equilibrium points in p-space:
Ppts:=seq(eval(invvarchange,pts[i]),i=1..3);
simplify([Ppts]);
evalf(Ppts);
seq(simplify(eval(RHS,Ppts[i])),i=1..3); #Ought to give zeros and does

#In helping someone else he gave me a link to a paper mentioning poincaré coordinates on p. 18 of the pdf-version obtained from
http://arxiv.org/abs/gr-qc/0612180

If you can live with the numbering, this is somewhat short:
sys:={x[1]-x[2]=1, x[1]+x[2]-x[3]-x[4]=5};
res:=solve(sys);
S,R:=selectremove(evalb,res);
f:=indets(S)=~subs(x=t,indets(S));
f union subs(f,R);

You could use Heaviside:

restart;
f:=Heaviside(t-3)-Heaviside(t-4);
plot(f,t=0..5,thickness=3,discont=true);
eval(f,t=3);
?Heaviside
_EnvUseHeavisideAsUnitStep:=true:
eval(f,t=3);
plot(f,t=0..5,thickness=3,discont=true);
#The same could be done by using piecewise:

g:=piecewise(t>=3 and t<4,1,0);
plot(g,t=0..5,thickness=3,discont=true);

You can set the environmental variable UseHardwareFloats to false.
Simple example:

restart;
Digits:=8:
UseHardwareFloats:=false:
LinearAlgebra:-GaussianElimination(Matrix([[1.2345,6.789],[11.1213,-9876]]));

See
?LinearAlgebra,details
at the bottom there is a link to a worksheet with examples: Examples/LA_NAG