@Rouben Rostamian  By de2 v cannot be zero on any interval no matter what value is assigned to signum(0).
If you look at an interval where v appears to be zero as it does on 3.25..3.9:
odeplot(dsol, [t, v(t)], t=3.25..3.9);

then we see rapid oscillation. Also the following plot shows that.
odeplot(dsol,[t,signum(v(t))],t=3.25..3.9,thickness=3,style=point);

But as I understand the physics the house must be at relative rest (i.e. v(t) = 0) on intervals.
#########################
Making changes in order to match the original problem I tried using g = 9.81, mu=0.8, and cos instead of sin:

X := t -> P*cos(omega*t);
V := D(X);
A := D(V);
ic := x(0)=0, v(0)=0;
m:=1: g:=9.81: P:=1: mu:=0.8: omega:=1:

The results where disturbing. This one in particular:
odeplot(dsol, [t, v(t)], t=0..6*Pi/omega);

Another observation: While the results when using stepsize and abserr looked much like each other for your choice of parameters, they didn't for the choice above. The graphs of x(t) were rather different, while both graphs of v(t) were rather wild (but different).

Thus I don't trust the results.

And I should add: The mathematical model consisting of just de2 (and de1).

@phil76600 I tested Rouben's code in Maple 12.02. It didn't complain about 'stepsize'. Which Maple do you have 12.01 or 12.0?

But try replacing the word 'stepsize' by 'abserr'.

maxfun = 0 means that there is no limit to the number of function evaluations when computing the solution. You could say that a more intuitive way to do this would have been maxfun = infinity.

@Rouben Rostamian  This is very interesting. I shall have a closer look.

I was surprised by the effect of giving the stepsize. I didn't know that option existed except for the classical methods!

@phil76600 Your revised problem:

ode1 := m*diff(x(t), t, t) = k*Xr-mu*m*g if mu*m*g<k*Xr else diff(x(t),t) = V(t) = 0

(if I understand it correctly) can be solved by using 'events'.
See the help for that.
?dsolve,events

The event [k*Xr=mu*m*g..infinity,diff(x(t),t)=0] means that if kXr falls outside of the interval mu*m*g..infinity at some time t1 then the action is to set diff(x(t),t) = 0 at that time t1. 

restart;
Xr := U*t-x(t);
k := 10; m := 1; g := 10; mu := .2; U:=1;
ode1 := m*diff(x(t), t, t) = k*Xr-mu*m*g;
res:=dsolve({ode1,x(0)=0,D(x)(0)=0},numeric,events=[[k*Xr=mu*m*g..infinity,diff(x(t),t)=0]]);
plots:-odeplot(res,[t,x(t)],0..15);
plots:-odeplot(res,[t,diff(x(t),t)],0..15);
plots:-odeplot(res,[t,mu*m*g-Xr],0..15,thickness=3);

@phil76600 Let us not mix the two discussions, but let me just point out that mu*m*g>m*ω^2*A only depends on the given parameters, so either this relation holds under the whole process or it doesn't.

@phil76600 Actually they are not the same. Take a look at the of the interval t=0..15 (e.g.).
I saved both plots and displayed them together:

Notice that for mu = 1000 the piecewise gives 0 for all t in the interval 0..15. This is not the case for mu=0.2.
Using
pw:=piecewise(mu*m*g<k*Xr,mu*m*g,0);
you can try
plots:-odeplot(res,[t,pw],0..15,thickness=3);
You will notice that for mu = 1000 the values are all zero (therefore I added thickness=3). For mu=0.2 this is definitely not the case.

@Rouben Rostamian  Using your notation a = diff(v(t),t), so from your equation de2 we get a+A = -mu*g*signum(v(t)). Thus abs(A+a) = mu*g*abs(signum(t)), which is equal to mu*g if v(t)<>0 and equal to 0 if v(t)=0 if we accept that signum(0) = 0.
Thus the first part of the piecewise right hand side only holds if v(t) = 0, whereas the second case (the otherwise,output 0) holds whenever v(t)<>0.
Surely, you didn't mean it like that.

@tobbie247 If you my code from above and try
res(0);

you will get very nearly the values for A and B you gave in a reply to Carl:
0.742408087440 and 0.943866213084

I don't know anything about the Differential Transform Method either.

@phil76600
I'm not happy with my solution either. I'll have to rethink this.

## Why do you have cos(t)^2 in the Maple code when cos(omega*t) is not squared in the image of the formula?
## Are you sure mu = 0.8 is correct?
## Is the model actually correct as stated (apart from the square mentioned above?)
As I understand it presently, the ground on which the house is standing is moving in a periodic way as in
xsol(t) = A*cos(omega*t).
x(t) is the position of the house relative to the ground.
Thus the abolute position of the house is X(t) = x(t) + xsol(t).
Using Newton's 2nd law we get
m*X'' = friction
so that
m*(x'' + xsol'') = friction
Thus
x'' = friction/m - xsol'' = friction/m + A*omega^2*cos(omega*t) .
The problem for me is what to set friction to be when there is no relative motion, i.e. when v(t) = 0.
Am I misunderstanding something?

@tobbie247 You can try dsolve as in the following, where I had success with infinity = inf = 28, but not with inf = 29.
For my convenience I used x instead of eta.

ode1:=diff(f(x),x$3)+3*f(x)*diff(f(x),x,x)-2*diff(f(x),x)^2+theta(x)-m*diff(f(x),x)=0;
ode2:=diff(theta(x),x,x)+3*Pr*f(x)*diff(theta(x),x)+s*theta(x)=0;
params:={Pr=1,s=1,m=2};
bcs:={f(0)=0,D(f)(0)=0,theta(0)=1,D(f)(inf)=1,theta(inf)=0};
sys:=eval({ode1,ode2},params);
inf:=28:
res:=dsolve(sys union bcs,numeric);
plots:-odeplot(res,[[x,diff(f(x),x)],[x,theta(x)]],0..inf);
plots:-odeplot(res,[[x,f(x)],[x,theta(x)]],0..inf);

##What is the meaning of A and B in the discussion above in relation to the system of odes?

@rlewis Did you not just copy and paste all the lines I gave?

I just did that myself directly from my answer above. It ran without problems in Maple 2015.1.

Which Maple version are you using?

## Note: It seems that you are using a version behaving like Maple 12, see below. ###############

##########
I just tried in Maple 15. There the third ball is white (!!), thus cannot be seen!
You can change that by setting the color, either a common color as in color=blue or 3 different colors as shown in the revised lines below. There I also add the option axes=box. In Maple 2015 this is the default, in Maple 15 the default is axes = none.

restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50,color=[blue,red,green]); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
animate(display,['pointplot3d'(pts,symbol=solidsphere,symbolsize=50,color=[blue,red,green],axes=box),'pointplot3d'(ptsL,connect,thickness=3,color=black)],t=0..50,frames=100);

###################### Maple 12 #######
In Maple 12 I see the error message you report. Instead of the previous method of animation you can do:
restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50,color=[blue,red,green]); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
p:=tt->display(pointplot3d(eval(pts,t=tt),symbol=solidsphere,symbolsize=50,color=[blue,red,green],axes=box),pointplot3d(eval(ptsL,t=tt),connect,thickness=3,color=black));
animate(p,[t],t=0..50,frames=100);

@Kanellopoulos There are a lot of things to vary in these kinds of numerical problems, so I wouldn't dare conclude much else than for the settings we are dealing with in Carl's worksheet, results don't seem reliable.
I tried (using sysS as described above) to look at relative residuals. By that I mean the difference between the left and right hand sides of the pdes divided by the sum of their absolute values.
Since sysS is a set as I defined it, I replaced it by
sysS1:=convert(sysS,list);
so that the order is known.
Using the relative residuals meant replacing (lhs-rhs) by
((lhs-rhs)/(abs@rhs+abs@lhs))
in the definition of Residuals.
The results were rather discouraging. I found that in 8 pairs of residuals one of them (always the second: pde2) the relative residual had absolute value 1.

@Carl Love I tried your worksheet with the following change.
Using your sys I did
sysS:=solve(value(sys),{diff(h(x,t),t),diff(u(x,t),t)});

Then I used that instead of sys in the rest of the worksheet, i.e. in defining sol and in defining Residuals.

I got residuals over 0.1 in absolute value in 7 of the 18 (X,T) cases. One even over 3.

@tomleslie I tried making a new document in document mode and 2d input, just 3 lines. I copied this and pasted it into a fresh worksheet (worksheet mode and 1d input).
There I converted the pasted code into 1d-input. That worked, but there are no prompts or execution group delimiters and the area covered by the lines seems to be a document block.
This can be changed by using the menu item Format/Remove Document Block.
Unfortunately this meant that (in my case) 3 lines of
print(??); # input placeholder
also appeared.
#############
I just had a look at the discussion about converting a file made in document mode and 2d input into a file in worksheet mode and 1d input.
http://www.mapleprimes.com/questions/203080-How-To-Convert-Document-In-Document

I had occasion to do this the other day
http://www.mapleprimes.com/questions/204487-Too-Lengthy-Parametric-Equations

It would have been easier to do the conversion by exporting to .mpl as suggested by Alejandro Jakubi in the first reference.

@cuongtd The output from allvalues(solve....) is not a numerical result, but is the exact result.
As it turns out all 32 solutions (16 real and 16 imaginary) can be expressed in terms of the 16 roots of a polynomial of degree 16. Once you use evalf on that exact result numerical rootfinding is used, but since the problem just lies with the polynomial surely all roots are found.
To see the result before using evalf just change the colon in the first line below into a semicolon: It is huge!
#I'm assuming that eqs and eqs2 are defined as I did in my answer. Still using Digits=30.

resA:=allvalues([solve(eqs union eqs2)]): ##Sequence of length 16 each consisting of a list of two solutions
nops(resA[1]); # A pair of 2 solutions
nops([resA]); #16 pairs
nops(ListTools:-Flatten([resA])); #32 solutions
indets(resA[1],specfunc(RootOf)); #Notice the two arguments to RootOf: a polynomial of degree 16 and an index
ROF:=indets([resA],specfunc(RootOf)); #All of the RootOf's
nops(ROF); #16
pol:=map2(op,1,ROF); #The same polynomial for all
map2(op,2,ROF); #The 16 indices
#Now using evalf as I originally did.
res:=evalf(resA):
resR:=remove(has,[res],I); #Looking at the real results only
##Comment: It so happens that each two in the 16 pairs are either both real or both imaginary. I checked that.
map(nops,resR); #8 pairs
resRflat:=ListTools:-Flatten(resR);
nops(%); #16
##Graphical illustration: Animation showing one solution at a time.
##The 4 points in each picture are (ci,si), i=1..4.
plts:=seq(plots:-pointplot(subs(resRflat[j],[seq([cat(c,i),cat(s,i)],i=1..4)]),color=[red,blue,green,maroon],symbolsize=25),j=1..nops(resRflat)):
plots:-display(plts,insequence);
##Note: Recall that my (c4,s4) is the point (cos(t4),sin(t4)), where t4 = t+t1.